# Parameter Estimation for 3-parameter Weibull Distribution

```(* Generalizes 2-parameter Weibull distribution to  *) (*     3
parameters                                 *) (*
WeibullDistribution[alpha, beta] =               *) (*
WeibullDistribution[0, alpha, beta]          *) (* Includes Moment,
Median, and Mode statistics for *) (* 2- and 3-parameter Weibull
distributions         *)

Needs["Statistics`ContinuousDistributions`"];
Needs["Statistics`DescriptiveStatistics`"];
Unprotect[WeibullDistribution];
WeibullDistribution/:
PDF[WeibullDistribution[a_,b_,m_],x_]/; x<a := 0;
WeibullDistribution/:
PDF[WeibullDistribution[a_, b_, m_], x_] :=
b/m * ((x-a)/m)^(b - 1) * Exp[-((x-a)/m)^b];  WeibullDistribution/:
Domain[WeibullDistribution[a_, b_, m_]] :=
Interval[{a, Infinity}]
WeibullDistribution/:
CDF[WeibullDistribution[a_, b_, m_], x_] :=
1 - Exp[-((x-a)/m)^b];
WeibullDistribution/:
Mean[WeibullDistribution[a_, b_, m_]] :=
a + m * Gamma[1 + 1/b];
WeibullDistribution/:
Variance[WeibullDistribution[a_, b_, m_]] :=
m^2 * (Gamma[1 + 2/b] - Gamma[1 + 1/b]^2); WeibullDistribution/:
StandardDeviation[
WeibullDistribution[a_, b_, m_]] :=
m * Sqrt[(Gamma[1 + 2/b] - Gamma[1 + 1/b]^2)];
WeibullDistribution/:
Moment[WeibullDistribution[a_, b_, m_], n_Integer] :=
Sum[Binomial[n,k]*a^(n-k)*m^k*Gamma[1+k/b],{k,0,n}];
WeibullDistribution/:
Moment[WeibullDistribution[b_, m_], n_Integer] :=
m^n * Gamma[1 + n/b];
WeibullDistribution/:
Mode[WeibullDistribution[a_:0, b_, m_]] :=
a + m * (1 - 1/b)^(1/b);
WeibullDistribution/:
Median[WeibullDistribution[a_:0, b_, m_]] :=
a + m * Log[2]^(1/b);
Protect[WeibullDistribution];

Likelihood[dist_, x_List] :=
Times @@ (PDF[dist, #]& /@ x);
LogLikelihood[dist_, x_List] :=
Log[Times @@ (PDF[dist, #]& /@ x)];

data = 3. + RandomArray[WeibullDistribution[3.,1.],100]; minData =
Min[data];
dist = WeibullDistribution[a,b,m];
llf = LogLikelihood[dist, data] // PowerExpand; eqnLLF = {D[llf, a] ==
0, D[llf, b] == 0, D[llf, m] == 0}; eqnEst = {Mean[dist] ==
Mean[data],Variance[dist] == Variance[data],
Median[dist] == Median[data]};

Initial estimate of parameters:

pEst = FindRoot[eqnEst, {a,minData}, {b, 1.}, {m, 1.}]

Maximum likelihood estimate of parameters:

pMLE = FindRoot[eqnLLF,{a,a/.pEst},{b,b/.pEst},{m,m/.pEst}] /.
x_Complex -> Re[x]

dist = WeibullDistribution[a, b, m] /. pMLE; xmin = a /. pMLE;
xmax = Mean[dist] + 3 StandardDeviation[dist]; xL = (xmin+3xmax)/4;
Plot[{PDF[dist,x], CDF[dist,x]}, {x,xmin,xmax},
AxesLabel->{"x", None},
PlotLabel->"PDF and CDF for Weibull Distribution",
Epilog->{
Text[StringForm[
"Mean = ``", NumberForm[Mean[dist]//N,3]], {xL,.8},{-1,0}],
Text[StringForm[
"StdDev = ``", NumberForm[StandardDeviation[dist]//N,3]],
{xL, .65},{-1,0}],
Text[StringForm[
"Mode = ``", NumberForm[Mode[dist]//N,3]], {xL, .5},{-1,0}],
Text[StringForm[
"Median = ``", NumberForm[Median[dist]//N,3]],
{xL, .35},{-1,0}]}];

Bob Hanlon

```

• Prev by Date: ANNOUNCE: New german book on Mathematica 3
• Next by Date: Re: Inserted Objects wont Print!!
• Prev by thread: ANNOUNCE: New german book on Mathematica 3
• Next by thread: Formatting and Style Sheet woes