simpi                  package:pscl                  R Documentation

_M_o_n_t_e _C_a_r_l_o _e_s_t_i_m_a_t_e _o_f _p_i (_3._1_4_1_5_9_2_6_5...)

_D_e_s_c_r_i_p_t_i_o_n:

     Monte Carlo estimation of pi

_U_s_a_g_e:

     simpi(n)

_A_r_g_u_m_e_n_t_s:

       n: integer, number of Monte Carlo samples, defaults to 1000

_D_e_t_a_i_l_s:

     A crude Monte Carlo estimate of pi can be formed as follows. 
     Sample from the unit square many times (i.e., each sample is
     formed with two independent draws from a uniform density on the
     unit interval).  Compute the proportion p of sampled points that
     lie inside a unit circle centered on the origin; such points (x,y)
     have distance from the origin d=sqrt(x^2 + y^2) less than 1.  Four
     times p is a Monte Carlo estimate of pi.  This function is a
     wrapper to a simple C function, bringing noticable speed gains and
     memory efficiencies over implementations in native R.

     Contrast this Monte Carlo method with Buffon's needle and
     refinements thereof (see the discussion in Ripley (1987, 193ff).

_V_a_l_u_e:

     the Monte Carlo estimate of pi

_A_u_t_h_o_r(_s):

     Simon Jackman jackman@stanford.edu

_R_e_f_e_r_e_n_c_e_s:

     Ripley, Brain D. 1987 [2006].  _Stochastic Simulation_. Wiley:
     Hoboken, New Jersey.

_E_x_a_m_p_l_e_s:

     seed <- round(pi*10000)  ## hah hah hah
     m <- 6
     z <- rep(NA,m)
     lim <- matrix(NA)
     for(i in 1:m){
       cat(paste("simulation for ",i,"\n"))
       set.seed(seed)
       timings <- system.time(z[i] <- simpi(10^i))
       print(timings)
       cat("\n")
       lim[i] <- qbinom(prob=pi/4,size=10^i,.975)/10^i * 4
     }

     ## convert to squared error
     z <-(z - pi)^2
     lim <- (lim - pi)^2

     plot(x=1:m,
          y=z,
          type="b",
          pch=16,
          log="y",
          axes=FALSE,
          ylim=range(z,lim),
          xlab="Monte Carlo Samples",
          ylab="Log Squared Error")
     lines(1:m,lim,col="blue",type="b",pch=1)
     legend(x="topright",
            legend=c("95% bound",
              "Realized"),
            pch=c(1,16),
            lty=c(1,1),
            col=c("blue","black"),
            bty="n")
     axis(1,at=1:m,
          labels=c(expression(10^{1}),
            expression(10^{2}),
            expression(10^{3}),
            expression(10^{4}),
            expression(10^{5}),
            expression(10^{6})))
     axis(2)

