search for: rbinom

Displaying 20 results from an estimated 118 matches for "rbinom".

2008 Sep 09
naive variance in GEE
Hi, The standard error from logistic regression is slightly different from the naive SE from GEE under independence working correlation structure. Shouldn't they be identical? Anyone has insight about this? Thanks, Qiong a<-rbinom(1000,1) b<-rbinom(1000,2,0.1) c<-rbinom(1000,10,0.5) summary(gee(a~b, id=c,family="binomial",corstr="independence"))$coef summary(glm(a~b,family="binomial"))
2012 Jan 03
nls and rbinom function: step factor 0.000488281 reduced below 'minFactor' of 0.000976562
I  am trying to learn nls using a simple simulation. I assumed that the binomial prob varies linearly as 0.2 + 0.3*x in  x {0,1}, and the objective is to recover the known parameters a=0.2, b=0.3 frame d has 1000 rows... d$x<-runif(0,1)               d$y<-rbinom(1000,1,0.2+0.3*d$x)  table(d$y,cut(d$x,breaks=5));   (-0.000585,0.199] (0.199,0.399] (0.399,0.599] (0.599,0.799] (0.799,0.999]   0               154           149           130           122           114   1                34            48            71            76           102  z <- nls(...
2011 Dec 14
hclust and ggplot2
...get an error when trying to use ggplot; Error: ggplot2 doesn't know how to deal with data of class phylo. Regular plot works fine but I can't get ggplot2 to work. see code below.... rows=100 columns=100 #create matrix all=matrix(nrow=rows, ncol=columns) #initialize first column all[,1]=rbinom(rows,1,.5) #set probability probv=.9 for (j in 2:columns) { for (i in 1:rows) { #set the probability based on the values in the columns to the left if (all[i,j-1]==1) {all[i,j]=rbinom(1,1,1-probv)} else {all[i,j]=rbinom(1,1,probv)} } } #make...
2008 May 28
"rbinom" not using probability of success right
I am trying to simulate a series of ones and zeros (1 or 0) and I am using "rbinom" but realizing that the number of successes expected is not accurate. Any advice out there. This is the example: N<-500 status<-rbinom(N, 1, prob = 0.15) count<-sum(status) 15 percent of 500 should be 75 but what I obtain from the "count" variable is 77 that gives the pro...
2013 May 21
problem with "transform" and "get" functions
...0's of variables). However, I get the following error message: " Error: unexpected '=' in: "{ test0 <- transform(test, (get(predictor.names[i])) =" > } Error: unexpected '}' in "}" > " Here is my code: n<-200 set.seed(111) X1 <- rbinom(n,1,0.4) X2 <- rbinom(n,1,0.5) X3 <- rbinom(n,1,0.5) X4 <- rbinom(n,1,(0.5 + 0.2*X1 - 0.3*X2 + 0.005*X3)) Y <- rnorm(n,(3 + X1 + 0.3*X2 + 0.1*X3 + 0.05*X4),.4) data.frame(Y,X1,X2,X3,X4) -> test var.names <- colnames(test) for (i in (1:5)) { test0 <- transform(test, (get(var.n...
2005 Feb 28
A problem about outer()
...bout function outer() that I can't understand. Just see the following example. The two NaNs are due to 0/0, but I can't figure out the cause of the last two errors. I wonder if some one can explain this for me. ___________________________________________________________________ > sx=rbinom(10,1,0.5);ot=rbinom(10,1,0.5);ag <- rbinom(10,100,0.3);ho <- rbinom(10,100,0.5) > dp <- function(s,a,h)sum((sx==s)&(ag==a)&(ho==h)&(ot==1))/sum((sx==s)&(ag==a)&(ho==h)) > (function(x,y)dp(1,x,y))(2,3) [1] NaN > (function(x,y)dp(0,x,y))(27,52) [1] NaN > d...
2012 Jan 30
mgcv bam() with grouped binomial data
...nction in the R mgcv package for a large set of grouped binary data. However, I have found that this function does not take data in the format of cbind(numerator, denominator) on the left hand side of the formula. As an example, consider the following dat1 <- data.frame(id=rep(1:6, each=3), num=rbinom(18, size=10, prob=0.8), den=rbinom(18, size=5, prob=0.5), x=rnorm(18)) m1.1 <- gam(cbind(num, den) ~ x+s(id, bs="re"), data=dat1, family=binomial) m2.1 <- bam(cbind(num, den) ~x+ s(id, bs="re"), data=dat1, family=binomial) Running the above results in > dat1 <- dat...
2007 Jul 06
How does the r-distribution function work
I am trying to understand what rbinom function does. Here is some sample code. Are both the invocations of bfunc effectively doing the same or I am missing the point? Thanks, Pieter bfunc <- function(n1,p1,sims) { c<-rbinom(sims,n1,p1) c } a=c() b=c() p1=.5 for (i in 1:10000){ a[i]=bfunc(30,p1,1) } b=bfunc(30,p1,10000)
2006 Feb 06
Evaluate output after each rep()
Hi R-Help, I'm trying a develop a test simulation where i evaluate the probability of not getting a value of 100 from the function rbinom(6000, 200, .5) [indeed, a very small probability]. At the end of each rep, I would like to evaluate the output, continue with the loop if the output contains the value 100, stop if the output lacks a 100. How do I get R to evaluate the output after each rep? >sim <- function(nn){ &gt...
2013 May 23
sample(c(0, 1)...) vs. rbinom
...than mathematics, to get insight into statistics, and R is the mandated tool.?? A student in the class recently inquired about different approaches to sampling from a binomial distribution.? I've appended some code that exhibits the idea, the gist of which is that using sample(c(0, 1), ...) and rbinom(...) should give equivalent results. The surprising (to me) result is that the two approaches DO give the same result, EXCEPT when the probability is exactly 0.5.? See Appendix A for the code and Appendix B for the output.? I don't think this issue is system-dependent, but I've put my sess...
1998 Feb 27
R-beta: is there a way to get rid of loop?
...However, for each ns I want an associated binomial variate with p=phit. Is there a better way to do this than by looping?<-function(nhit,nmiss,nfa,ncr) { #monte carlo simulated 95% CI nsignal<-nhit+nmiss nnoise<-nfa+ncr phit<-nhit/nsignal pfa<-nfa/nnoise ns<-rbinom(2000,nsignal,.5) nn<-nsignal+nnoise-ns nh<-NULL nf<-NULL for(i in 1:2000) { nh<-c(nh,rbinom(1,ns[i],phit)) nf<-c(nf,rbinom(1,nn[i],pfa)) } g<-Gsens(nh,ns-nh,nf,nn-nf) quantile(g,c(.025,.975), na.rm=TRUE) } Thanks very much for any help! Bill Simpson -.-.-.-.-.-.-.-.-.-.-.-....
2001 Apr 08
Random Number Testing...
Hi, Suppose I want to test a set of random numbers (using the Binormal random number generator), input the following command to generate a set of random numbers: test <- rbinom( 10000000, 1, 0.5 ) How can I, after obtaining the result, export the vector "test" into an text file? Then, suppose I want to do something like: > x <- 1e10 > x [1] 1e+10 > test <- rbinom( x, 1, 0.5 ) Error in rbinom(n, size, prob) : invalid argu...
2006 Oct 09
testing for error
...result)<1 #and so choose not use the result but if the procedure does not produce an error, then "if grep('Error', result)<1" now produces the result logical(0) and the loop then fails with the message set.seed(1) cumulator=rep(0,100) for (i in 1:100){ y1=rnorm(100) x0=rbinom(100,1,0.02) x1=rbinom(100,1,0.5) x2=rbinom(100,1,0.5) x1[x0]=NA dat=data.frame(y1,x1) result=try(lm(y1~x1,, data=dat),T); print(result) x1=x2; dat2=data.frame(x1) if (grep('Error',result)<1) cumulator=cumulator+predict(result,x2) } The above runs and rejects the 'r...
2008 Dec 16
simulate binary markov chain
...ulating a large binary sequence (length >10 million) using a (1st order) markov model with known (2x2) transition matrix. It needs to be reasonably fast. I have tried the following; mc<-function(sq,P){ s<-c() x<-row.names(P) n<-length(sq) p1<-sum(sq)/n s[1] <- rbinom(1,1,p1); for ( i in 2:n){ s[i] <- rbinom( 1, 1, P[s[i-1]+1] ) } return(s) } P<-c(0.63,0.27) x<-rbinom(500,1,0.5) new<-mc(x,P) thanks in advance! Chris
2010 Sep 24
Standard Error for difference in predicted probabilities
Is there a way to estimate the standard error for the difference in predicted probabilities obtained from a logistic regression model? For example, this code gives the difference for the predicted probability of when x2==1 vs. when x2==0, holding x1 constant at its mean: y=rbinom(100,1,.4) x1=rnorm(100, 3, 2) x2=rbinom(100, 1, .7) mod=glm(y ~ x1 + x2, family=binomial) pred=predict(mod, newdata=data.frame(cbind(x1=rep(mean(x1), 100), x2)), type="response") diff=unique(pred)[1]-unique(pred)[2] diff I know that predict() will output SE's for each predicted val...
2015 Jan 22
Simulación de modelo logit con interacción
...racias Emilio logisticsimulation <- function(n){ dat <- data.frame(x1=sample(0:1, n,replace=TRUE), x2=sample(0:1, n,replace=TRUE)) odds <- exp(-1 - 4 * dat$x1 + 7*dat$x2 - 1 *dat$x1* dat$x2 ) pr <- odds/(1+odds) res <- replicate(100, { dat$y <- rbinom(n,1,pr) coef(glm(y ~ x1*x2, data = dat, family = binomial())) }) t(res) } res <- logisticsimulation(100) apply(res,2,median) ## (Intercept) x1 x2 x1:x2 ## -1.0986123 -18.4674562 20.4823593 -0.0512933 Deberían salir -1, -4, 7, 1
2008 Jul 14
how to correlate nominal variables?
Dear R-Users, I need functions to calculate Yule's Y or Cram?rs Index, in order to correlate variables that are nominally scaled? Am I wrong? Are such functions existing? Sincerely, Timo
2012 Aug 27
How to generate a matrix of Beta or Binomial distribution
...ectly generate X without any loops. I tried the following code X<- rbeta(2,shape1,shape2) but it looks not right. So I was wondering if there is an R code doing this. For Binomial distribution, I have a similar question. I hope to generate a 10000 x n matrix Y whose i th row is a sample from rbinom(n,2,p[i]), where p is a vector of binomial probabilities. We can also do it by loops for(i in 1:10000) { Y[i,]=rbinom(n,2,p[i]) } But it would be nice to do it without using any loops. Is this possible in R? Many thanks for great help and suggestions! Jeff -- View this message in contex...
2012 Sep 27
Generating an autocorrelated binary variable
Hi R-fellows, I am trying to simulate a multivariate correlated sample via the Gaussian copula method. One variable is a binary variable, that should be autocorrelated. The autocorrelation should be rho = 0.2. Furthermore, the overall probability to get either outcome of the binary variable should be 0.5. Below you can see the R code (I use for simplicity a diagonal matrix in rmvnorm even if it
2012 Oct 25
How to extract auc, specificity and sensitivity
I am running my code in a loop and it does not work but when I run it outside the loop I get the values I want. n <- 1000; # Sample size fitglm <- function(sigma,tau){ x <- rnorm(n,0,sigma) intercept <- 0 beta <- 0 ystar <- intercept+beta*x z <- rbinom(n,1,plogis(ystar)) xerr <- x + rnorm(n,0,tau) model<-glm(z ~ xerr, family=binomial(logit)) int<-coef(model)[1] slope<-coef(model)[2] pred<-predict(model) result<-ifelse(pred>.5,1,0) accuracy<-length(which(result==z))/length(z) accuracy...