How to find estimates and bias from simulated data using method of least square estimation

Question

    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.

Answer 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` *)