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-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
> Behalf
> Of Anamika Chaudhuri
> Sent: Friday, March 09, 2012 11:48 AM
> To: r-help@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@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.
______________________________________________
R-help@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.
