Hi, there
Below is my code to one Homework question. I couldn't come up with the
reasonable answer.
could you please help me to figure out what is the problem with my code?
thank you
Question is Coding P{X=j} =(1/2)^(j+1) + (1/2) *2^(j-1)/3^j
my code is
sim <- function(n.gen){
urandom <- runif(n.gen)
sim.vector <- rep(0,n.gen)
for(j in 1:n.gen){
i <- 1
p <- 5/12
F <- p
while(urandom[j] >= F){
p <- p*((1/2)^(i+1)+1/3*(2/3)^i)/((1/2)^i+(1/2)*(2/3)^i)
F <- F+p
i<-i+1
}
sim.vector[j] <- i
}
# output
sim.vector
}
result is
1 2 3 4 5 6 7 8 11
0.37 0.22 0.16 0.13 0.05 0.02 0.03 0.01 0.01
always, there are some numbers missing, it should be continuous.
why 9 and 10 are missing
thank you
sophia
_________________________________________________________________
[[elided Hotmail spam]]
[[alternative HTML version deleted]]
Ssophia wrote:> > Hi, there > > > > Below is my code to one Homework question. I couldn't come up with the reasonable answer.^^^^^^^^ PLEASE do read the posting guide http://www.R-project.org/posting-guide.html> could you please help me to figure out what is the problem with my code? > thank you > > > Question is Coding P{X=j} =(1/2)^(j+1) + (1/2) *2^(j-1)/3^j > my code is > sim <- function(n.gen){ > urandom <- runif(n.gen) > sim.vector <- rep(0,n.gen) > for(j in 1:n.gen){ > i <- 1 > p <- 5/12 > F <- p > while(urandom[j] >= F){ > p <- p*((1/2)^(i+1)+1/3*(2/3)^i)/((1/2)^i+(1/2)*(2/3)^i) > F <- F+p > i<-i+1 > } > sim.vector[j] <- i > } > # output > sim.vector > } > > > > result is > > 1 2 3 4 5 6 7 8 11 > 0.37 0.22 0.16 0.13 0.05 0.02 0.03 0.01 0.01 > > > > always, there are some numbers missing, it should be continuous. > > why 9 and 10 are missing > > thank you > > > > sophia > > > > _________________________________________________________________ > [[elided Hotmail spam]] > > [[alternative HTML version deleted]] > > ______________________________________________ > 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.
It looks like what you did was something like:> tmp <- sim(100) > table(tmp)/100To get your results (it is easier for us to help if you tell us what you actually did, also setting a seed and telling us that seed helps us reproduce what you did exactly). The reason that you did not see 9 and 10 is just due to chance, the probabilities for those values is small and it is not surprising that you did not see a couple of unlikely values. You can either run the sim function more times or tell the computer exactly what values you want counted so it knows to include them, for example:> tmp <- factor(tmp, levels= 1:max(tmp) ) > table(tmp)/100Hope this helps, -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare greg.snow at imail.org 801.408.8111> -----Original Message----- > From: r-help-bounces at r-project.org [mailto:r-help-bounces at r- > project.org] On Behalf Of Ssophia > Sent: Wednesday, February 25, 2009 7:28 PM > To: r-help at r-project.org > Subject: [R] Hi, Coding problem > > > > Hi, there > > > > Below is my code to one Homework question. I couldn't come up with > the reasonable answer. > could you please help me to figure out what is the problem with my > code? > thank you > > > Question is Coding P{X=j} =(1/2)^(j+1) + (1/2) *2^(j-1)/3^j > my code is > sim <- function(n.gen){ > urandom <- runif(n.gen) > sim.vector <- rep(0,n.gen) > for(j in 1:n.gen){ > i <- 1 > p <- 5/12 > F <- p > while(urandom[j] >= F){ > p <- p*((1/2)^(i+1)+1/3*(2/3)^i)/((1/2)^i+(1/2)*(2/3)^i) > F <- F+p > i<-i+1 > } > sim.vector[j] <- i > } > # output > sim.vector > } > > > > result is > > 1 2 3 4 5 6 7 8 11 > 0.37 0.22 0.16 0.13 0.05 0.02 0.03 0.01 0.01 > > > > always, there are some numbers missing, it should be continuous. > > why 9 and 10 are missing > > thank you > > > > sophia > > > > _________________________________________________________________ > [[elided Hotmail spam]] > > [[alternative HTML version deleted]] > > ______________________________________________ > 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.