5
$\begingroup$

I am wondering how you would obtain Scale and Shape parameter values on a Weibull Distribution's Confidence Interval bands (95% CI).

The following post nicely illustrates confidence interval bands with the PDF and CDF plots with bootstrap: Weibull distribution parameters $k$ and $c$ for wind speed data

Is there a simple way to extract the Scale and Shape parameter values from the confidence interval bands (plotted with 2 red bands)?

$\endgroup$
7
  • $\begingroup$ Choosing some values, plotting the cdf, and see if it falls within the band? $\endgroup$ Sep 20, 2015 at 11:57
  • $\begingroup$ I guess that would be an option if there is not a way to get the Scale and Shape values through R-code. You would think that there would be a better way but i am not sure. $\endgroup$
    – esh88
    Sep 21, 2015 at 13:15
  • $\begingroup$ Maybe I misunderstood your question. Do you have origina data? $\endgroup$ Sep 21, 2015 at 19:33
  • $\begingroup$ The dataset has 15 million records. It would be difficult to get it in here. I determined the Weibull distribution parameters to be Shape= 1.781096 and Scale = 33.669511 with the the fitdistr package. rw.small<-rweibull(100,shape=1.781096,scale=33.669511) $\endgroup$
    – esh88
    Sep 28, 2015 at 15:36
  • $\begingroup$ with 15 million observations the estimation uncertainty should be negligible ... $\endgroup$ Sep 28, 2015 at 15:37

1 Answer 1

5
$\begingroup$

Here is one way of finding confidence interval, using R and the CRAN package fitdistrplus (extending fitdist function from package mass). It uses maximum likelihood for the estimation (default method in fitdist) and likelihood profiling for the confidence intervals (this is implemented in function confint):

install.packages("fitdistrplus", dep=TRUE)
> set.seed(1234)   # For reproducibility
> x  <-  rweibull(10000,shape=1.6, scale=33)
> xfit  <-  fitdist(x,"weibull")
> xfit
Fitting of the distribution ' weibull ' by maximum likelihood 
Parameters:
      estimate Std. Error
shape  1.61569 0.01259297
scale 32.94673 0.21474443
> confint(xfit)
          2.5 %    97.5 %
shape  1.591008  1.640372
scale 32.525835 33.367618

Poking a little bit around, I don't think this is based on profiling, but is simply the asymptotic interval:

> class(xfit)
[1] "fitdist"
> methods(class="fitdist")
[1] coef     logLik   plot     print    quantile summary  vcov    
see '?methods' for accessing help and source code
> methods(confint)
[1] confint.default       confint.glm*          confint.lm           
[4] confint.nls*          confint.polr*         confint.profile.glm* 
[7] confint.profile.nls*  confint.profile.polr*
see '?methods' for accessing help and source code

so, wee see ... confint doesn't have a method for objects of class "fitdist", so the default method is used. That gives the asymptotic interval. Profiling is done by the profile function (in mass), which do not have a method for "fitdist" objects:

> methods(profile)
[1] profile.glm*  profile.nls*  profile.polr*
see '?methods' for accessing help and source code
> prof  <-  profile(xfit)
Error in UseMethod("profile") : 
  no applicable method for 'profile' applied to an object of class "fitdist"
> 

But we can use bootstrapping (with your 15 million data, will take a very long time ...):

> xfit.boot  <-  bootdist(xfit)   # From fitdistrplus
> xfit.boot
Parameter values obtained with parametric bootstrap 
     shape    scale
1 1.620763 32.59287
2 1.609556 32.61217
3 1.620512 33.53286
4 1.612359 32.96155
5 1.606430 33.03867
6 1.600184 32.60450
> plot(xfit.boot)
> summary(xfit.boot)
Parametric bootstrap medians and 95% percentile CI 
         Median      2.5%     97.5%
shape  1.615281  1.592227  1.638437
scale 32.953338 32.536000 33.394774

The plot function above simple gives a scatterplot of the bootstrapped parameter values, shown below.

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.