Parameter Estimation for 3-parameter Weibull Distribution
- To: mathgroup@smc.vnet.net
- Subject: [mg11914] Parameter Estimation for 3-parameter Weibull Distribution
- From: Bob Hanlon <BobHanlon@aol.com>
- Date: Thu, 9 Apr 1998 00:33:16 -0400
(* 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