Varying values of parameters and sample size won't work

Question

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] 

Answer 1

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"}]}}]

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).