[R] Counterintuitive Simulation Results

Spencer Graves spencer.graves at pdf.com
Thu Aug 4 20:20:08 CEST 2005


Hi, Robert:

	  Yes, I understand now.  Thanks very much for your insight.

	  Best Wishes,
	  Spencer Graves

McGehee, Robert wrote:

> Spencer,
> On the first iteration of your simulation, all of the Qp.t + z.t < 0, so
> you're adding a vector of rep(4.5, 20000) to a random distribution with
> mean -0.1. So one would expect on iteration 2, your mean would have
> dropped by about 0.1 (which it does). This process continues until about
> the 20th iteration when we start seeing that a large number of our
> initial starting points are floored at zero (because of the pmax). For
> points greater than zero, we continue to subtract an average of 0.1
> (actually less than this), but for those points already at zero, we're
> actually adding a mean of 0.348 (since we can never subtract from a zero
> number in this case), which starts the trajectory moving upward towards
> its asymptote.
> 
> #That is, for those paths far above 0.1, we are subtracting
> 
>>mean(rnorm(10000, mean = -0.1))
> 
> [1] -0.1059246
> 
> #And for those paths already at zero, we are adding
> 
>>mean(pmax(0, rnorm(10000, mean = -0.1)))
> 
> [1] 0.3482376
> 
> To see a simulation a bit closer to what you were expecting, replace the
> starting values with a random distribution with mean Qp0.
> 
> i.e. replace
> 
>>Qp.t <- rep(Qp0, nSims)
> 
> with
> 
>>Qp.t <- rnorm(nSims, Qp0, sd = 3.7)
> 
> 
> Robert
> 
> 
> -----Original Message-----
> From: Spencer Graves [mailto:spencer.graves at pdf.com] 
> Sent: Thursday, August 04, 2005 12:16 PM
> To: r-help at stat.math.ethz.ch
> Subject: [R] Counterintuitive Simulation Results
> 
> 
> 	  I wonder if someone can help me understand some
> counterintuitive 
> simulation results.  Below please find 12 lines of R code that 
> theoretically, to the best of my understanding, should produce 
> essentially a flat line with no discernable pattern.  Instead, I see an 
> initial dramatic drop followed by a slow rise to an asymptote.
> 
> 	  The simulation computes the mean of 20,000 simulated
> trajectories of 
> 400 observations each of a one-sided Cusum of independent normal 
> increments with mean EZ[t] = (-0.1) and unit variance.  Started with any
> 
> initial value, the mean of the Cusum should approach an asymptote as the
> 
> number of observations increases;  when started at that asymptote, it 
> should theoretically stay flat, unlike what we see here.
> 
> 	  I would think this could be an artifact of the simulation 
> methodology, but I've gotten essentially this image with several 
> independently programmed simulations in S-Plus 6.1, with all six 
> different random number generators in R 1.9.1 and 2.1.1 and with MS 
> Excel.  For modest changes in EZ[t] < 0, I get a different asymptote but
> 
> pretty much the same image.
> 
> #################################################
> simCus5 <- function(mu=-0.1, Qp0=4.5, maxTime=400, nSims=20000){
>    Qp.mean <- rep(NA, maxTime)
>    Qp.t <- rep(Qp0, nSims)
>    for(i in 1:maxTime){
>      z.t <- (mu + rnorm(nSims))
>      Qp.t <- pmax(0, Qp.t+z.t)
>      Qp.mean[i] <- mean(Qp.t)
>    }
>    Qp.mean
> }
> set.seed(1)
> plot(simCus5(Qp0=4.5))
> #################################################
> 
> 	  Thanks for your time in reading this.
> 	  Best Wishes,
> 	  Spencer Graves
> 
> Spencer Graves, PhD
> Senior Development Engineer
> PDF Solutions, Inc.
> 333 West San Carlos Street Suite 700
> San Jose, CA 95110, USA
> 
> spencer.graves at pdf.com
> www.pdf.com <http://www.pdf.com>
> Tel:  408-938-4420
> Fax: 408-280-7915
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
> 

-- 
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves at pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915




More information about the R-help mailing list