I am a new R user. As a test, I want to write a simple code that does the following simulation: 1. Randomly generate a number from a distribution, say, Poisson. Let's say that number is 3. 2. Randomly generate 3 numbers from another distribution, say, Normal. 3. Compute the sum of the numbers generated in step 2 and read it into a vector, V. 4. Repeat steps 1 through 3 for 100,000 times. 5. Derive quantiles (e.g., 0.95th, 0.99th) of V. Any help in getting me going would be greatly appreciated. [[alternative HTML version deleted]]
Try this ... poisNums <- rpois(100000, 3) # Calculate 100,000 poissons (mean 3) normNums <- sapply(poisNums, rnorm, mean=3) # Calculate N random normals # for each simulated poisson. Normal distribution is mean 3, sd 1 sumNums <- sapply(normNums, sum) mySeq <- seq(0.85, 1, length=50) plot(seq(0.85, 1, length=50)*100, quantile(sumNums, mySeq), type="l", main="Quantiles of simulated data", xlab="Quantiles", ylab="") abline(h=quantile(sumNums, c(0.95, 0.99)), col=2, lwd=2) R does this for me in about 20 seconds ... Rich. Mango Solutions Tel : (01628) 418134 Mob : (07967) 808091 -----Original Message----- From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of jonathan_wang at sbcglobal.net Sent: 22 February 2004 17:18 To: r-help at stat.math.ethz.ch Subject: [R] Simulation help I am a new R user. As a test, I want to write a simple code that does the following simulation: 1. Randomly generate a number from a distribution, say, Poisson. Let's say that number is 3. 2. Randomly generate 3 numbers from another distribution, say, Normal. 3. Compute the sum of the numbers generated in step 2 and read it into a vector, V. 4. Repeat steps 1 through 3 for 100,000 times. 5. Derive quantiles (e.g., 0.95th, 0.99th) of V. Any help in getting me going would be greatly appreciated. [[alternative HTML version deleted]] ______________________________________________ R-help at stat.math.ethz.ch mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html --- Incoming mail is certified Virus Free.
How about the following:
> set.seed(5)
> N <- 8 # later 100000
> (nPois <- rpois(N, 2))
[1] 1 3 4 1 0 3 2 3
> z <- rnorm(sum(nPois))
> V <- tapply(z, rep(1:N, nPois), sum)
> quantile(V, c(0, .05, .25, .5, .75, .95, 1))
0% 5% 25% 50% 75% 95%
100%
-2.7812799 -2.4600296 -1.5393631 -1.0803926 0.5800796 1.4626007
1.7114409
The same code works in S-Plus 6.2 and R 1.8.1 under Windows 2000,
though the answers different as S-Plus and R use different random number
generators.
hope this helps. spencer graves
jonathan_wang at sbcglobal.net wrote:
>I am a new R user. As a test, I want to write a simple code that does the
following simulation:
>
>1. Randomly generate a number from a distribution, say, Poisson. Let's
say that number is 3.
>2. Randomly generate 3 numbers from another distribution, say, Normal.
>3. Compute the sum of the numbers generated in step 2 and read it into a
vector, V.
>4. Repeat steps 1 through 3 for 100,000 times.
>5. Derive quantiles (e.g., 0.95th, 0.99th) of V.
>
>Any help in getting me going would be greatly appreciated.
>
> [[alternative HTML version deleted]]
>
>______________________________________________
>R-help at stat.math.ethz.ch mailing list
>https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
>
>
In a reply to my inquiry, a respondent offered the following simulation code:> lens <- rpois(100000, 3)> V <- numeric(100000)> for(i in 1:length(lens)) V[i] <- sum(rnorm(lens[i]))> quantile(V, c(.95, .99))The code worked exactly the way I wanted. I had an Excel model that does the same thing. The quantiles from both models were almost identical. Two follow-up questions: 1. I also did a mean(V), but the result I got in R was 5 times greater than the result I got using the Excel model. Any idea why that might be? 2. When I tried using the random generation from the Generalized Pareto Distribution, I kept getting an error as follows: "Error in rgpd(lens[i]) : Argument "xi" is missing, with no default" My code was:> for(i in 1:length(lens)) V[i] <- sum(rgpd(lens[i]),0.5,mu=0.1,beta=0.8)Note: rgpd() is part of library(evir); a Windows version of which can be downloaded at http://www.maths.lancs.ac.uk/~stephena/software/win.html Many thanks! [[alternative HTML version deleted]]