On Jun 23, 2011, at 9:44 AM, Adriana Bejarano wrote:
> Dear R gurus,
>
> I have the following code, but I still not know how to estimate and
> extract
> confidence intervals (95%CI) from resampling.
>
If you have a distribution of values, say "resamp.stat", of a
statistic from a properly performed resampling operation you can
extract and display easily the 5th and 95th percentiles.
CI.stat <- quantile(resamp.stat, c(0.05, 0.95) )
CI.stat
Note: I do not think that 100 replications would generally be
sufficient for final work, although its probably acceptable for testing.
Your code as posted initially threw a bunch of errors since you did
not include library(MASS), but fixing that fairly obvious problem
shows that you can draw a nice plot. However, it remains unclear what
statistic of what distribution you desire to assess. Mean, median, ...
of what?
I do not think the error that arose on my machine from the wrapped
text here:
#draw random numbers from a weibull distribution 100 times with
... shape=xwei.shape, scale = xwei.scale -> error
..... was causing any problem, but there were a bunch of warnings that
ought to be investigated:
Right after the loop I see ten of these:
Warning messages:
1: In dweibull(x, shape, scale, log) : NaNs produced
--
David
> Thanks!
>
> ~Adriana
>
> #data
> penta<-
> c
>
(770,729,640,486,450,410,400,340,306,283,278,260,253,242,240,229,201,198,190,186,180,170,168,151,150,148,147,125,117,110,107,104,85,83,80,74,70,66,54,46,45,43,40,38,10
> )
> x<-log(penta+1)
> plot(ecdf(x), ylab="Probability", xlab="Concentration
(Ln+1)")
>
> x.wei<-fitdistr(x,"weibull")
> x.wei
> shape scale
> 6.7291685 5.3769965
> (0.7807718) (0.1254696)
> xwei.shape <- x.wei$estimate[[1]]
> xwei.scale <- x.wei$estimate[[2]]
>
> x.wei<-fitdistr(x,"weibull")
> x.wei
> xwei.shape <- x.wei$estimate[[1]]
> xwei.scale <- x.wei$estimate[[2]]
> curve(pweibull(x, shape=xwei.shape, scale =
> xwei.scale,lower.tail=TRUE,
> log.p=FALSE), add=TRUE,col='green',lwd=3)
>
> #draw random numbers from a weibull distribution 100 times with
> shape=xwei.shape, scale = xwei.scale
> draw <- lapply(1:100, function(.x){
> out<-rweibull(x, shape=xwei.shape, scale = xwei.scale)
> })
> newx<- data.frame(draw)
>
> colnames(newx)<-paste("x", 1:100, sep = "")
> newmat<-data.matrix(newx)
>
> # matrix of coefficients
> rownum=2
> colnum=100
> ResultMat<-matrix(NA, ncol=colnum, nrow=rownum)
> rownum2=45
> colnum2=100
> ResultMat2<-matrix(NA, ncol=colnum2, nrow=rownum2)
>
> #loop through each column in the source matrix
> for (i in 1:100)
> {
> sel_col<-newmat[col(newmat)==i]
> {ResultMat[,i]<-coef(fitdistr(sel_col,"weibull"))}
> xwei.shape<- ResultMat[1,i]
> xwei.scale<- ResultMat[2,i]
> curve(pweibull(x, shape=xwei.shape, scale=xwei.scale, lower.tail=TRUE,
> log.p = FALSE), add=TRUE, col='grey',lwd=0.5)
> ResultMat2[,i]<-pweibull(x, shape=xwei.shape, scale >
xwei.scale,lower.tail=TRUE, log.p=FALSE)
> }
>
> ## convert dataframe to numeric
> MatOut<- as.matrix(ResultMat2)
> mode(MatOut) <- "numeric"
>
> # initiate variables
> mm<-ml<-mu<-rep(0,length(MatOut[,1]))
>
> # mean and upper/lower quantiles
> for(i in 1:length(MatOut[,1])){
> mm[i]<- mean(MatOut[i,])
> ml[i]<- quantile(MatOut[i,], 0.025, na.rm=TRUE)
> mu[i]<- quantile(MatOut[i,], 0.975, na.rm=TRUE)
> }
> #lines(x, mm, col="black")
> lines(x, ml, col="blue", lwd=2)
> lines(x, mu, col="blue", lwd=2)
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
David Winsemius, MD
West Hartford, CT