Modifying the Parameter Estimation Method of wblfit Function

15 views (last 30 days)
By default, the wblfit function returns the 'maximum likelihood estimate' as describied below:
parmhat = wblfit(data) returns the maximum likelihood estimates , parmhat, of the parameters of the Weibull distribution given the values in the vector data, which must be positive. parmhat is a two-element row vector: parmhat(1) estimates the Weibull parameter a, and parmhat(2) estimates the Weibull parameter b, in the pdf.
How can I change the parameter estimation method from maximum likelihood estimate to rank regression method of the least square curve fit? I've spent 2 days trying to figure this out and have had no luck. The 'rank regression method' is very common when it comes to weibull analysis on small sample sizes. Any help would be greatly appreciated.

Answers (2)

Tom Lane
Tom Lane on 4 Apr 2015
I think you should not modify wblfit, but just do it directly or write a separate function. Here's an example of how to do what I think you want:
% Random sample with paramaters 1000 and 1.3
n = 1000;
x = wblrnd(1000,1.3,n,1);
% Fit by maximum likelihood
p = wblfit(x)
% Create probability plot on log scale
z = wblinv(((1:n)-.5)/n,1,1)'; % idealized sample from standard Weibull
plot(log(z),sort(log(x)),'bx')
% Fit points on plot by least squares
b = polyfit(log(z),sort(log(x)),1)
line(log(z),b(1)*log(z)+b(2),'LineStyle',':')
% Transform to Weibull parameter scale
p1 = 1/b(1)
p2 = exp(b(2))
  1 Comment
Ali
Ali on 5 Apr 2015
Tom,
Thank you for your feedback! The result was close but not exactly the right answer.
For a sample of failure times of [30 49 82 90 96]. The Eta would be 79.8 and the Beta would be 2.165 using the rank regression model. Using your algorithm, I arrive at an Eta 78.723 and a Beta of 2.465, which is very close.
I am trying to go through and determine if I can make the necessary adjustments to arrive at the proper solution.

Sign in to comment.


Dutchie
Dutchie on 1 Nov 2016
Solution, is the calculation of the Z-variable.
z = wblinv(((1:n)-.3)/(n+0.4),1,1)'
results are correct.. as per the source.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!