-1

I used bootstrapping to obtain confidence intervals of a Weibull distribution. Then I plotted the Confidence Bands in a plot.

Code is below:

set.seed(123)
rw.small<-rweibull(100,shape=1.781096,scale=33.669511)
xs <- seq(0,100, len=500)

boot.pdf <- sapply(1:100, function(i) {
     xi <- sample(rw.small, size=length(rw.small), replace=TRUE)
     MLE.est <- suppressWarnings(fitdist(xi, distr="weibull",lower=0))  
     dweibull(xs, shape=MLE.est$estimate["shape"],  scale = MLE.est$estimate["scale"])
 })

par(bg="white",las=1,cex=1.2)

plot(xs, boot.pdf[, 1], type="l", col=rgb(.6, .6, .6, .1), ylim=range(boot.pdf),
      xlab="Note Life (months)", ylab="Probability density",main= "Probability Distribution")

for(i in 2:ncol(boot.pdf)) lines(xs, boot.pdf[, i], col=rgb(.6, .6, .6, .1))
quants <- apply(boot.pdf, 1, quantile, c(0.025, 0.5, 0.975))
min.point <- apply(boot.pdf, 1, min, na.rm=TRUE)
max.point <- apply(boot.pdf, 1, max, na.rm=TRUE)
lines(xs, quants[1, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[3, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[2, ], col="darkred", lwd=2)

Two questions: 1. How do i obtain the Shape and Scale Parameter values from the lower and upper bounds of the confidence interval plotted (i.e. the two red lines)?

lines(xs, quants[1, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[3, ], col="red", lwd=1.5, lty=2)
  1. I would like to plot the same chart with R's ggplot package. Any ideas as to how i might do this with ggplot syntax?
StupidWolf
  • 45,075
  • 17
  • 40
  • 72
esh88
  • 61
  • 1
  • 6

1 Answers1

2

I'm not sure about question 1.

Here is some code for question 2 though.

library (fitdistrplus) #you forgot to mention this
library (ggplot2)
library (tidyr)

dat <- gather(as.data.frame(boot.pdf), i, y, -1)
dat2 <- gather(as.data.frame(t(quants)), quantile, y)
dat$x <- dat2$x <- xs

ggplot(dat, aes(x, y, group = i)) +
  geom_line(col = 'grey') +
  geom_line(data = dat2, aes(group = quantile, col = quantile), size = 1) +
  scale_color_manual(values  = c('2.5%' = 'red', '50%' = 'darkred', '97.5%' = 'red')) +
  theme_bw() + xlab("Note Life (months)") + ylab("Probability density") + 
  ggtitle("Probability Distribution")

enter image description here

Axeman
  • 32,068
  • 8
  • 81
  • 94
  • Thanks Axeman. For question 1, i am really trying to find the Weibull Shape and Scale parameter values of the 2.5% and 97.5% curves. For example, the original model had shape=1.781096,scale=33.669511. – esh88 Sep 17 '15 at 18:49
  • Yes I understand, but I'm not sure how to do that. You could have a look at the `fitdistr` function in the `MASS` package. – Axeman Sep 18 '15 at 08:18
  • Or look here: http://stackoverflow.com/questions/29374915/data-series-how-can-i-fit-a-distribution-in-r – Axeman Sep 18 '15 at 08:21
  • `fits <- apply(quants[, -1], 1, fitdist, 'weibull')` This runs, but the estimates are very different from what you'd expect. I guess something is wrong in your code, but I'm not sure what. – Axeman Sep 18 '15 at 08:28