1 蒙特卡洛法（Monte Carlo method）计算\(\pi\)的值

用随机函数runif()随机生成N个在1个单位正方形内的散点，根据三角公式筛选出在圆内的点,数据如下:

## 生成基础数据

set.seed(2021)

N <-  5000

df01 <- tibble(
  x = runif(N, -0.5, 0.5),
  y = runif(N, -0.5, 0.5)
) %>% 
  mutate(
    z = if_else(
      x^2 + y ^2 < 0.25, 1, 0),
    p = runif(N,0,2 * pi),
    x2 = cos(p) * 0.5,
    y2 = sin(p) * 0.5
    )

knitr::kable(head(df01), caption = "随机数据", col.names = c("散点X", "散点Y", "是否在圆内",
                                                          "随机生成0-180度", "圆X", "圆Y"), digits = 3, align = "c" )

1.1 绘图

画出散点的分布图。

## 绘图
df01 %>% ggplot() +
  geom_point(aes(x, y, col= as.factor(z)),
             alpha = 0.5, size = 0.5, show.legend = F) + 
  theme_bw() +
  geom_point(aes(x2,y2), alpha = 0.5, size = 0.5)

图1.1: 1个单位正方形内的散点

1.2 计算\(\hat\pi\)

根据面积比可得\(\hat\pi=4\frac{^{ }S_c}{S_s}\)。当分布5000个点的时候,\(\hat\pi=\frac{3925\times4}{5000}=\) 3.14。以实验投掷散点数量为X轴，以估的\(\hat\pi\)为Y轴绘图如下。很明显地可以看出随着点的数量增多\(\hat\pi\)逐渐向\(\pi\)真实值（蓝线）逼近。

estmatepi <- tibble(y = cumsum(df01$z)/(1:N)*4,
                    x = 1:N)
estmatepi %>% ggplot(aes(x,y))+
  geom_line()+
  geom_hline(yintercept = pi, col = "blue") +
  theme_classic() +
  xlab("N") +
  ylab(expression(hat(pi)))

图1.2: 估计圆周率与实验次数关系