MathGroup Archive 1998

[Date Index] [Thread Index] [Author Index]

Search the Archive

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