Re: [R] rbinom

2011-12-30 Thread Scott Raynaud
This makes sense.  Guess I should have put a pencil to it.

Further investigation revealed that it is indeed a possibility 
that the relation between x and y is nonlinear:

ax+bx^2+c

where a, b and c are to be determined.  My question is 
how to code this in my simulated data.  I could do 
something like this after appropriately defining beta.
meanpred and varpred:

  x[,2]-rnorm(length,meanpred[2],sqrt(varpred[2]))
  x[,3]-rnorm(length,meanpred[3],sqrt(varpred[3]))
   fixpart-x%*%beta
    binomprob-exp(fixpart)/(1+exp(fixpart))
 data$y-rbinom(n1,1,binomprob)

but I'd need to square my x[,3] values before multiplying 
them by beta.  Can I say:

x[,3]-(rnorm(length,meanpred[3],sqrt(varpred[3])))^2 in
lieu of x[,3]-rnorm(length,meanpred[3],sqrt(varpred[3]))?


- Original Message -
From: peter dalgaard pda...@gmail.com
To: Scott Raynaud scott.rayn...@yahoo.com
Cc: r-help@r-project.org r-help@r-project.org
Sent: Tuesday, December 27, 2011 9:15 AM
Subject: Re: [R] rbinom


On Dec 27, 2011, at 15:47 , Scott Raynaud wrote:

 I have the following code (which I did not write) that generates 
 data based on a logistic model.  I'm only getting a single record 
 with y=1.  It seems implausible that in 50k cases that have a 
 single y=1.  Does that ring alarm bells for anyone else?
  

Not really. As far as I can tell, fixpart is roughly -10.5 (= -1.5 - .25*36), 
so binomprob is around 2.75e-5, which - nonlinearity notwithstanding - suggests 
that the expected number of positives out of 50K is something like 1.4.

To do this more precisely, just compute and print sum(binomprob) in the code 
you gave.

 beta-c(-1.585600,-0.246900)
 betasize-length(beta)
 meanpred-c(0,35.90)
 varpred-c(0,1.00)
 #loop code
 x-matrix(1,length,betasize) #length set to 50k
 #loop code
  x[,2]-rnorm(length,meanpred[2],sqrt(varpred[2])) #length set to 50k
    fixpart-x%*%beta
    binomprob-exp(fixpart)/(1+exp(fixpart))
      data$y-rbinom(n1,1,binomprob)
 #more loop 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.

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd@cbs.dk  Priv: pda...@gmail.com

__
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.


Re: [R] rbinom

2011-12-27 Thread peter dalgaard

On Dec 27, 2011, at 15:47 , Scott Raynaud wrote:

 I have the following code (which I did not write) that generates 
 data based on a logistic model.  I'm only getting a single record 
 with y=1.  It seems implausible that in 50k cases that have a 
 single y=1.  Does that ring alarm bells for anyone else?
  

Not really. As far as I can tell, fixpart is roughly -10.5 (= -1.5 - .25*36), 
so binomprob is around 2.75e-5, which - nonlinearity notwithstanding - suggests 
that the expected number of positives out of 50K is something like 1.4.

To do this more precisely, just compute and print sum(binomprob) in the code 
you gave.

 beta-c(-1.585600,-0.246900)
 betasize-length(beta)
 meanpred-c(0,35.90)
 varpred-c(0,1.00)
 #loop code
 x-matrix(1,length,betasize) #length set to 50k
 #loop code
   x[,2]-rnorm(length,meanpred[2],sqrt(varpred[2])) #length set to 50k
fixpart-x%*%beta
 binomprob-exp(fixpart)/(1+exp(fixpart))
  data$y-rbinom(n1,1,binomprob)
 #more loop 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.

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd@cbs.dk  Priv: pda...@gmail.com

__
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] rbinom

2011-12-27 Thread Scott Raynaud
I have the following code (which I did not write) that generates 
data based on a logistic model.  I'm only getting a single record 
with y=1.  It seems implausible that in 50k cases that have a 
single y=1.  Does that ring alarm bells for anyone else?
 
beta-c(-1.585600,-0.246900)
betasize-length(beta)
meanpred-c(0,35.90)
varpred-c(0,1.00)
#loop code
x-matrix(1,length,betasize) #length set to 50k
#loop code
  x[,2]-rnorm(length,meanpred[2],sqrt(varpred[2])) #length set to 50k
   fixpart-x%*%beta
    binomprob-exp(fixpart)/(1+exp(fixpart))
 data$y-rbinom(n1,1,binomprob)
#more loop 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.


[R] rbinom and probability

2011-02-04 Thread Mcmahon, Kwyatt
Hello compadRes,

I'm developing a script that selects cells over a certain metabolic rate to 
kill them.  A rate between 9 and 12 means that the cells are candidates for 
death.

I'll show you what I mean:

# a would be a vector of cell metabolic rates.
a-c(8, 7, 9, 6, 10, 11, 4, 5, 6)

#now identify which cells will be candidates for death, namely those cells with 
metabolic rates9
b-a[a6]
#This would select all cells with rates 6
rbinom(length(b), size=1, prob=1)

#This will select cells using a an exponential probability distribution
rbinom(length(b), size=1, prob=(seq(0, 1, by=0.1)^2))

I have two questions:

1)  Am I correct in my interpretations?

2)  I'd actually like to have this probability scaled so that cells with 
rates of 9 are unlikely to be selected and those near 12 are highly likely.  
How can I code this?

Thanks in advance,

Wyatt

 sessionInfo()
R version 2.11.1 (2010-05-31)
i386-pc-mingw32

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

loaded via a namespace (and not attached):
[1] tools_2.11.1



K. Wyatt McMahon, Ph.D.
Postdoctoral Research Associate
Department of Internal Medicine
3601 4th St. - Lubbock, TX - 79430
P: 806-743-4072
F: 806-743-3148


[[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.


Re: [R] rbinom and probability

2011-02-04 Thread David Winsemius


On Feb 4, 2011, at 4:24 PM, Mcmahon, Kwyatt wrote:


Hello compadRes,

I'm developing a script that selects cells over a certain  
metabolic rate to kill them.  A rate between 9 and 12 means that the  
cells are candidates for death.


I'll show you what I mean:

# a would be a vector of cell metabolic rates.
a-c(8, 7, 9, 6, 10, 11, 4, 5, 6)

#now identify which cells will be candidates for death, namely those  
cells with metabolic rates9

b-a[a6]
#This would select all cells with rates 6
rbinom(length(b), size=1, prob=1)


No, you have an overly complex way of generating a vector of all 1's  
with length(b). And it's not clear what you mean by select. I see no  
caode that would do any selection or indexing.




#This will select cells using a an exponential probability  
distribution

rbinom(length(b), size=1, prob=(seq(0, 1, by=0.1)^2))


Again no selection appears to be taking place.



I have two questions:

1)  Am I correct in my interpretations?

2)  I'd actually like to have this probability scaled so that  
cells with rates of 9 are unlikely to be selected and those near 12  
are highly likely.  How can I code this?


a[ sample(1:length(a), prob=0.1+(a9) ,replace=TRUE) ]
[1] 10 10 11 10 10 11 10 11 10
 a[sample(1:length(a), prob=0.1+(a9) ,replace=TRUE) ]
[1] 10 10 11 11 11 11  5 10 10

The rates over 9 are 11 times  (= 1.1/0.1) as likely to be selected.  
You can adjust the ratio by altering the 0.1 value.+ with the first  
argument numeric act to coerce the logical vector to 1's and 0's.


--
David.




Thanks in advance,

Wyatt


sessionInfo()

R version 2.11.1 (2010-05-31)
i386-pc-mingw32

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

loaded via a namespace (and not attached):
[1] tools_2.11.1



K. Wyatt McMahon, Ph.D.
Postdoctoral Research Associate
Department of Internal Medicine
3601 4th St. - Lubbock, TX - 79430
P: 806-743-4072
F: 806-743-3148


[[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.


David Winsemius, MD
West Hartford, CT

__
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.


Re: [R] rbinom for a matrix

2008-07-10 Thread Ben Bolker
ACroske Audy3272 at yahoo.com writes:

 
 
 Ben:
 Thanks for the reply. One further question, and this is where my novice
 status at R shows through. The code makes sense, but what would I put it for
 m? Is it the same number for all three (that was my first thought since it
 was the same placeholder for all three). Number of rows in the matrix is 56.
 There are 2576 total cells (observations). Each cell has its own
 probability.
 Thanks for your help!
 
  
  matrix(rbinom(length(m),prob=m,size=1),nrow=nrow(m))
  
or (perhaps marginally more efficiently?)
  
  y - (runif(m)m)
  storage.mode(y) - double
  

  m is your matrix of probabilities; for
example, you might do

 m = as.matrix(read.table(mystuff.txt))

(or import it from ArcGIS, or whatever).

length(m) will be computed automatically (=2576),
as will nrow(m) (=56).

  Erik Iverson's solution should work too, but
I think mine will be much more efficient -- picking
a whole slug of random numbers at once is much
faster than looping and picking them one at a time.

  Ben Bolker

__
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] rbinom for a matrix

2008-07-09 Thread ACroske

I have a large matrix full of probabilities; I would like to convert each
probability to a 1 or a 0 using rbinom.
How can I do this on the entire matrix? The matrix was converted from a
raster ArcMap dataset, so the matrix is essentially a map. Because of this,
I have no column headings.
Thanks! 
-- 
View this message in context: 
http://www.nabble.com/rbinom-for-a-matrix-tp18366867p18366867.html
Sent from the R help mailing list archive at Nabble.com.

__
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.


Re: [R] rbinom for a matrix

2008-07-09 Thread Ben Bolker
ACroske Audy3272 at yahoo.com writes:

 
 
 I have a large matrix full of probabilities; I would like to convert each
 probability to a 1 or a 0 using rbinom.
 How can I do this on the entire matrix? The matrix was converted from a
 raster ArcMap dataset, so the matrix is essentially a map. Because of this,
 I have no column headings.
 Thanks! 

  How about

matrix(rbinom(length(m),prob=m,size=1),nrow=nrow(m))

  or (perhaps marginally more efficiently?)

y - (runif(m)m)
storage.mode(y) - double

  Ben Bolker

__
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.


Re: [R] rbinom for a matrix

2008-07-09 Thread Erik Iverson

Is this what you're looking for?

test - matrix(runif(100, 0, 1), nrow = 20)
nr - nrow(test)
matrix(sapply(test, rbinom, n = 1, size = 1), nrow = nr)

ACroske wrote:

I have a large matrix full of probabilities; I would like to convert each
probability to a 1 or a 0 using rbinom.
How can I do this on the entire matrix? The matrix was converted from a
raster ArcMap dataset, so the matrix is essentially a map. Because of this,
I have no column headings.
Thanks!


__
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.


Re: [R] rbinom for a matrix

2008-07-09 Thread Dylan Beaudette
On Wednesday 09 July 2008, Ben Bolker wrote:
 ACroske Audy3272 at yahoo.com writes:
  I have a large matrix full of probabilities; I would like to convert each
  probability to a 1 or a 0 using rbinom.
  How can I do this on the entire matrix? The matrix was converted from a
  raster ArcMap dataset, so the matrix is essentially a map. Because of
  this, I have no column headings.
  Thanks!

   How about

 matrix(rbinom(length(m),prob=m,size=1),nrow=nrow(m))

   or (perhaps marginally more efficiently?)

 y - (runif(m)m)
 storage.mode(y) - double

   Ben Bolker


Wait a second. Are you trying to convert each probability into the most 
likely '1' or '0' through rounding? The code example above will give you a 
different answer every time you run it. Is that what you are looking for?

Just curious,

Dylan


-- 
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341

__
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.


Re: [R] rbinom for a matrix

2008-07-09 Thread Ben Bolker

-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Dylan Beaudette wrote:
| On Wednesday 09 July 2008, Ben Bolker wrote:
| ACroske Audy3272 at yahoo.com writes:
| I have a large matrix full of probabilities; I would like to convert
each
| probability to a 1 or a 0 using rbinom.
| How can I do this on the entire matrix? The matrix was converted from a
| raster ArcMap dataset, so the matrix is essentially a map. Because of
| this, I have no column headings.
| Thanks!
|   How about
|
| matrix(rbinom(length(m),prob=m,size=1),nrow=nrow(m))
|
|   or (perhaps marginally more efficiently?)
|
| y - (runif(m)m)
| storage.mode(y) - double
|
|   Ben Bolker
|
|
| Wait a second. Are you trying to convert each probability into the most
| likely '1' or '0' through rounding? The code example above will give
you a
| different answer every time you run it. Is that what you are looking for?
|
| Just curious,
|
| Dylan
|

~   I assumed that since the original poster said convert each
probability to a 1 or 0 using rbinom that they did indeed want
a random assignment, rather than rounding ...
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.6 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org

iD8DBQFIdRiZc5UpGjwzenMRAuLkAKCA6ZhjrVsg5RJGYFGwq+6vz2pWxgCfetYk
9V/ZPGi2jfVGV5jx+LGPs8s=
=d8c9
-END PGP SIGNATURE-

__
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.


Re: [R] rbinom for a matrix

2008-07-09 Thread Dylan Beaudette
On Wednesday 09 July 2008, Ben Bolker wrote:
 Dylan Beaudette wrote:
 | On Wednesday 09 July 2008, Ben Bolker wrote:
 | ACroske Audy3272 at yahoo.com writes:
 | I have a large matrix full of probabilities; I would like to convert

 each

 | probability to a 1 or a 0 using rbinom.
 | How can I do this on the entire matrix? The matrix was converted from a
 | raster ArcMap dataset, so the matrix is essentially a map. Because of
 | this, I have no column headings.
 | Thanks!
 |
 |   How about
 |
 | matrix(rbinom(length(m),prob=m,size=1),nrow=nrow(m))
 |
 |   or (perhaps marginally more efficiently?)
 |
 | y - (runif(m)m)
 | storage.mode(y) - double
 |
 |   Ben Bolker
 |
 | Wait a second. Are you trying to convert each probability into the most
 | likely '1' or '0' through rounding? The code example above will give

 you a

 | different answer every time you run it. Is that what you are looking for?
 |
 | Just curious,
 |
 | Dylan

 ~   I assumed that since the original poster said convert each
 probability to a 1 or 0 using rbinom that they did indeed want
 a random assignment, rather than rounding ...

Sorry about not including a reference to whom I was addressing. Meant that 
comment to go to the original poster.

-- 
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341

__
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.


Re: [R] rbinom for a matrix

2008-07-09 Thread ACroske

Ben:
Thanks for the reply. One further question, and this is where my novice
status at R shows through. The code makes sense, but what would I put it for
m? Is it the same number for all three (that was my first thought since it
was the same placeholder for all three). Number of rows in the matrix is 56.
There are 2576 total cells (observations). Each cell has its own
probability.
Thanks for your help!




Ben Bolker wrote:
 
 ACroske Audy3272 at yahoo.com writes:
 
 
 
 I have a large matrix full of probabilities; I would like to convert each
 probability to a 1 or a 0 using rbinom.
 How can I do this on the entire matrix? The matrix was converted from a
 raster ArcMap dataset, so the matrix is essentially a map. Because of
 this,
 I have no column headings.
 Thanks! 
 
   How about
 
 matrix(rbinom(length(m),prob=m,size=1),nrow=nrow(m))
 
   or (perhaps marginally more efficiently?)
 
 y - (runif(m)m)
 storage.mode(y) - double
 
   Ben Bolker
 
 __
 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.
 
 

-- 
View this message in context: 
http://www.nabble.com/rbinom-for-a-matrix-tp18366867p18369250.html
Sent from the R help mailing list archive at Nabble.com.

__
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.


Re: [R] rbinom for a matrix

2008-07-09 Thread ACroske

Yes I do want a random assignment, instead of rounding. (From what I
understand of the rbinom command, it will randomly assign 1 or 0, and the
higher the given probability, the higher the likelihood of a 1... Feel free
to correct me if I'm wrong!)



Ben Bolker wrote:
 
 -BEGIN PGP SIGNED MESSAGE-
 Hash: SHA1
 
 Dylan Beaudette wrote:
 | On Wednesday 09 July 2008, Ben Bolker wrote:
 | ACroske Audy3272 at yahoo.com writes:
 | I have a large matrix full of probabilities; I would like to convert
 each
 | probability to a 1 or a 0 using rbinom.
 | How can I do this on the entire matrix? The matrix was converted from
 a
 | raster ArcMap dataset, so the matrix is essentially a map. Because of
 | this, I have no column headings.
 | Thanks!
 |   How about
 |
 | matrix(rbinom(length(m),prob=m,size=1),nrow=nrow(m))
 |
 |   or (perhaps marginally more efficiently?)
 |
 | y - (runif(m)m)
 | storage.mode(y) - double
 |
 |   Ben Bolker
 |
 |
 | Wait a second. Are you trying to convert each probability into the most
 | likely '1' or '0' through rounding? The code example above will give
 you a
 | different answer every time you run it. Is that what you are looking
 for?
 |
 | Just curious,
 |
 | Dylan
 |
 
 ~   I assumed that since the original poster said convert each
 probability to a 1 or 0 using rbinom that they did indeed want
 a random assignment, rather than rounding ...
 -BEGIN PGP SIGNATURE-
 Version: GnuPG v1.4.6 (GNU/Linux)
 Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org
 
 iD8DBQFIdRiZc5UpGjwzenMRAuLkAKCA6ZhjrVsg5RJGYFGwq+6vz2pWxgCfetYk
 9V/ZPGi2jfVGV5jx+LGPs8s=
 =d8c9
 -END PGP SIGNATURE-
 
 __
 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.
 
 

-- 
View this message in context: 
http://www.nabble.com/rbinom-for-a-matrix-tp18366867p18370010.html
Sent from the R help mailing list archive at Nabble.com.

__
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] rbinom not using probability of success right

2008-05-29 Thread Kyle Matoba

Message: 24
Date: Wed, 28 May 2008 05:53:26 -0700 (PDT)
From: Philip Twumasi-Ankrah [EMAIL PROTECTED]
Subject: [R] rbinom not using probability of success right
To: r-help@r-project.org
Message-ID: [EMAIL PROTECTED]
Content-Type: text/plain

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 probability of success to be 0.154. Not very good.

Is there another way beyond using sample and rep together?

I understand you correctly you want there to be exactly 75 ones.  If this 
is what you are trying to do then using pseudorandom variables is the 
incorrect way of going about it.  Your suggestion of sample(c(rep
(0,545),rep(1,75))) seems to me to be the best way of going about it since 
conceptually this is what you are doing: taking permutation of a fixed set 
of numbers.

Best,

Kyle

__
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] rbinom not using probability of success right

2008-05-28 Thread Philip Twumasi-Ankrah
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 probability of success to be 0.154. Not very good.

Is there another way beyond using sample and rep together?


A Smile costs Nothing  

 But Rewards Everything

Happiness is not perfected until it is shared
  -Jane Porter  
  


   
[[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.


Re: [R] rbinom not using probability of success right

2008-05-28 Thread Prof Brian Ripley
You asked for each of 500 to be included with probability 0.15, not for 
15% of 500.  If you want the latter, use sample, e.g.


sample(c(rep(1,75), rep(0,425)))

And to see if your 77 is reasonable for binomial sampling:


binom.test(77, 500, 0.15)


Exact binomial test

data:  77 and 500
number of successes = 77, number of trials = 500, p-value = 0.8022
alternative hypothesis: true probability of success is not equal to 0.15
95 percent confidence interval:
 0.1234860 0.1886725
sample estimates:
probability of success
 0.154

so it certainly is.

On Wed, 28 May 2008, Philip Twumasi-Ankrah wrote:


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 probability of success to be 0.154. Not very good.

Is there another way beyond using sample and rep together?


[...]

--
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

__
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.


Re: [R] rbinom not using probability of success right

2008-05-28 Thread Ted Harding
On 28-May-08 12:53:26, Philip Twumasi-Ankrah wrote:
 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 probability of success to be 0.154. Not
 very good.

The difference (77 - 75 =2) is well within the likely sampling
variation when 500 values are sampled independently with
P(1)=0.15:

The standard deviation of the resulting number of 1s is
sqrt(500*0.15*0.85) =  7.98, so the difference of 2 is only 1/4 of
a standard deviation, hence very likely to be equalled or exceeded.

Your chance of getting exactly 75 by this method is quite small:

  dbinom(75,500,0.15)
  [1] 0.04990852

and your chance of being 2 or more off your target is

  1 - sum(dbinom((74:76),500,0.15))
  [1] 0.8510483

 Is there another way beyond using sample and rep together?

It looks as though you are seeking to obtain exactly 75 1s,
randomly situated, the rest being 0s, so in effect you do need
to do something on the lines of sample and rep. Hence,
something like

  status - rep(0,500)
  status[sample((1:500),75,replace=FALSE)] - 1

Hoping this helps,
Ted.


E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 094 0861
Date: 28-May-08   Time: 14:19:24
-- XFMail --

__
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.


Re: [R] rbinom not using probability of success right

2008-05-28 Thread Charles Annis, P.E.
Do it again.  What did you get this time?  Then do it another time.  Do you
see what is happening?

Charles Annis, P.E.

[EMAIL PROTECTED]
phone: 561-352-9699
eFax:  614-455-3265
http://www.StatisticalEngineering.com
 

-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Philip Twumasi-Ankrah
Sent: Wednesday, May 28, 2008 8:53 AM
To: r-help@r-project.org
Subject: [R] 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 probability of success to be 0.154. Not very good.

Is there another way beyond using sample and rep together?


A Smile costs Nothing  

 But Rewards Everything

Happiness is not perfected until it is shared
  -Jane Porter



   
[[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.


Re: [R] rbinom : Does randomness preclude precision?

2008-05-28 Thread Philip Twumasi-Ankrah
Teds reply is a bit comforting and as indicated in my post, I am resorting to 
using sample but as an academic issue, does randomness preclude precision? 

Randomness should be in the sequence of zeros and ones and how they are 
simulated at each iteration of the process but not in the eventual nature of 
the distribution. 

I mean if I simulated a Normal (0, 1) and got a Normal(1.5, 2) these would be 
very different distributions. It is the same with simulating a Binomial(1, 
p=0.15) and getting Binomial(1, 0.154)

[EMAIL PROTECTED] wrote: On 28-May-08 12:53:26, Philip Twumasi-Ankrah wrote:
 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 probability of success to be 0.154. Not
 very good.

The difference (77 - 75 =2) is well within the likely sampling
variation when 500 values are sampled independently with
P(1)=0.15:

The standard deviation of the resulting number of 1s is
sqrt(500*0.15*0.85) =  7.98, so the difference of 2 is only 1/4 of
a standard deviation, hence very likely to be equalled or exceeded.

Your chance of getting exactly 75 by this method is quite small:

  dbinom(75,500,0.15)
  [1] 0.04990852

and your chance of being 2 or more off your target is

  1 - sum(dbinom((74:76),500,0.15))
  [1] 0.8510483

 Is there another way beyond using sample and rep together?

It looks as though you are seeking to obtain exactly 75 1s,
randomly situated, the rest being 0s, so in effect you do need
to do something on the lines of sample and rep. Hence,
something like

  status - rep(0,500)
  status[sample((1:500),75,replace=FALSE)] - 1

Hoping this helps,
Ted.


E-Mail: (Ted Harding) 
Fax-to-email: +44 (0)870 094 0861
Date: 28-May-08   Time: 14:19:24
-- XFMail --



A Smile costs Nothing  

 But Rewards Everything

Happiness is not perfected until it is shared
  -Jane Porter  
  


   
[[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.


Re: [R] rbinom : Does randomness preclude precision?

2008-05-28 Thread Charles Annis, P.E.
What do you mean by ... *eventual* nature of the distribution?  If you
simulated 100 samples, would you expect to see 1.5 successes?  Or 1?  Or 2?
How many, in your thinking, is eventual?



Charles Annis, P.E.

[EMAIL PROTECTED]
phone: 561-352-9699
eFax:  614-455-3265
http://www.StatisticalEngineering.com
 

-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Philip Twumasi-Ankrah
Sent: Wednesday, May 28, 2008 9:52 AM
To: [EMAIL PROTECTED]
Cc: r-help@r-project.org
Subject: Re: [R] rbinom : Does randomness preclude precision?

Teds reply is a bit comforting and as indicated in my post, I am resorting
to using sample but as an academic issue, does randomness preclude
precision? 

Randomness should be in the sequence of zeros and ones and how they are
simulated at each iteration of the process but not in the eventual nature of
the distribution. 

I mean if I simulated a Normal (0, 1) and got a Normal(1.5, 2) these would
be very different distributions. It is the same with simulating a
Binomial(1, p=0.15) and getting Binomial(1, 0.154)

[EMAIL PROTECTED] wrote: On 28-May-08 12:53:26, Philip
Twumasi-Ankrah wrote:
 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 probability of success to be 0.154. Not
 very good.

The difference (77 - 75 =2) is well within the likely sampling
variation when 500 values are sampled independently with
P(1)=0.15:

The standard deviation of the resulting number of 1s is
sqrt(500*0.15*0.85) =  7.98, so the difference of 2 is only 1/4 of
a standard deviation, hence very likely to be equalled or exceeded.

Your chance of getting exactly 75 by this method is quite small:

  dbinom(75,500,0.15)
  [1] 0.04990852

and your chance of being 2 or more off your target is

  1 - sum(dbinom((74:76),500,0.15))
  [1] 0.8510483

 Is there another way beyond using sample and rep together?

It looks as though you are seeking to obtain exactly 75 1s,
randomly situated, the rest being 0s, so in effect you do need
to do something on the lines of sample and rep. Hence,
something like

  status - rep(0,500)
  status[sample((1:500),75,replace=FALSE)] - 1

Hoping this helps,
Ted.


E-Mail: (Ted Harding) 
Fax-to-email: +44 (0)870 094 0861
Date: 28-May-08   Time: 14:19:24
-- XFMail --



A Smile costs Nothing  

 But Rewards Everything

Happiness is not perfected until it is shared
  -Jane Porter



   
[[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.


Re: [R] rbinom : Does randomness preclude precision?

2008-05-28 Thread Charles Annis, P.E.
I think I see the rub:  You would like to see the distribution of a sample
be identical to the distribution from which it was sampled.  But if it is
random then that can happen only in the long run, not on every sample.  That
is why samples from a normal density are *not* themselves normal - they're
t.  When the sample size is large enough the differences between a random
sample's density and its parent density become vanishingly small.  Thus the
differences you observe from repeated random samples from the binomial.
Repeated sampling produces slightly different numbers of successes.  How
could it be otherwise?

Charles Annis, P.E.

[EMAIL PROTECTED]
phone: 561-352-9699
eFax:  614-455-3265
http://www.StatisticalEngineering.com
 

From: Philip Twumasi-Ankrah [mailto:[EMAIL PROTECTED] 
Sent: Wednesday, May 28, 2008 10:36 AM
To: [EMAIL PROTECTED]
Subject: RE: [R] rbinom : Does randomness preclude precision?

Charles,
When you simulate data from a distribution, what you effect are doing is
generating a sequence of values that would correspond to that distribution.
So you can generate 1000 values from a normal distribution and expect that
when you check on the distribution of your sample (what you do with your
qqnorm or Q-Q plot), it should be a close fit with the theoretical
distribution with the assigned parameter values. It will be difficult to
explain why a simulated data may be different from the distribution it is
was generated from . I think you can not blame it on randomness. 

I hope you understand what I am trying to determine.

Charles Annis, P.E. [EMAIL PROTECTED] wrote:
What do you mean by ... *eventual* nature of the distribution? If you
simulated 100 samples, would you expect to see 1.5 successes? Or 1? Or 2?
How many, in your thinking, is eventual?



Charles Annis, P.E.

[EMAIL PROTECTED]
phone: 561-352-9699
eFax: 614-455-3265
http://www.StatisticalEngineering.com


-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Philip Twumasi-Ankrah
Sent: Wednesday, May 28, 2008 9:52 AM
To: [EMAIL PROTECTED]
Cc: r-help@r-project.org
Subject: Re: [R] rbinom : Does randomness preclude precision?

Teds reply is a bit comforting and as indicated in my post, I am resorting
to using sample but as an academic issue, does randomness preclude
precision? 

Randomness should be in the sequence of zeros and ones and how they are
simulated at each iteration of the process but not in the eventual nature of
the distribution. 

I mean if I simulated a Normal (0, 1) and got a Normal(1.5, 2) these would
be very different distributions. It is the same with simulating a
Binomial(1, p=0.15) and getting Binomial(1, 0.154)

[EMAIL PROTECTED] wrote: On 28-May-08 12:53:26, Philip
Twumasi-Ankrah wrote:
 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 probability of success to be 0.154. Not
 very good.

The difference (77 - 75 =2) is well within the likely sampling
variation when 500 values are sampled independently with
P(1)=0.15:

The standard deviation of the resulting number of 1s is
sqrt(500*0.15*0.85) = 7.98, so the difference of 2 is only 1/4 of
a standard deviation, hence very likely to be equalled or exceeded.

Your chance of getting exactly 75 by this method is quite small:

dbinom(75,500,0.15)
[1] 0.04990852

and your chance of being 2 or more off your target is

1 - sum(dbinom((74:76),500,0.15))
[1] 0.8510483

 Is there another way beyond using sample and rep together?

It looks as though you are seeking to obtain exactly 75 1s,
randomly situated, the rest being 0s, so in effect you do need
to do something on the lines of sample and rep. Hence,
something like

status - rep(0,500)
status[sample((1:500),75,replace=FALSE)] - 1

Hoping this helps,
Ted.


E-Mail: (Ted Harding) 
Fax-to-email: +44 (0)870 094 0861
Date: 28-May-08 Time: 14:19:24
-- XFMail --



A Smile costs Nothing 

But Rewards Everything

Happiness is not perfected until it is shared
-Jane Porter




[[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.


A Smile costs Nothing  
But Rewards Everything

Happiness is not perfected until it is shared
  -Jane Porter  

  

__
R-help@r-project.org mailing

Re: [R] rbinom : Does randomness preclude precision?

2008-05-28 Thread Daniel Nordlund
 -Original Message-
 From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf
 Of Philip Twumasi-Ankrah
 Sent: Wednesday, May 28, 2008 6:52 AM
 To: [EMAIL PROTECTED]
 Cc: r-help@r-project.org
 Subject: Re: [R] rbinom : Does randomness preclude precision?
 
 Teds reply is a bit comforting and as indicated in my post, I am resorting to 
 using
 sample but as an academic issue, does randomness preclude precision?

I would say yes, at least in the manner you seem to be thinking of it.  Think 
about it for a minute.  Say you use the following code

status - sample(rep(c(0,1),c(425,75)))

You will get your 75 ones randomly distributed throughout the vector status.  
What would you expect the value to be for

sum(status[1:100])

I suggest that you will see the same kind of variation in the sum as you would 
see for 

sum(rbinom(100, 1, p=.15))

So, you can randomly arrange a certain percentage of ones in a sequence, but 
you should not expect to see that exact percentage in any subset of the 
sequence.

One other example, if you flip a fair coin 100 times, do you expect to get 
exactly 50 heads?

Hope this is helpful,

Dan

Daniel Nordlund
Bothell, WA USA


 
 Randomness should be in the sequence of zeros and ones and how they are 
 simulated
 at each iteration of the process but not in the eventual nature of the 
 distribution.
 
 I mean if I simulated a Normal (0, 1) and got a Normal(1.5, 2) these would be 
 very
 different distributions. It is the same with simulating a Binomial(1, p=0.15) 
 and getting
 Binomial(1, 0.154)
 
 [EMAIL PROTECTED] wrote: On 28-May-08 12:53:26, Philip Twumasi-
 Ankrah wrote:
  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 probability of success to be 0.154. Not
  very good.
 
 The difference (77 - 75 =2) is well within the likely sampling
 variation when 500 values are sampled independently with
 P(1)=0.15:
 
 The standard deviation of the resulting number of 1s is
 sqrt(500*0.15*0.85) =  7.98, so the difference of 2 is only 1/4 of
 a standard deviation, hence very likely to be equalled or exceeded.
 
 Your chance of getting exactly 75 by this method is quite small:
 
   dbinom(75,500,0.15)
   [1] 0.04990852
 
 and your chance of being 2 or more off your target is
 
   1 - sum(dbinom((74:76),500,0.15))
   [1] 0.8510483
 
  Is there another way beyond using sample and rep together?
 
 It looks as though you are seeking to obtain exactly 75 1s,
 randomly situated, the rest being 0s, so in effect you do need
 to do something on the lines of sample and rep. Hence,
 something like
 
   status - rep(0,500)
   status[sample((1:500),75,replace=FALSE)] - 1
 
 Hoping this helps,
 Ted.

__
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.


Re: [R] rbinom with computed probability

2007-11-24 Thread Berwin A Turlach
G'day Sigalit,

On Fri, 23 Nov 2007 11:25:30 +0200
sigalit mangut-leiba [EMAIL PROTECTED] wrote:

 Hello,
 I have a loop with probability computed from a logistic model like
 this: for (i in 1:300){
 
 p[i]-exp(-0.834+0.002*x[i]+0.023*z[i])/(1+exp(-0.834+0.002*x[i]+0.023
 +z[i]))
 
 x and z generated from normal distribution.
 I get 300 different probabilities And I want to generate variables
 from bernulli distribution
 with P for every observation:
 
 T[i]-rbinom(1,1,p[i])
 But i get missing values for T.
 What I'm doing wrong?

I guess the problem is that in the numerator you have 0.023*z[i] but
0.023+z[i] in the denominator.  Thus, some p[i] can be outside of
[0,1] which would produce NAs in T.

But why a for loop? This code is readily vectorised:

p - exp(-0.834+0.002*x+0.023*z)/(1+exp(-0.834+0.002*x+0.023*z))

HTH.

Cheers,

Berwin

=== Full address =
Berwin A TurlachTel.: +65 6515 4416 (secr)
Dept of Statistics and Applied Probability+65 6515 6650 (self)
Faculty of Science  FAX : +65 6872 3919   
National University of Singapore
6 Science Drive 2, Blk S16, Level 7  e-mail: [EMAIL PROTECTED]
Singapore 117546http://www.stat.nus.edu.sg/~statba

__
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] rbinom with computed probability

2007-11-23 Thread sigalit mangut-leiba
Hello,
I have a loop with probability computed from a logistic model like this:
for (i in 1:300){

p[i]-exp(-0.834+0.002*x[i]+0.023*z[i])/(1+exp(-0.834+0.002*x[i]+0.023
+z[i]))

x and z generated from normal distribution.
I get 300 different probabilities And I want to generate variables from
bernulli distribution
with P for every observation:

T[i]-rbinom(1,1,p[i])
But i get missing values for T.
What I'm doing wrong?
Thank you,
Sigalit.

[[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.


Re: [R] rbinom with computed probability

2007-11-23 Thread Ingmar Visser
I don't run into problems doing this:

x=rnorm(300)
z=rnorm(300)
for (i in 1:300){
  p[i]-exp(-0.834+0.002*x[i]+0.023*z[i])/(1+exp(-0.834+0.002*x[i]+0.023
+z[i])); T[i]-rbinom(1,1,p[i])
}

hth, Ingmar

On 23 Nov 2007, at 10:25, sigalit mangut-leiba wrote:

 Hello,
 I have a loop with probability computed from a logistic model like  
 this:
 for (i in 1:300){

 p[i]-exp(-0.834+0.002*x[i]+0.023*z[i])/(1+exp(-0.834+0.002*x[i]+0.023
 +z[i]))

 x and z generated from normal distribution.
 I get 300 different probabilities And I want to generate variables  
 from
 bernulli distribution
 with P for every observation:

 T[i]-rbinom(1,1,p[i])
 But i get missing values for T.
 What I'm doing wrong?
 Thank you,
 Sigalit.

   [[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.

Ingmar Visser
Department of Psychology, University of Amsterdam
Roetersstraat 15
1018 WB Amsterdam
The Netherlands
t: +31-20-5256723



[[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.