bootstrapconfidence intervalprobabilityrweibull distribution
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)?
There is no good reason not to do such a generalization, which would unite the Weibull and inverse Weibull distributions. So reasons must be historical or accidental.
Also see the comments thread.
Here's the fitted pdf and cdf (Weibull) for each of locations 1 to 3:
Let's break down what we need to do here, keeping in mind that the end goal is to estimate the cumulative proportion of area planted with a certain crop at some value for the random variable time $X$:
The first step is to fit a distribution (e.g. $Weibull\left(a:= \text{shape},b := \text{scale}\right)$ or $Beta\left(\alpha,\beta\right)$).
You made an error in fitting the data on a Weibull distribution because the function fitdist (which uses MLE by default) expects a vector of observations, which would be 'time $X_i$ when some area $i$ was planted with a certain crop'.
Here is some R that turns the cumulative data into a vector of observations, which can then be used to fit the distribution with fitdist:
dat <- dat[order(dat$loc.id,dat$year.id),]
dat$loc.id <- factor(dat$loc.id)
dat$year.id <- factor(dat$year.id)
library(plyr)
#hacking the original data from the cumulative
to.density <- function(cdf) { #given input (a cdf vector), returns respective densities
y <- c(); cdf <- cdf[,'cum.per.plant']
for (cum in seq(1,NROW(cdf))) {y <- c(y, ifelse(cum==1, cdf[cum], cdf[cum]-cdf[cum-1]))}
return(y)
}
#apply to each location X year in dat, then kbind density data to dat
densities <- dlply(dat, .(loc.id, year.id), to.density)
dat <- cbind(dat,dens=unlist(densities))
#forming a dataset of the variable planting time
#--because the proportions go to the thousandths place, a pseudo-sample of time.id's with
#--n=1000 should allow usage of fitdist
rep.vec <- function(df) {
y <- c()
for (row in rownames(df)) {y <- c(y, rep(df[row,'time.id'], df[row,'dens']*1000))}
return(y)
}
#apply to each location in dat - may want to include year.id in the split (excluded)
time.sample <- dlply(dat, .(loc.id), rep.vec)
Be sure to assess the 'fit' of your estimated parameters (perhaps by calculating the mean square error) on the data.
Your code does not indicate that each location and year will be fitted with its own distribution, although you suggest you want to calculate the cumulative area planted for a given location or year.
Here is some R for fitting each location:
#fitting separate Weibull distributions for each loc.id (may want to include year.id in
#the split)
library(fitdistrplus)
fit.weibull <- function(loc) {
y <- summary(fitdist(loc,'weibull'))[[1]]
return(y)
}
params <- lapply(time.sample, fit.weibull) #apply to each element in time sample
Finally, consider the inclusion of a location parameter, which shifts the graph of the pdf in a negative or positive direction along the x-axis; this should be appropriate because in many locations, no area gets plotted until the $X=2^{nd}$ week. The inclusion of a location parameter would likely improve the fit of your distributions.
Use the fitted cdf (with the parameters informed by the previous step) to predict the cumulative proportion of area planted on a certain day for a given location.
It is important to understand what you are doing when you want to rescale the random variable time $X$ to represent days rather than weeks, which simply involves dividing the data (a vector of time observations) by 7. This can be seen in the equation for the Weibull pdf itself.
I've demonstrated it through a few lines of R:
#creating predictive model - on a daily rather than weekly domain
predict1.cum.plant <- function(day, loc.id, params) {
#'scaling' 7x the scale parameter
pweibull(day, shape=params[[loc.id]][[1]], scale=params[[loc.id]][[2]]*7)
}
predict2.cum.plant <- function(day, loc.id, params) {
#'scaling' 1/7x the random variable time
pweibull(day/7, shape=params[[loc.id]][[1]], scale=params[[loc.id]][[2]])
}
#example - the models are equal, which is obvious from the cdf equation!
predict1.cum.plant(35, 2, params)
predict2.cum.plant(35, 2, params)
pweibull(5, shape=params[[2]][[1]], scale=params[[2]][[2]]) #by week
Some supplemental code of mine can be found here.
Very similar methods can be used to fit a Beta distribution.
Best Answer
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):
Poking a little bit around, I don't think this is based on profiling, but is simply the asymptotic interval:
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:
But we can use bootstrapping (with your 15 million data, will take a very long time ...):
The plot function above simple gives a scatterplot of the bootstrapped parameter values, shown below.