FindDistributionParameters
FindDistributionParameters[data,dist]
finds the parameter estimates for the distribution dist from data.
FindDistributionParameters[data,dist,{{p,p0},{q,q0},…}]
finds the parameters p, q, … with starting values p0, q0, ….
Details and Options
- FindDistributionParameters returns a list of replacement rules for the parameters in dist.
- The data must be a list of possible outcomes from the given distribution dist.
- The distribution dist can be any parametric univariate, multivariate, or derived distribution with unknown parameters.
- The following options can be given:
-
AccuracyGoal Automatic the accuracy sought ParameterEstimator "MaximumLikelihood" what parameter estimator to use PrecisionGoal Automatic the precision sought WorkingPrecision Automatic the precision used in internal computations - The following basic settings can be used for ParameterEstimator:
-
"MaximumLikelihood" maximize the log‐likelihood function "MethodOfMoments" match raw moments "MethodOfCentralMoments" match central moments "MethodOfCumulants" match cumulants "MethodOfFactorialMoments" match factorial moments - The maximum likelihood method attempts to maximize the log-likelihood function
, where 
are the distribution parameters and 
is the PDF of the distribution. - The method of moments solves
, 
, 
where 
is the 
sample moment and 
is the 
moment of the distribution with parameters 
. - Method-of-moment-based estimators may not satisfy all restrictions on parameters.
Examples
open all close allBasic Examples (3)
Obtain the maximum likelihood parameter estimates assuming a Laplace distribution:
FindDistributionParameters[{1.5, 2.8, 4.3, .5}, LaplaceDistribution[μ, σ]]Obtain the method of moments estimates:
FindDistributionParameters[{1.5, 2.8, 4.3, .5}, LaplaceDistribution[μ, σ], ParameterEstimator -> "MethodOfMoments"]Estimate parameters for a multivariate distribution:
data = RandomVariate[DirichletDistribution[{3, 4, 5}], 100];params = FindDistributionParameters[data, DirichletDistribution[{Subscript[α, 1], Subscript[α, 2], Subscript[α, 3]}]]Compare the difference between the original and estimated PDFs:
Plot3D[PDF[DirichletDistribution[{Subscript[α, 1], Subscript[α, 2], Subscript[α, 3]}] /. params, {x, y}] - PDF[DirichletDistribution[{3, 4, 5}], {x, y}], {x, 0, 1}, {y, 0, 1}, PlotRange -> All]Estimate parameters from quantity data:
data = RandomVariate[LaplaceDistribution[Quantity[100, "Grams"], Quantity[500, "Milligrams"]], 100]FindDistributionParameters[data, LaplaceDistribution[m, s]]Scope (15)
Basic Uses (5)
Estimate both parameters for a binomial distribution:
bdata = {23, 16, 24, 17, 16, 17, 18, 17, 19, 23, 19, 13, 19, 18, 22, 21, 19, 22, 19, 23};FindDistributionParameters[bdata, BinomialDistribution[n, p]]Estimate p, assuming n is known:
FindDistributionParameters[bdata, BinomialDistribution[50, p]]Estimate n, assuming p is known:
FindDistributionParameters[bdata, BinomialDistribution[n, 2 / 5]]Get the distribution with maximum likelihood parameter estimate for a particular family:
vals = RandomReal[ChiSquareDistribution[20], 1000];res = FindDistributionParameters[vals, ChiSquareDistribution[n]]Check goodness of fit by comparing a histogram of the data and the estimate's PDF:
Show[Histogram[vals, Automatic, "ProbabilityDensity"], Plot[PDF[ChiSquareDistribution[n] /. res, x], {x, 0, 50}]]Perform goodness-of-fit tests with null distribution from res:
DistributionFitTest[vals, ChiSquareDistribution[n] /. res, {"TestDataTable", All}]Perform tests correcting for estimation of the parameter:
DistributionFitTest[vals, ChiSquareDistribution[n], {"TestDataTable", All}]Estimate parameters by maximizing the log‐likelihood:
data = RandomVariate[GeometricDistribution[.4], 100];phat = FindDistributionParameters[data, GeometricDistribution[p]]Plot the log‐likelihood function to visually check that the solution is optimal:
Show[Plot[LogLikelihood[GeometricDistribution[p], data], {p, 0, 1}], Graphics[{PointSize[Large], Red, Point[{p, LogLikelihood[GeometricDistribution[p], data]} /. phat]}]]Visualize a log‐likelihood surface to find rough values for the parameters:
data = BlockRandom[SeedRandom[100];RandomVariate[WeibullDistribution[4, 5], 50]];cp = ContourPlot[LogLikelihood[WeibullDistribution[α, β], data], {α, 2, 6}, {β, 3, 8}]Supply those rough values as starting values for the estimation:
params = FindDistributionParameters[data, WeibullDistribution[α, β], {{α, 3}, {β, 5.5}}]Mark the optimal point on the contour plot:
Show[cp, Graphics[{PointSize[Large], Red, Point[{α, β} /. params]}]]Estimate the normal approximation of Poisson data:
data = RandomVariate[PoissonDistribution[20], 100];FindDistributionParameters[data, NormalDistribution[μ, σ]]FindDistributionParameters[data, NormalDistribution[μ, σ], WorkingPrecision -> 20]Univariate Parametric Distributions (2)
Estimate parameters for a continuous distribution:
data = {0.8, 0.98, 0.66, 0.09, 0.03, 0.16, 0.88, 0.32, 0.56, 0.3, 0.82, 0.34, 0.93, 0.33, 0.25, 0.74, 0.3, 0.74, 0.16, 0.02};FindDistributionParameters[data, BetaDistribution[α, β]]Estimate parameters for a discrete distribution:
data = {16, 1, 2, 5, 3, 5, 1, 6, 4, 14, 2, 9, 1, 1, 2, 2, 8, 1, 3, 1};param = FindDistributionParameters[data, LogSeriesDistribution[θ]]Compare the fitted and empirical CDFs:
Plot[{CDF[LogSeriesDistribution[θ] /. param, x], CDF[EmpiricalDistribution[data], x]}, {x, 0, 20}, Exclusions -> None]Multivariate Parametric Distributions (2)
Estimate parameters for a discrete multivariate distribution:
mdata = RandomVariate[MultinomialDistribution[20, {1 / 3, 1 / 6, 1 / 10, 2 / 5}], 1000];FindDistributionParameters[mdata, MultinomialDistribution[n, {Subscript[p, 1], Subscript[p, 2], Subscript[p, 3], Subscript[p, 4]}]]Estimate parameters for a continuous multivariate distribution:
data = RandomVariate[BinormalDistribution[{1, 2}, {1 / 3, 4}, 3 / 4], 1000];params = FindDistributionParameters[data, BinormalDistribution[{Subscript[μ, 1], Subscript[μ, 2]}, {Subscript[σ, 1], Subscript[σ, 2]}, ρ]]Visualize the density functions for the marginal distributions:
Table[Plot[PDF[NormalDistribution[Subscript[μ, i], Subscript[σ, i]] /. params, x], {x, Min[data[[All, i]]], Max[data[[All, i]]]}, PlotRange -> All], {i, 2}]//RowObtain the covariance matrix from the formula:
{{Subscript[σ, 1]^2, Subscript[σ, 1] Subscript[σ, 2] ρ}, {Subscript[σ, 1] Subscript[σ, 2] ρ, Subscript[σ, 2]^2}} /. params//MatrixFormDerived Distributions (6)
Estimate parameters for a truncated normal:
SeedRandom[1];
data = RandomVariate[TruncatedDistribution[{-1, 2}, NormalDistribution[1, 1 / 3]], 100];FindDistributionParameters[data, TruncatedDistribution[{min, max}, NormalDistribution[μ, σ]]]Estimate parameters for a constructed distribution:
dist = ProbabilityDistribution[λ Exp[λ(1 - x)], {x, 1, Infinity}, Assumptions -> λ > 0]data = {1.4, 5.1, 1.7, 1.6, 1.1, 3.9, 2.2, 1.3, 2., 1.5};param = FindDistributionParameters[data, dist]Plot[LogLikelihood[dist, data], {λ, .1, 2}, Epilog -> {PointSize[Medium], Red, Point[{λ, LogLikelihood[dist, data]} /. param]}]Estimate parameters for a product distribution:
data = RandomVariate[ProductDistribution[BetaDistribution[2, 3], LaplaceDistribution[10, 1]], 100];FindDistributionParameters[data, ProductDistribution[BetaDistribution[α, β], LaplaceDistribution[μ, σ]]]Estimate parameters for a copula distribution:
data = RandomVariate[CopulaDistribution[{"Frank", 1 / 2}, {GammaDistribution[2, 3], ChiSquareDistribution[10]}], 100];FindDistributionParameters[data, CopulaDistribution[{"Frank", c}, {GammaDistribution[α, β], ChiSquareDistribution[n]}]]Estimate parameters for a component mixture:
data = RandomVariate[MixtureDistribution[{2 / 3, 1 / 3}, {GammaDistribution[2, 3], NormalDistribution[1, 1 / 4]}], 100];𝒟1 = MixtureDistribution[{p, 1 - p}, {GammaDistribution[α, β], NormalDistribution[μ, σ]}];param1 = FindDistributionParameters[data, 𝒟1]//QuietEstimate the mixture probabilities assuming the component distributions are known:
𝒟2 = MixtureDistribution[{p, 1 - p}, {GammaDistribution[2, 3], NormalDistribution[1, 1 / 4]}];param2 = FindDistributionParameters[data, 𝒟2]Visualize the two estimates against the data:
Show[Histogram[data, {0, 20, .5}, "PDF"], Plot[{PDF[𝒟1 /. param1, x], PDF[𝒟2 /. param2, x]}, {x, 0, 20}, PlotRange -> All, PlotLegends -> {"𝒟1", "𝒟2"}]]Estimate parameters for a distribution in specified units:
data = RandomVariate[QuantityDistribution[ParetoDistribution[2, 3], "Mega" * "USDollars"], 100]FindDistributionParameters[data, QuantityDistribution[ParetoDistribution[minSalary, α], "USDollars"]]Options (4)
ParameterEstimator (3)
Estimate parameters by matching cumulants:
data = RandomVariate[ParetoDistribution[2, 8, 5], 1000];FindDistributionParameters[data, ParetoDistribution[α, β, μ], ParameterEstimator -> "MethodOfCumulants"]Other moment‐based methods typically give similar results:
FindDistributionParameters[data, ParetoDistribution[α, β, μ], ParameterEstimator -> "MethodOfCentralMoments"]FindDistributionParameters[data, ParetoDistribution[α, β, μ], ParameterEstimator -> "MethodOfFactorialMoments"]Estimate parameters based on default moments:
data = RandomVariate[NormalDistribution[2, 3], 1000];FindDistributionParameters[data, NormalDistribution[μ, σ], ParameterEstimator -> "MethodOfMoments"]Estimate parameters from the first and fourth moments:
FindDistributionParameters[data, NormalDistribution[μ, σ], ParameterEstimator -> {"MethodOfMoments", "MomentOrders" -> {1, 4}}]Obtain the maximum likelihood estimates using the default method:
data = RandomVariate[BetaDistribution[2, 3], 1000];FindDistributionParameters[data, BetaDistribution[α, β]]Use FindMaximum to obtain the estimates:
FindDistributionParameters[data, BetaDistribution[α, β], ParameterEstimator -> {"MaximumLikelihood", Method -> "FindMaximum"}]Use EvaluationMonitor to extract the points sampled:
{params, {points}} = Reap[FindDistributionParameters[data, BetaDistribution[α, β], ParameterEstimator -> {"MaximumLikelihood", Method -> {"FindMaximum", EvaluationMonitor :> Sow[{α, β}]}}]];Visualize the sequences of sampled 
and 
values:
ListLogPlot[Transpose[points], PlotLegends -> {"α", "β"}]WorkingPrecision (1)
Use machine precision for continuous parameters by default:
data = RandomVariate[GumbelDistribution[2, 5], 100, WorkingPrecision -> 50];FindDistributionParameters[data, GumbelDistribution[α, β]]Obtain a higher-precision result:
FindDistributionParameters[data, GumbelDistribution[α, β], WorkingPrecision -> 25]Applications (17)
Use One Parameter Estimator to Get Starting Values for Another (1)
data = {11.3, 19.2, 8.7, 17.5, 7., 17.1, 8.5, 9.8, 15.2, 10.4, 17.9, 20.3, 6.5, 7.9, 11.8, 8.5, 12.8, 14.2, 6.6, 9.};Get the method of moments estimate:
mm = FindDistributionParameters[data, ChiSquareDistribution[n], ParameterEstimator -> "MethodOfMoments"]Use the method of moments estimate as the starting value for ml estimation:
FindDistributionParameters[data, ChiSquareDistribution[n], mm /. Rule -> List]Obtain ml estimates for a gamma distribution:
mle = FindDistributionParameters[data, GammaDistribution[α, β]]Use those as starting values for the method of moments:
FindDistributionParameters[data, GammaDistribution[α, β], mle /. Rule -> List, ParameterEstimator -> "MethodOfMoments"]Obtain Starting Values for Another Estimation (1)
Estimate Laplace parameters for data from an ExponentialPowerDistribution:
data = RandomVariate[ExponentialPowerDistribution[1.4, 2, 1 / 3], 1000];lparams = FindDistributionParameters[data, LaplaceDistribution[μ, σ]]Use the Laplace estimate as a starting point for estimating exponential power parameters:
exppow = FindDistributionParameters[data, ExponentialPowerDistribution[k, μ, σ], Join[{{k, 1}}, lparams /. Rule -> List]]Compare the data with the Laplace and exponential power estimates:
Show[Histogram[data, Automatic, "PDF"], Plot[{PDF[LaplaceDistribution[μ, σ] /. lparams, x], PDF[ExponentialPowerDistribution[k, μ, σ] /. exppow, x]}, {x, 0, 5}, PlotRange -> All, PlotLegends -> {LaplaceDistribution, ExponentialPowerDistribution}]]Parameter Estimation of Similarly Shaped Distributions (1)
Model lognormal distributed data with a gamma distribution:
lnorm = LogNormalDistribution[2, 0.3];sample = RandomVariate[lnorm, 10 ^ 4];edist = FindDistributionParameters[sample, GammaDistribution[α, β]]Compare the distributions of the simulation and estimated distributions:
Plot[{PDF[lnorm, x], PDF[GammaDistribution[α, β] /. edist, x]}, {x, 0, 25}, PlotLegends -> {LogNormalDistribution, GammaDistribution}]Accident Claims (1)
The number of accident claims per policy per year from an insurance company:
AccidentsPerPolicy = {ConstantArray[0, {81714}], ConstantArray[1, {11306}], ConstantArray[2, {1618}], ConstantArray[3, {250}], ConstantArray[4, {40}], ConstantArray[5, {7}]}//Flatten;Estimate the parameter 
for a logarithmic series distribution for policy claims shifted by 1:
edist = FindDistributionParameters[AccidentsPerPolicy + 1, LogSeriesDistribution[θ]]See that the estimate gives a maximal result:
Show[Plot[LogLikelihood[LogSeriesDistribution[θ], AccidentsPerPolicy + 1], {θ, 0.2, .3}], Graphics[{PointSize[Medium], Red, Point[{θ, LogLikelihood[LogSeriesDistribution[θ], AccidentsPerPolicy + 1]} /. edist]}]]Word Lengths in Different Languages (1)
Get word length data for several languages:
languages = {"Dutch", "English", "French", "German", "Hebrew", "Polish", "Russian", "Spanish", "Swedish"};worddata = Table[StringLength /@ DictionaryLookup[{l, All}], {l, languages}];Model the word lengths for each language as binomially distributed with 
:
binom = Table[FindDistributionParameters[i, BinomialDistribution[40, p]], {i, worddata}]Compare the actual and estimated distributions:
Partition[Table[Show[Histogram[worddata[[i]], {Range[25] - 1 / 2}, "ProbabilityDensity", PlotLabel -> languages[[i]]], DiscretePlot[PDF[BinomialDistribution[40, p /. binom[[i]]], x], {x, 0, 25}, PlotRange -> All, PlotStyle -> PointSize[.02]]], {i, Length[languages]}], 3]//GridBootstrap the distribution of p values based on these 9 results:
bootstrap = Mean[RandomChoice[p /. binom, {9, 10000}]];Histogram[bootstrap, Automatic, "ProbabilityDensity"]Estimate the expected value of p and a standard deviation for the estimate:
FindDistributionParameters[bootstrap, NormalDistribution[p, Subscript[σ, p]]]Text Frequency (1)
The word count in a text follows a Zipf distribution:
text = ExampleData[{"Text", "OriginOfSpecies"}, "Words"];wordCount = Tally[text][[All, 2]];Fit a ZipfDistribution to the word frequency data:
rhohat = FindDistributionParameters[wordCount, ZipfDistribution[ρ]]Fit a truncated ZipfDistribution to counts at most 50 using rhohat as a starting value:
rhohat50 = FindDistributionParameters[Select[wordCount, # ≤ 50&], ZipfDistribution[50, ρ], rhohat /. Rule -> List]Visualize the CDFs up to the truncation value:
Plot[{CDF[d1 = ZipfDistribution[ρ] /. rhohat, x], CDF[d2 = ZipfDistribution[50, ρ] /. rhohat50, x]}, {x, 0, 50}, PlotLegends -> {d1, d2}]Estimate the proportion of the original data not included in the truncated model:
Probability[x > 50, xZipfDistribution[ρ]]% /. rhohatEarthquake Magnitudes (1)
Find estimates for a multimodal MixtureDistribution model:
ExampleData[{"Statistics", "USEarthquakes"}, "ColumnDescriptions"]magnitudes = Select[ExampleData[{"Statistics", "USEarthquakes"}], #[[1]] ≥ 1935&][[All, 7]];The magnitudes of earthquakes in the United States in the selected years have two modes:
h = Histogram[magnitudes, 20, "ProbabilityDensity"]Fit distribution from possible mixtures of one NormalDistribution with another:
dist = MixtureDistribution[{p, 1 - p}, {NormalDistribution[a, b], NormalDistribution[c, d]}];params = FindDistributionParameters[magnitudes, dist]//QuietExtract the means of the components:
means = {a, c} /. paramsThe components' means are far enough apart that they are still the modes:
Plot[PDF[dist /. params, x], {x, 2, 8}, PlotRange -> All, Epilog -> {Red, PointSize[Medium], Point[Table[{i, PDF[dist /. params, i]}, {i, means}]]}]Wind Speed Analysis (1)
Model monthly maximum wind speeds in Boston:
maxWinds = WeatherData["Boston", "MaxWindSpeed", {{1950, 1, 1}, {2009, 12, 31}, "Month"}, "Value"]//QuantityMagnitude;Fit the data to a RayleighDistribution:
rparam = FindDistributionParameters[maxWinds, RayleighDistribution[θ]]evparams = FindDistributionParameters[maxWinds, ExtremeValueDistribution[a, b]]Compare the empirical and fitted quantiles to see where the models deviate from the data:
Row[Table[QuantilePlot[maxWinds, i, PlotLabel -> Head[i]], {i, {RayleighDistribution[θ] /. rparam, ExtremeValueDistribution[a, b] /. evparams}}]]Distribution of Incomes (1)
Model incomes at a large state university:
ExampleData[{"Statistics", "UniversitySalaries"}, "ColumnDescriptions"]universityA = Cases[ExampleData[{"Statistics", "UniversitySalaries"}], {dept_, perc_, salary_, campus : "A"} :> {perc, salary}];salaries = Cases[universityA, {p_ ? Positive, s_ ? Positive} :> s / p];Assume the salaries are Dagum distributed:
dparams = FindDistributionParameters[salaries, DagumDistribution[p, a, b]]Assume they follow a more general Pareto distribution:
pparams = FindDistributionParameters[salaries, ParetoDistribution[k, α, γ, μ], ParameterEstimator -> {"MaximumLikelihood", "Method" -> "SimulatedAnnealing"}]Compare the subtle differences in the estimated distributions:
Show[Histogram[salaries, {0, 200000, 10000}, "PDF"], Plot[{PDF[DagumDistribution[p, a, b] /. dparams, x], PDF[ParetoDistribution[k, α, γ, μ] /. pparams, x]}, {x, 0, 200000}, PlotLegends -> {DagumDistribution, ParetoDistribution}]]Market Change in Stock Values (1)
Use a beta distribution to model the proportion of Dow Jones Industrial stocks that increase in value on a given day:
dji = EntityList[EntityClass["Financial", "DowJonesIndustrialAverageComponents"]]Find daily change for Dow Jones Industrial stocks:
Clear[djichange];djichange = DeleteMissing @Table[FinancialData[tick, "FractionalChange", {{2009, 3, 1}, {2013, 3, 1}, "Day"}, "Value"], {tick, dji}];Number of days for each financial entity:
days = First[djichange]["PathLength"]Extract values from time series for each entity and normalize numeric quantities:
djichangevals = Normal /@ Values /@ djichange;Check if each entity has the same length of data:
ArrayQ[djichangevals]MatrixPlot[djichangevals, AspectRatio -> 1, FrameLabel -> {"Company", "Day"}]Calculate the daily ratio of companies with an increase in value:
ratios = Table[Count[djichangevals[[All, i]], _ ? Positive] / Length[dji], {i, days}];Find parameter estimates, excluding days with zero or all companies having an increase in value:
posratios = Select[ratios, 0 < # < 1&];
bparams = FindDistributionParameters[posratios, BetaDistribution[α, β]]Visualize the likelihood contours and mark the optimal point:
ContourPlot[LogLikelihood[BetaDistribution[α, β], posratios], {α, 1, 1.5}, {β, 1, 1.4}, Epilog -> {Red, PointSize[Large], Point[{α, β} /. bparams]}]Automobile Fuel Efficiency (1)
The average city and highway mileage for midsize cars follows a binormal distribution:
midsizeMPG2009 = {{16, 22}, {18, 26}, {17, 25}, {16, 23}, {16, 23}, {14, 19}, {13, 19}, {10, 14}, {10, 17}, {18, 28}, {18, 27}, {17, 25}, {17, 25}, {17, 26}, {17, 26}, {16, 25}, {17, 25}, {15, 22}, {15, 23}, {11, 17}, {11, 17}, {17, 28}, {16, 24}, {16, 25}, {16, 25}, {17, 26}, {18, 26}, {17, 26}, {17, 25}, {17, 26}, {13, 19}, {15, 24}, {17, 26}, {15, 22}, {22, 33}, {22, 30}, {18, 29}, {17, 26}, {26, 34}, {21, 30}, {13, 20}, {19, 27}, {16, 27}, {21, 30}, {13, 20}, {19, 27}, {16, 27}, {24, 30}, {23, 27}, {23, 29}, {21, 25}, {19, 27}, {10, 15}, {9, 16}, {17, 25}, {20, 29}, {20, 28}, {18, 26}, {24, 33}, {25, 33}, {16, 25}, {15, 23}, {22, 32}, {22, 32}, {20, 28}, {23, 30}, {24, 32}, {19, 27}, {19, 26}, {18, 25}, {17, 24}, {16, 24}, {16, 23}, {16, 24}, {16, 23}, {20, 22}, {17, 25}, {20, 29}, {20, 28}, {18, 26}, {17, 24}, {18, 28}, {20, 29}, {21, 30}, {17, 25}, {18, 25}, {23, 32}, {17, 24}, {16, 22}, {15, 22}, {13, 19}, {13, 20}, {20, 27}, {16, 25}, {23, 32}, {23, 31}, {18, 27}, {19, 26}, {35, 33}, {19, 26}, {24, 30}, {21, 29}, {26, 31}, {27, 33}, {24, 32}, {11, 18}, {24, 30}, {24, 32}, {22, 33}, {17, 26}, {26, 34}, {21, 31}, {21, 31}, {19, 28}, {33, 34}, {48, 45}, {19, 29}, {15, 23}, {15, 22}, {16, 25}};Assume city and highway miles per gallon are normally distributed and correlated:
bparams = FindDistributionParameters[midsizeMPG2009, BinormalDistribution[{μ1, μ2}, {σ1, σ2}, ρ]]Extract the estimated average city and highway mileages:
{μ1, μ2} /. bparamsExtract the estimated correlation between city and highway mileages:
ρ /. bparamsVisualize the joint density on a logarithmic scale with the mean mileage marked with a blue point:
Show[ContourPlot[Log[PDF[BinormalDistribution[{μ1, μ2}, {σ1, σ2}, ρ] /. bparams, {x, y}]], {x, 0, 50}, {y, 0, 50}, Contours -> 50], Graphics[{PointSize[Medium], Red, Point[{μ1, μ2} /. bparams]}]]Earthquake Waiting Times (1)
The data contains waiting times in days between serious (magnitude at least 7.5 or over 1000 fatalities) earthquakes worldwide, recorded from 12/16/1902 to 3/4/1977:
earthquakesWaitingTimes = ExampleData[{"Statistics", "EarthquakeWaitingTimes"}];Model waiting times by an ExponentialDistribution:
param = FindDistributionParameters[earthquakesWaitingTimes, ExponentialDistribution[λ]]Estimate the average and median number of days between major earthquakes:
{Mean[ExponentialDistribution[λ]], Median[ExponentialDistribution[λ]]} /. paramEarthquake Frequency (1)
The number of earthquakes per year can be modeled by SinghMaddalaDistribution:
ExampleData[{"Statistics", "USEarthquakes"}, "Description"]earthquakes = Tally[Select[ExampleData[{"Statistics", "USEarthquakes"}], #[[1]] ≥ 1935&][[All, 1]]][[All, 2]]Fit the distribution to the data:
params = FindDistributionParameters[earthquakes, SinghMaddalaDistribution[q, a, b]]Compute the maximized log‐likelihood:
ll = LogLikelihood[SinghMaddalaDistribution[q, a, b], earthquakes] /. paramsVisualize the log‐likelihood profiles near the optimal parameter values:
Table[With[{var = params[[i, 1]], est = params[[i, 2]]}, Plot[LogLikelihood[SinghMaddalaDistribution[q, a, b] /. Drop[params, {i}], earthquakes], {var, .8est, 1.2est}, Epilog -> {PointSize[Medium], Point[{est, ll}]}, PlotRange -> {Floor[ll], Ceiling[ll]}, PlotLabel -> params[[i]]]], {i, 3}]//RowTime between Geyser Eruptions (1)
Mixtures can be used to model multimodal data:
ExampleData[{"Statistics", "OldFaithful"}, "ColumnDescriptions"]waiting = ExampleData[{"Statistics", "OldFaithful"}][[All, 2]];A histogram of waiting times for eruptions of the Old Faithful geyser exhibits two modes:
h = Histogram[waiting, 25, "ProbabilityDensity"]Fit a mixture of gamma and normal distributions to the data:
mdist = MixtureDistribution[{p, 1 - p}, {GammaDistribution[α, β], NormalDistribution[μ, σ]}];mparams = FindDistributionParameters[waiting, mdist, {{p, 1 / 3}, {α, 90}, {β, .5}, {μ, 80}, {σ, 5}}]Compare the histogram to the PDF of the estimated distribution:
Show[h, Plot[PDF[mdist /. mparams, x], {x, 0, 100}, PlotStyle -> Thick]]Stock Price Distribution (1)
Lognormal distribution can be used to model stock prices:
googleStock = QuantityMagnitude@FinancialData["GOOGL", {{2006, 3, 10}, {2018, 3, 15}, "Day"}, "Value"];Fit the distribution to the data:
params = FindDistributionParameters[googleStock, LogNormalDistribution[μ, σ]]Visualize the profile likelihoods, fixing one parameter at the fitted value:
Block[{ml = LogLikelihood[LogNormalDistribution[μ, σ] /. params, googleStock], llh1 = LogLikelihood[LogNormalDistribution[mu, σ /. params], googleStock],
llh2 = LogLikelihood[LogNormalDistribution[μ /. params, s], googleStock]
}, GraphicsRow[{Plot[llh1, {mu, 5.5, 6.5}, Epilog -> {PointSize[Large], Point[{μ /. params, ml}]}, AxesLabel -> {μ, None}, PlotLabel -> First[params]], Plot[llh2, {s, .1, 1}, Epilog -> {PointSize[Large], Point[{σ /. params, ml}]}, AxesLabel -> {σ, None}, PlotLabel -> Last[params]]}]]Water Flow Rates (1)
Consider the annual minimum daily flows given in cubic meters per second for the Mahanadi river:
minFlow = ExampleData[{"Statistics", "MahanadiRiverFlow"}];Model the annual minimum mean daily flows as a MinStableDistribution:
params = FindDistributionParameters[minFlow, MinStableDistribution[a, b, c]]Simulate annual minimum mean daily flows for the next 30 years:
ListPlot[RandomVariate[MinStableDistribution[a, b, c] /. params, 30], Filling -> 0, AxesOrigin -> {0, 0}]Population Sizes (1)
Use a Pareto distribution to model Australian city population sizes:
cities = CityData[{All, "Australia"}];populations = QuantityMagnitude[CityData[#, "Population"]& /@ cities];params = FindDistributionParameters[populations, ParetoDistribution[k, α, γ, μ]]Get the probability that a city has a population at least 10000 under a Pareto distribution:
Probability[x ≥ 10000, xParetoDistribution[k, α, γ, μ]]Compute the probability given the parameter estimates:
% /. paramsCompute the probability based on the original data:
Probability[x ≥ 10000, xpopulations]//NProperties & Relations (8)
FindDistributionParameters gives estimates as replacement rules:
FindDistributionParameters[{5, 8, 3, 4, 9, 6}, PoissonDistribution[μ]]EstimatedDistribution gives a distribution with parameter estimates inserted:
EstimatedDistribution[{5, 8, 3, 4, 9, 6}, PoissonDistribution[μ]]FindProcessParameters returns a list of parameter estimates for a random process:
data = RandomFunction[BernoulliProcess[0.4], {0, 10 ^ 5}];FindProcessParameters[data, BernoulliProcess[p]]FindDistributionParameters returns a list of parameter estimates for a distribution:
data = RandomVariate[BernoulliDistribution[0.4], 10 ^ 5];FindDistributionParameters[data, BernoulliDistribution[p]]Estimate distribution parameters by maximum likelihood:
data = RandomVariate[ExponentialDistribution[3], 1000];FindDistributionParameters[data, ExponentialDistribution[λ]]Use DistributionFitTest to test quality of the fit:
ℋ = DistributionFitTest[data, ExponentialDistribution[λ], "HypothesisTestData"]Extract the fitted distribution parameter:
ℋ["FittedDistributionParameters"]Obtain a table of relevant test statistics and p‐values:
ℋ[{"TestDataTable", All}]Estimate parameters in a parametric distribution:
data = RandomVariate[NormalDistribution[0, 1], 50];pars = FindDistributionParameters[data, NormalDistribution[μ, σ]]Get a nonparametric kernel density estimate using SmoothKernelDistribution:
skdist = SmoothKernelDistribution[data]Compare the PDFs for the nonparametric and parametric distributions:
Plot[{PDF[skdist, x], PDF[NormalDistribution[μ, σ] /. pars, x]}, {x, -3, 3}, PlotLegends -> {SmoothKernelDistribution, NormalDistribution}]Visualize the nonparametric density using SmoothHistogram:
SmoothHistogram[data, PlotStyle -> Thick]Get a maximum likelihood estimate of parameters:
data = RandomVariate[GammaDistribution[4, 10], 25];edist = FindDistributionParameters[data, GammaDistribution[α, β]]Compute the likelihood using Likelihood:
Likelihood[GammaDistribution[α, β] /. edist, data]Compute the log‐likelihood using LogLikelihood:
LogLikelihood[GammaDistribution[α, β] /. edist, data]Estimate parameters by matching raw moments:
data = RandomVariate[BetaDistribution[2, 5], 100];edist = FindDistributionParameters[data, BetaDistribution[α, β], ParameterEstimator -> "MethodOfMoments"]Compute raw moments from the data using Moment:
{Moment[data, 1], Moment[data, 2]}Compute the same moments from the beta distribution for the estimated parameters:
{Moment[BetaDistribution[α, β], 1], Moment[BetaDistribution[α, β], 2]} /. edistEstimate parameters for a Weibull distribution:
data = RandomVariate[WeibullDistribution[3, 15], 100];edist = FindDistributionParameters[data, WeibullDistribution[α, β]]Use QuantilePlot to visualize the empirical quantiles versus the theoretical quantiles:
QuantilePlot[data, WeibullDistribution[α, β] /. edist]Obtain the same visualization when the estimation is done within QuantilePlot:
QuantilePlot[data, WeibullDistribution[α, β]]FindDistributionParameters ignores time stamps in TimeSeries and EventSeries:
ts = TemporalData[TimeSeries, {{{1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1,
0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0,
0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0,
1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1}}, {{0, 99, 1}}, 1, {"Continuous", 1}, {"Discrete", 1}, 1,
{ResamplingMethod -> {"Interpolation", InterpolationOrder -> 1}}}, False, 10.1];FindDistributionParameters[ts, BernoulliDistribution[p]]FindDistributionParameters[ts["Values"], BernoulliDistribution[p]]For TemporalData, all the path structure is ignored:
td = TemporalData[«4»];FindDistributionParameters[td, NormalDistribution[a, b]]FindDistributionParameters[td["ValueList"]//Flatten, NormalDistribution[a, b]]Possible Issues (3)
Solutions of method-of-moment equations can give parameters that are not valid:
FindDistributionParameters[{10, 12, 8, 15}, BinomialDistribution[n, p], ParameterEstimator -> "MethodOfMoments"]DistributionParameterAssumptions[BinomialDistribution[n, p]] /. %For a continuous distribution:
data = RandomVariate[NormalDistribution[-1, .1], 10 ^ 3];FindDistributionParameters[data, NormalDistribution[m, m], ParameterEstimator -> "MethodOfMoments"]DistributionParameterQ[NormalDistribution[m, m] /. %]
Good starting values may be needed to obtain a good solution:
data = {4.1, 5.8, 4.3, 5.6, 3.8, 2.5, 3.1, 5.3, 5.4, 3.3};FindDistributionParameters[data, RiceDistribution[.01, α, β], ParameterEstimator -> "MethodOfMoments"]
FindDistributionParameters[data, RiceDistribution[.01, α, β], {{α, 4}, {β, .02}}, ParameterEstimator -> "MethodOfMoments"]Good starting values may result in quicker results:
data = RandomVariate[DirichletDistribution[Range[25]], 10 ^ 4];FindDistributionParameters[data, DirichletDistribution[Array[g, 25]]]//AbsoluteTimingFindDistributionParameters[data, DirichletDistribution[Array[g, 25]], Transpose[{Array[g, 25], Range[25]}]]//AbsoluteTimingText
Wolfram Research (2010), FindDistributionParameters, Wolfram Language function, https://reference.wolfram.com/language/ref/FindDistributionParameters.html.
CMS
Wolfram Language. 2010. "FindDistributionParameters." Wolfram Language & System Documentation Center. Wolfram Research. https://reference.wolfram.com/language/ref/FindDistributionParameters.html.
APA
Wolfram Language. (2010). FindDistributionParameters. Wolfram Language & System Documentation Center. Retrieved from https://reference.wolfram.com/language/ref/FindDistributionParameters.html