0
    f[\[Alpha]_, \[Lambda]_] := 
     ProbabilityDistribution[(
      4*\[Lambda]*\[Alpha]*
       x^(\[Alpha] - 1)*(1 - x)*(2 - x)^(\[Alpha] - 
        1)*(1 - x^\[Alpha]*(2 - x)^\[Alpha])^(\[Lambda] - 
        1))/(1 + (1 - x^\[Alpha]*(2 - x)^\[Alpha])^\[Lambda])^2, {x, 0, 
       1}, Assumptions -> \[Alpha] > 0 && \[Lambda] > 0]
    n = 10
    \[Kappa] = 1000
    x = Table[RandomVariate[f[2, 0.5], n], {i, \[Kappa]}]
    logl = Table[Sum[((
        1 - (1 - x[[i]]^\[Alpha]*(2 - x[[i]])^\[Alpha])^\[Lambda])/(
        1 + (1 - x[[i]]^\[Alpha]*(2 - x[[i]])^\[Alpha])^\[Lambda]) - i/(
        n + 1))^2, {i, n}],{i, \[Kappa]}]

FindMinimum[{logl, \[Alpha] > 0, \[Lambda] > 0}, {{\[Alpha], 2}, {\[Lambda], 0.5}}]

I tried the above code to find the estimates and bias after simulating the data. For that ,I generated the sample from pdf. Then specified the function which is to be minimised, but i am unable to find the estimates. Could anyone please help me to correct the code to find the bias and estimates. As the sample size increases, bias to be decreased.

Akhilraj N S
  • 115
  • 4
  • You don't need that table - you can write logl = LogLikelihood[f[α, λ], x]. Also why are you creating κ many batches of n samples ? – flinty Apr 11 '21 at 14:13
  • Also you're minimizing the log likelihood which is wrong - you're supposed to be maximizing it. – flinty Apr 11 '21 at 14:21
  • 1
    A few other things: (1) you run into problems when x is used in the definition of the probability distribution and x is used to store the random samples if you don't clear x each time you run the code and (2) FindMaximum (good catch @flinty) takes a single function rather than a list of functions. So using FindMaximum[{#, ....]& /@ logl would get you the separate results. – JimB Apr 11 '21 at 14:25
  • @flinty i had given k=1000 for simulating "n" 1000 times. – Akhilraj N S Apr 12 '21 at 05:46

1 Answers1

1

Using @flinty 's suggestions and a smaller number of simulations, here is how you can capture the estimates for $\alpha$ and $\lambda$ to estimate bias in those parameter estimators:

(* Define probability distribution *)
f[α_, λ_] := 
 ProbabilityDistribution[(4*λ*α*xx^(α - 1)*(1 - xx)*(2 - xx)^(α - 1)*(1 - xx^α*(2 - xx)^α)^(λ - 1))/
  (1 + (1 - xx^α*(2 - xx)^α)^λ)^2, {xx, 0, 1}, 
  Assumptions -> α > 0 && λ > 0]

(* Generate κ samples each of size n *)
n = 10;
κ = 100;
SeedRandom[12345];
α0 = 2;
λ0 = 0.5;
x = RandomVariate[f[α0, λ0], {n, κ}];


(* Construct log likelihood for each set of samples  *)
logl = LogLikelihood[f[α, λ], #] & /@ x;


(* Find maximum likelihood estimates *)
results = {α, λ} /. FindMaximum[{#, α > 0, λ > 0}, {{α, α0}, {λ, λ0}}][[2]] & /@ logl;


(* α Summary )
meanα = Mean[results[[All, 1]]] ( 2.2894658167785593*) biasα = meanα - α0 (* 0.2894658167785593 )
StandardDeviation[results[[All, 1]]]   ( 0.5252458945336623` *)


(* λ Summary )
meanλ = Mean[results[[All, 2]]]  ( 0.5350471371576843*) biasλ = meanλ - λ0 (* 0.03504713715768426 )
StandardDeviation[results[[All, 2]]]  ( 0.06487409567245328` *)
JimB
  • 41,653
  • 3
  • 48
  • 106