1
$\begingroup$

1.Hello , I have customized a distribution, and then trying to estimate the parameters for smaller values of sample size say n=10,30,40,50 and 60.It runs for n=10 and n=1000 and then it does not run for any other values of either n or for other values of parameters. any suggestion? Thanks. 

    Clearall[myDist] ;
    Clearall[sigma, alpha, theta];
    Clear[W]
    ClearSystemCache[]
    myDist[sigma_, alpha_, theta_] = 
      ProbabilityDistribution[
       alpha theta/
         sigma (1 - Exp[-(x/sigma)^alpha])^(theta - 1) Exp[-  (x/sigma)^
       alpha ]*(x/sigma)^(alpha - 1), {x, 0, Infinity}, 
       Assumptions -> sigma > 0 && alpha > 0 && theta > 0];
    DistributionParameterAssumptions[myDist[sigma, alpha, theta]];
     multiparams = 
     Monitor[Table[W = RandomVariate[myDist[3, 2, 1], 1000];
       params = 
        FindDistributionParameters[W, myDist[a, b, c], AccuracyGoal -> 5,PrecisionGoal -> 5, WorkingPrecision -> 5], {i, 5}], i]
Improve this question
$\endgroup$
4
  • $\begingroup$ Sorry the last line of the code got cut while posting, you can complete it with "PrecisionGoal -> 5, WorkingPrecision -> 5], {i, 5}], i] " $\endgroup$
    – nutan
    Oct 3, 2017 at 16:02
  • $\begingroup$ You have a few issues. For first two lines it is ClearAll not Clearall. You should be using SetDelayed for myDist instead of Set (see answer 18487). $\endgroup$
    – Edmund
    Oct 3, 2017 at 23:29
  • $\begingroup$ @Edmund, Thanks for response. I made the suggested changes and corrections. i still face the same problem. any other suggestions? $\endgroup$
    – nutan
    Oct 5, 2017 at 11:23
  • $\begingroup$ @Edmund: Hello I made th changes as per suggestions and ran the program. I am still having the same problem . it does not run frothier values of n, thanks in advance. $\endgroup$
    – nutan
    Oct 13, 2017 at 9:55

1 Answer 1

Reset to default
1
$\begingroup$

Update

I've increased the number of simulations so that the high correlations among the parameter estimators based on the observed Fisher Information matrix can be justified. (And I've deleted my inappropriately sarcastic comment from below the answer.)

I think the main problem is that for small sample sizes, the maximum likelihood estimators don't converge all of the time. I'd also suggest not using FindDistributionParameters if for no other reason one can't obtain standard errors of the estimates using that function.

Below is some code to perform the maximum likelihood estimation (although even 5,000 iterations is not always enough) and obtain the covariance and correlation matrices. Note that the estimators are all very highly correlated with each other and that many times goes hand-in-hand with lack of convergence.

myDist[sigma_, alpha_, theta_] := 
  ProbabilityDistribution[
   alpha theta/
     sigma (1 - Exp[-(x/sigma)^alpha])^(theta - 1) Exp[-(x/sigma)^
       alpha]*(x/sigma)^(alpha - 1), {x, 0, Infinity}, 
   Assumptions -> sigma > 0 && alpha > 0 && theta > 0];

n = 60;  (* Sample size *)
nSimulations = 50; (* Number of simulations *)
myD = myDist[a, b, c];  (* Get distribution in terms of parameters a, b, and c *)

(* Loop through the estimation procedure *)
SeedRandom[12345];
estimates = ConstantArray[0, nSimulations];
Do[data = RandomVariate[myDist[3, 2, 1], n];

 (* Log of the likelihood *)
 logL = LogLikelihood[myD, data];

 (* Maximum likelihood estimates *)
 mle = FindMaximum[{logL, 
    DistributionParameterAssumptions[myDist[a, b, c]]}, {{a, 3}, {b, 2}, {c, 1}}, 
    MaxIterations -> 5000];

 (* Parameter estimator covariance and correlation matrices *)
 cov = -Inverse[(D[logL, {{a, b, c}, 2}]) /. mle[[2]]];
 cor = Table[cov[[i, j]]/(cov[[i, i]]^0.5 cov[[j, j]]^0.5), {i, 3}, {j, 3}];

 (* Save results *)
 estimates[[i]] = Join[mle, {cov}, {cor}],
 {i, nSimulations}]

A typical correlation matrix for the parameter estimators is

estimates[[3, 4]] // MatrixForm

$$\left( \begin{array}{ccc} 1. & 0.932378 & -0.953094 \\ 0.932378 & 1. & -0.958199 \\ -0.953094 & -0.958199 & 1. \\ \end{array} \right)$$

If we plot the estimates from the 50 simulations we see that the estimated correlations are high and of the same sign as estimated by the observed Fisher Information matrix.

e = {a, b, c} /. # & /@ estimates[[All, 2]];
GraphicsGrid[{{ListPlot[e[[All, {1, 2}]], PlotRange -> All, 
    Frame -> True,
    FrameLabel -> {"Estimate of a", "Estimate of b"}],
   ListPlot[e[[All, {1, 3}]], PlotRange -> All, Frame -> True,
    FrameLabel -> {"Estimate of a", "Estimate of c"}]},
  {ListPlot[e[[All, {2, 3}]], PlotRange -> All, Frame -> True,
    FrameLabel -> {"Estimate of b", "Estimate of c"}]}}]

Estimates from simulations

Because of the nonlinear association of estimators of $a$ and $b$ with $c$, using the correlation (a measure of linear fit) is probably not the best measure to use here. However, the "sign" and "strength" of the relationships is certainly similar to using the correlation coefficient: The estimators of $a$ and $b$ are positively related and both are negatively related with the estimator of $c$ (at least for the true values of $a$, $b$, and $c$ being 3, 2, and 1, respectively).

Improve this answer
$\endgroup$
8
  • $\begingroup$ wondering if it is posible to know the name of "automatic" method in FindMaximum. $\endgroup$
    – nutan
    Oct 16, 2017 at 12:21
  • 1
    $\begingroup$ Normally I would say that you should try the list of methods given in the documentation but here if one lists any of the methods other than Automatic, one gets the error message FindMaximum::ucmtd: "Method ->Newton can only be used for unconstrained problems." Removing the constraints either results in non-convergence or the same exact estimates as Automatic for more than one of the methods. So my response is "I don't know and I don't know how I might find out in this case." But I think your bigger concern should be the small sample size and the large correlations. $\endgroup$
    – JimB
    Oct 16, 2017 at 14:44
  • $\begingroup$ Thanks. High correlation among the estimates is okay, this is a distribution with one scale and two shape parameters. What is not correct is the positive correlation between the first two parameters. I checked with my density function and things look alright. One can notice that for third parameter c=1 this is a Weibull and it is known that the scale and shape parameters of Weibull are highly negatively correlated. $\endgroup$
    – nutan
    Oct 23, 2017 at 20:01
  • $\begingroup$ For normal distribution, yes the two parameters are independent but for example if we look at the MLE's of Kumaraswamy distribution, if we look at their Fisher information matrix and one can see the correlation between the parameters. Reference: Improved point estimation for the Kumaraswamy distribution Artur J. Lemonte. Journal of Statistical Computation and Simulation Vol. 81 , Iss. 12,2011 (end page 1973 and beginning of page 1974) the estimators of the two parameters are not independent. Thanks for your comment. $\endgroup$
    – nutan
    Oct 23, 2017 at 23:25
  • 1
    $\begingroup$ I don't have immediate access to that article but I'll be able to get it in a few days. In the meantime I get a positive (and not negative) correlation of around 0.3 for fitting a two-parameter Weibull with the following statements: x = RandomVariate[WeibullDistribution[3, 2], 1000]; logL = LogLikelihood[WeibullDistribution[a, b], x]; sol = FindMaximum[logL, {{a, 3}, {b, 2}}]; cov = -Inverse[(D[logL, {{a, b}, 2}]) /. sol[[2]]]; cor = Table[cov[[i, j]]/(cov[[i, i]] cov[[j, j]])^0.5, {i, 2}, {j, 2}]; cor // MatrixForm. $\endgroup$
    – JimB
    Oct 24, 2017 at 0:26

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.