R help - Jun 2013 - Add constraints to rbinom

Hi:

I am trying to generate Beta-Binomial random variables and then
calculate Binomial exact confidence intervals to the rate..I was
wondering if I needed to add a constraint such that x<=n to it, how do
I go about it. The reason is I am getting data where x>n which is
giving a rate>1. Heres my code:

set.seed(111)
k<-63
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(k,min,max))
for(i in 1:k)
x<-cbind(x,rbinom(300,n,p[i]))
x<-t(x)
rate<-t(t(x)/n)
se_rate<-sqrt(rate*(1-rate)/n)



# Exact Confidence Interval

l_cl_exact<-qbeta(.025,x,n-x+1)
u_cl_exact<-qbeta(.975,x+1,n-x)

for (i in 1:63){

for (j in 1:300)
{

if (x[i,j]==0)
{
l_cl_exact[i,j]<-0
u_cl_exact[i,j]<-u_cl_exact[i,j]
}
else if (x[i,j]==n[i])
{
l_cl_exact[i,j]<-l_cl_exact[i,j]
u_cl_exact[i,j]<-1
}
else
l_cl_exact[i,j]<-l_cl_exact[i,j]
u_cl_exact[i,j]<-u_cl_exact[i,j]

#print(c(i,j))

}
}


Really appreciate any help.
Thanks
Anamika

You don't need a constraint (rbinom won't give x>n), but you need to
make sure you are using the n you want to use: try
 x <- cbind(x,rbinom(300,n[i],p[i]))	# mind the "[i]" after the n
at the respective line. Furthermore, you need to remove one transformation of x
to make sure you divide by the right n, try
rate <- t(x/n)

I think that should do the trick. I didn't check the remainder of the code,
though.

HTH, Michael

> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at
r-project.org]
> On Behalf Of Anamika Chaudhuri
> Sent: Freitag, 28. Juni 2013 04:06
> To: r-help at r-project.org
> Subject: [R] Add constraints to rbinom
> 
> Hi:
> 
> I am trying to generate Beta-Binomial random variables and then calculate
> Binomial exact confidence intervals to the rate..I was wondering if I
needed
> to add a constraint such that x<=n to it, how do I go about it. The
reason is I
> am getting data where x>n which is giving a rate>1. Heres my code:
> 
> set.seed(111)
> k<-63
> 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(k,min,max))
> for(i in 1:k)
> x<-cbind(x,rbinom(300,n,p[i]))
> x<-t(x)
> rate<-t(t(x)/n)
> se_rate<-sqrt(rate*(1-rate)/n)
> 
> 
> 
> # Exact Confidence Interval
> 
> l_cl_exact<-qbeta(.025,x,n-x+1)
> u_cl_exact<-qbeta(.975,x+1,n-x)
> 
> for (i in 1:63){
> 
> for (j in 1:300)
> {
> 
> if (x[i,j]==0)
> {
> l_cl_exact[i,j]<-0
> u_cl_exact[i,j]<-u_cl_exact[i,j]
> }
> else if (x[i,j]==n[i])
> {
> l_cl_exact[i,j]<-l_cl_exact[i,j]
> u_cl_exact[i,j]<-1
> }
> else
> l_cl_exact[i,j]<-l_cl_exact[i,j]
> u_cl_exact[i,j]<-u_cl_exact[i,j]
> 
> #print(c(i,j))
> 
> }
> }
> 
> 
> Really appreciate any help.
> Thanks
> Anamika
> 
> ______________________________________________
> 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.