3
$\begingroup$

Say I have data that I want to fit to a distribution while applying assumptions on one or more of the parameters. For example, for ParetoDistribution[k, α, μ]

SeedRandom[567];
vals = RandomVariate[ParetoDistribution[22 10^6, .5, 5 10^6], 15];

For example, I might like to use FindDistributionParameters with the constraint .4 < α < .5 while leaving k and μ free to be selected over their domain.

Without this assumption

FindDistributionParameters[vals, ParetoDistribution[k, α, μ]]
{k -> 25.1961, α -> 0.0714327, μ -> 1.01492*10^7}

with α far from the constrained range I am seeking.

Is there a generalised method to inject parameter constraints into FindDistributionParameters?

I tried playing with ParameterMixtureDistribution using UniformDistribution on the parameter range.

FindDistributionParameters[vals,
 ParameterMixtureDistribution[
  ParetoDistribution[k, α, μ],
  {
   k \[Distributed] UniformDistribution[{0, ∞}],
   α \[Distributed] UniformDistribution[{0, ∞}],
   μ \[Distributed] UniformDistribution[{0, Min[vals]}]
   },
  Assumptions -> .4 < α < .5
  ]
 ]

However, I have not gotten too far with this.

The use case is mapping the fit of data to a collection of distributions each with its own set of parameter constraints based on the observation of the data. This is required for some distributions to obtain a fit with parameters such that mean and standard deviation are defined. For example in version 11.3 Statistics`Library`ParetoPickandsDistribution[μ, σ, ξ] a fit with ξ > 1 is a valid distribution but mean and standard deviation is undefined.

$\endgroup$
  • 1
    $\begingroup$ I know your question is about the mechanics of performing the estimation with the restrictions, but I just have to ask: If your data says that $\alpha$ is much closer to 0.07, how are you justifying a restriction between 0.4 and 0.5? Your start on this with the random distributions of the parameters is almost Bayesian. Why not go full Bayesian where your "restrictions" are your prior beliefs (in terms of probability distributions) and are modified by your data? $\endgroup$ – JimB Sep 6 '18 at 1:54
  • $\begingroup$ Sorry, one more sermon: probability distributions without means and standard deviations are not second class distributions. The parameters can still be estimated with associated estimates of precision. $\endgroup$ – JimB Sep 6 '18 at 2:14
  • $\begingroup$ Using starting values for the parameters comes closer. FindDistributionParameters[vals, ParetoDistribution[k, α, μ], {{k, 1*^7}, {α, 0.45}, {μ, 1*^6}}] evaluates to {k -> 7.20165*10^6, α -> 0.380515, μ -> 1.01492*10^7} $\endgroup$ – Bob Hanlon Sep 6 '18 at 2:37
  • $\begingroup$ @JimB No chance in increasing the sample size as I am working with real data where 15 points is considered a larger sample due to availability of the data. However, by applying knowledge on the behavior of the source of data and the impact of certain distribution parameters on the PDF of the distribution then assumptions can be introduced to produce a "reasonable" fit. $\endgroup$ – Edmund Sep 7 '18 at 12:21
  • $\begingroup$ @BobHanlon I initially considered this approach before posting but it is less than ideal because it cannot constrain a particular set of parameters and an initial value must be applied for all parameters. I only want to target certain parameters and apply specific constraints to each one. $\endgroup$ – Edmund Sep 7 '18 at 12:25
1
$\begingroup$
amin = .4;
amax = .5;
EstimatedDistribution[vals, ParetoDistribution[k, #, μ]] &@ 
 ArgMax[{LogLikelihood[EstimatedDistribution[vals, ParetoDistribution[k, a, μ]] , vals], 
   amin <= a <= amax}, a]

ParetoDistribution[8.61328*10^6, 0.4, 1.01492*10^7]

Append[FindDistributionParameters[vals, ParetoDistribution[k, #, μ]], α -> #] & @ 
 ArgMax[{LogLikelihood[EstimatedDistribution[vals, 
     ParetoDistribution[k, a, μ]] , vals], amin <= a <= amax}  , a]

{k -> 8.61328*10^6, μ -> 1.01492*10^7, α -> 0.4}

$\endgroup$
  • $\begingroup$ I like the that this approach essentially bootstraps the "MaximumLikelihood" method that FindDistributionParameters uses as its default with ArgMax to inject the parameter constraints. Also more than one parameter can be constrained with this method. (+1). $\endgroup$ – Edmund Sep 7 '18 at 12:19
  • $\begingroup$ @Edmund More sermons: There is no need to attempt to create new statistical theory. As Scotty says: "Cap'n! You can't change the laws of statistics!". A bootstrap that will give you a confidence interval is what is needed (as opposed to a point estimate) to judge whether any restrictions or re-examining the data collection process is necessary. $\endgroup$ – JimB Sep 7 '18 at 13:28
1
$\begingroup$

The answer is not to restrict the parameter(s). The answer is to take a larger sample size. The Type 2 Pareto Distribution you've chosen has 3 parameters (with values of very different scales) and expecting 15 random samples to do a good job is just not a reasonable expectation.

The estimated standard errors (which FindDistributionParameters doesn't provide) would likely be so large to cover the restricted range.

Using a sample size of 150 provides much better estimates.

SeedRandom[567];
vals = RandomVariate[ParetoDistribution[22 10^6, 0.5, 5 10^6], 150];
FindDistributionParameters[vals, ParetoDistribution[k, α, μ]]
(* {k -> 2.23607*10^7, α -> 0.51351, μ -> 5.26508*10^6} *)

Edit

I need to eat a bit of crow. While the sample size of 15 is a bit low, I now think the numerical instability induced by the wide range of scale for the parameters is more likely the dominating issue (although there's still no reason to restrict the parameters).

By scaling the parameters, more appropriate maximum likelihood estimators are found:

SeedRandom[567];
vals = RandomVariate[ParetoDistribution[22 10^6, 0.5, 5 10^6], 15];
sol = {k 10^6, α, μ 10^6} /.
  FindDistributionParameters[vals, ParetoDistribution[k 10^6, α, μ 10^6]]
(* {7.20165*10^6, 0.380515, 1.01492*10^7} *)

With a sample size of 150 the estimates are much better:

SeedRandom[567];
vals = RandomVariate[ParetoDistribution[22 10^6, 0.5, 5 10^6], 150];
sol = {k 10^6, α, μ 10^6} /.
  FindDistributionParameters[vals, ParetoDistribution[k 10^6, α, μ 10^6]]
(* {2.29388*10^7, 0.517994, 5.26508*10^6} *)

Also, the maximum likelihood estimator of $\mu$ is the minimum of the data: Min[vals].

2nd edit

To illustrate the issue of small sample size here are 1,000 simulations from the Pareto Type 2 distribution that shows with just a sample size of 15 how far (or how close) the estimates will be to the true values.

SeedRandom[567];
n = 15;
nSimulations = 1000;
parameter = {"k", "α", "μ"};
estimates = ConstantArray[{0, 0, 0}, nSimulations];
Do[x = RandomVariate[ParetoDistribution[22 10^6, 0.5, 5 10^6], n];
 estimates[[i]] = {k 10^6 , α, μ 10^6} /. 
   FindDistributionParameters[x, 
    ParetoDistribution[k 10^6, α, μ 10^6]]
 , {i, nSimulations}]
summary = 
  Table[{Mean[estimates[[All, i]]], Median[estimates[[All, i]]], 
    StandardDeviation[estimates[[All, i]]], 
    Quantile[estimates[[All, i]], 0.025],
    Quantile[estimates[[All, i]], 0.975], Min[estimates[[All, i]]], 
    Max[estimates[[All, i]]]}, {i, 3}];
TableForm[summary, 
 TableHeadings -> {{"k", "α", "μ"}, {"Mean", "Median", 
    "Std. Dev.", "Lower 2.5%", "Upper 2.5%", "Min", "Max"}}]

$$\left( \begin{array}{cccccccc} & \text{Mean} & \text{Median} & \text{Std. Dev.} & \text{Lower 2.5$\%$} & \text{Upper 2.5$\%$} & \text{Min} & \text{Max} \\ \text{k} & 4.46103\times 10^8 & 1.85451\times 10^7 & 7.94984\times 10^9 & 1.52004\times 10^6 & 1.51231\times 10^8 & 45484.7 & 2.17695\times 10^{11} \\ \alpha & 7.07993 & 0.491031 & 141.696 & 0.234873 & 1.67767 & 0.156161 & 4269.15 \\ \mu & 8.33604\times 10^6 & 7.11037\times 10^6 & 3.90537\times 10^6 & 5.04443\times 10^6 & 1.91982\times 10^7 & 5.00698\times 10^6 & 5.0886\times 10^7 \\ \end{array} \right)$$

So with a sample size of 15 one would expect the estimates of $alpha$ to be between 0.235 and 1.68 about 95% of the time.

$\endgroup$
  • $\begingroup$ No chance in increasing the sample size as I am working with real data where 15 points is considered a larger sample due to availability of the data. However, by applying knowledge on the behavior of the source of data and the impact of certain distribution parameters on the PDF of the distribution then assumptions can be introduced to produce a "reasonable" fit. $\endgroup$ – Edmund Sep 7 '18 at 11:59
  • $\begingroup$ Scaling the parameters looks promising. I will experiment with that. (+1) $\endgroup$ – Edmund Sep 7 '18 at 11:59
  • $\begingroup$ We both know that theoretically that scaling is irrelevant. But numerical methods can go wacko when parameters have wildly different magnitudes (which is the case here). I would not restrict the parameters but produce confidence intervals for the parameters as those intervals rather than the point estimates contain the information to check on preconceived notions. For this case (because the maximum likelihood estimate for $\mu$ is Min[vals], the usual Delta Method doesn't work. But a bootstrap approach would be your best bet to construct confidence intervals. $\endgroup$ – JimB Sep 7 '18 at 13:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.