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-boun...@r-project.org [mailto:r-help-boun...@r-project.org] > On Behalf Of Anamika Chaudhuri > Sent: Freitag, 28. Juni 2013 04:06 > To: r-help@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@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.