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

    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.

• You don't need that table - you can write logl = LogLikelihood[f[α, λ], x]. Also why are you creating κ many batches of n samples ? Apr 11, 2021 at 14:13
• Also you're minimizing the log likelihood which is wrong - you're supposed to be maximizing it. Apr 11, 2021 at 14:21
• 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, 2021 at 14:25
• @flinty i had given k=1000 for simulating "n" 1000 times. Apr 12, 2021 at 5:46

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

• Thank You for your support. Apr 12, 2021 at 5:56