HI All: Does anyone know the code behind the qbeta function in R? I am using it to calculate exact confidence intervals and I am getting 'NaN' at places I shouldnt be. Heres the simple code I am using: k<-3> x<-NULL > p<-rbeta(k,3,3)# so that the mean nausea rate is alpha/(alpha+beta) > min<-10 > max<-60 > n<-as.integer(runif(3,min,max)) > for(i in 1:k)+ x<-cbind(x,rbinom(5,n[i],p[i]))> > # Exact Confidence Interval > > l_cl_exact<-qbeta(.025, x, n-x+1)Warning message: In qbeta(p, shape1, shape2, lower.tail, log.p) : NaNs produced> u_cl_exact<-qbeta(.975, x+1, n-x)Warning message: In qbeta(p, shape1, shape2, lower.tail, log.p) : NaNs produced> x[,1] [,2] [,3] [1,] 8 12 14 [2,] 5 15 13 [3,] 5 12 12 [4,] 8 21 12 [5,] 8 14 12> n[1] 10 36 31> l_cl_exact[,1] [,2] [,3] [1,] 0.44390454 0.2184996 0.2314244 [2,] 0.04667766 NaN 0.2454760 [3,] 0.05452433 0.1855618 NaN [4,] 0.44390454 0.4862702 0.1855618 [5,] 0.10115053 NaN 0.2184996 Thanks for your help. Anamika [[alternative HTML version deleted]]
On Mar 9, 2012, at 2:48 PM, Anamika Chaudhuri wrote:> HI All: > > Does anyone know the code behind the qbeta function in R?Well, yes, but don't you think it would be wise to question whether your code might be the problem rather than the R code?> I am using it to calculate exact confidence intervals and I am getting > 'NaN' at places I shouldnt be. Heres the simple code I am using: > > k<-3 >> x<-NULL >> p<-rbeta(k,3,3)# so that the mean nausea rate is alpha/(alpha+beta) >> min<-10 >> max<-60 >> n<-as.integer(runif(3,min,max)) >> for(i in 1:k) > + x<-cbind(x,rbinom(5,n[i],p[i]))Isn't this going to make x get longer with each pass through the loop? I think your parameter are then going to be interpreted as values of "x". Looks like "user error" to me. -- David>> >> # Exact Confidence Interval >> >> l_cl_exact<-qbeta(.025, x, n-x+1) > Warning message: > In qbeta(p, shape1, shape2, lower.tail, log.p) : NaNs produced >> u_cl_exact<-qbeta(.975, x+1, n-x) > Warning message: > In qbeta(p, shape1, shape2, lower.tail, log.p) : NaNs produced >> x > [,1] [,2] [,3] > [1,] 8 12 14 > [2,] 5 15 13 > [3,] 5 12 12 > [4,] 8 21 12 > [5,] 8 14 12 >> n > [1] 10 36 31 >> l_cl_exact > [,1] [,2] [,3] > [1,] 0.44390454 0.2184996 0.2314244 > [2,] 0.04667766 NaN 0.2454760 > [3,] 0.05452433 0.1855618 NaN > [4,] 0.44390454 0.4862702 0.1855618 > [5,] 0.10115053 NaN 0.2184996 > > Thanks for your help. > Anamika > > [[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
Take a look at n-x+1, the second parameter to the beta distribution:> n <- c(10, 45, 38) > x <- rbind(c( 7, 45, 31),+ c(10, 40, 35), + c( 9, 44, 33), + c( 8, 44, 31), + c( 8, 45, 36))> n - x + 1[,1] [,2] [,3] [1,] 4 -6 15 [2,] 36 -29 4 [3,] 30 2 -22 [4,] 3 -5 15 [5,] 38 -34 3 You probably intended> sweep(1-x, MAR=2, n, `+`)[,1] [,2] [,3] [1,] 4 1 8 [2,] 1 6 4 [3,] 2 2 6 [4,] 3 2 8 [5,] 3 1 3 If you had been unlucky, none of the entries in n-x+1 would have been negative and you would have received no warning from qbeta to give a hint that n-x+1 was not working as expected. Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com> -----Original Message----- > From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf > Of Anamika Chaudhuri > Sent: Friday, March 09, 2012 11:48 AM > To: r-help at r-project.org > Subject: [R] qbeta function in R > > HI All: > > Does anyone know the code behind the qbeta function in R? > I am using it to calculate exact confidence intervals and I am getting > 'NaN' at places I shouldnt be. Heres the simple code I am using: > > k<-3 > > x<-NULL > > p<-rbeta(k,3,3)# so that the mean nausea rate is alpha/(alpha+beta) > > min<-10 > > max<-60 > > n<-as.integer(runif(3,min,max)) > > for(i in 1:k) > + x<-cbind(x,rbinom(5,n[i],p[i])) > > > > # Exact Confidence Interval > > > > l_cl_exact<-qbeta(.025, x, n-x+1) > Warning message: > In qbeta(p, shape1, shape2, lower.tail, log.p) : NaNs produced > > u_cl_exact<-qbeta(.975, x+1, n-x) > Warning message: > In qbeta(p, shape1, shape2, lower.tail, log.p) : NaNs produced > > x > [,1] [,2] [,3] > [1,] 8 12 14 > [2,] 5 15 13 > [3,] 5 12 12 > [4,] 8 21 12 > [5,] 8 14 12 > > n > [1] 10 36 31 > > l_cl_exact > [,1] [,2] [,3] > [1,] 0.44390454 0.2184996 0.2314244 > [2,] 0.04667766 NaN 0.2454760 > [3,] 0.05452433 0.1855618 NaN > [4,] 0.44390454 0.4862702 0.1855618 > [5,] 0.10115053 NaN 0.2184996 > > Thanks for your help. > Anamika > > [[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.