[R] qbeta function in R

[William Dunlap]

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.