SAS中文论坛

标题: Random seeds [打印本页]

作者: shiyiming 时间: 2012-1-13 07:51
标题: Random seeds
From Dapangmao's blog on sas-analysis

A footnote toward Rick  Wilkin’s recent article “<a href="http://blogs.sas.com/content/iml/2012/01/11/how-to-lie-with-a-simulation/">How to Lie with a Simulation</a>”. (Sit in front of a laptop w/o SAS; have to port all SAS/IML codes into R)<br />
<div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-f5Ni1mnG7TA/Tw9XYtJhs6I/AAAAAAAAA5k/JoMLm3XZQys/s1600/plot1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="250" src="http://1.bp.blogspot.com/-f5Ni1mnG7TA/Tw9XYtJhs6I/AAAAAAAAA5k/JoMLm3XZQys/s400/plot1.png" width="400" /></a></div><br />
Generated 10 seeds randomly to run Stochastic simulation of Buffon's needle experiment by <a href="http://blogs.sas.com/content/iml/2012/01/04/simulation-of-buffons-needle-in-sas-2/">Rick's method</a>. Hardly converge any of them …. <br />
<pre style="background-color: #ebebeb; border: 1px dashed rgb(153, 153, 153); color: #000001; font-size: 14px; line-height: 14px; overflow: auto; padding: 5px; width: 100%;"><code>
# Replicate Rick Wicklin's SAS/IML codes for Buffon's needle experiment
simupi <- function(N, seed) {
  set.seed(seed)
  z <- matrix(runif(N*2, 0, 1), N, 2)
  theta <- pi*z[, 1]
  y <- z[, 2] / 2
  P <- sum(y < sin(theta)/2) / N
  piEst <- 2/P
  Trials <- 1:N
  Hits <- (y < sin(theta)/2)
  Pr <- cumsum(Hits)/Trials
  Est <- 2/Pr
  cbind(Est, Trials, seed)
}

# Generated 10 seeds randomly
seed_vector <- floor(runif(1:10)*10000)

# Each simulation with 50000 iterations
N <- 50000

# Run these 10 simulations
rt <- list()
for (i in 1:length(seed_vector)) {
  rt[[i]] <- simupi(N, seed_vector[i])
}
results <- as.data.frame(do.call("rbind", rt))
results$seed <- as.factor(results$seed)

# Visualize individual simulation results
library(ggplot2)
p <- qplot(x = Trials, y = Est, data = results, geom = "line",
         color = seed, ylim = c(2.9, 3.5))
p + geom_line(aes(x = Trials, y = pi), color = "Black")
ggsave("c:/plot1.png")
</code></pre><br />
<div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-QgiXx9qaaFY/Tw9XnFUYFbI/AAAAAAAAA5w/YOpTHLt6aj8/s1600/plot2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="250" src="http://2.bp.blogspot.com/-QgiXx9qaaFY/Tw9XnFUYFbI/AAAAAAAAA5w/YOpTHLt6aj8/s400/plot2.png" width="400" /></a></div><br />
<br />
Averaging all results of the 10 simulations out. Then the curve converges easily. The application of this Monte Carlo simulation in Buffon's needle experiment is explained <a href="http://blogs.sas.com/content/iml/2012/01/11/how-to-lie-with-a-simulation/#comments">here</a> by Rick Wicklin. <br />
<pre style="background-color: #ebebeb; border: 1px dashed rgb(153, 153, 153); color: #000001; font-size: 14px; line-height: 14px; overflow: auto; padding: 5px; width: 100%;"><code>
# Visuazlie the average result
rtmean <- aggregate(Est ~ Trials, data = results, mean)
p <- qplot(x = Trials, y = Est, data = rtmean, geom = "line")
p + geom_line(aes(x = Trials, y = pi), color = "Red")
ggsave("c:/plot2.png")
</code></pre><div class="blogger-post-footer"><img width='1' height='1' src='https://blogger.googleusercontent.com/tracker/3256159328630041416-6280498129719129021?l=www.sasanalysis.com' alt='' /></div><img src="http://feeds.feedburner.com/~r/SasAnalysis/~4/66FM09kx2wk" height="1" width="1"/>

欢迎光临 SAS中文论坛 (http://www.mysas.net/forum/)