R语言由Monte Carlo方法计算积分

98 阅读 0 评论 65 点赞

我是靠谱客的博主乐观心锁，这篇文章主要介绍R语言由Monte Carlo方法计算积分，现在分享给大家，希望可以做个参考。

蒙特·卡罗方法（Monte Carlo method），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。

给大家分享三道用Monte Carlo方法解决的概率统计问题：

设每日向保险公司申请理赔的人数服从均值为 10 的Poisson 分布。每次理赔的数量的大小服从均值为 1000 的指数分布。另一方面保险公司每天收到 11000 的保险付费。假设保险公司的最初的资金为 25000。使用 Monte Carlo 方法计算在第一个 365 天，公司的资产恒为正的概率。

复制代码

1
2
3
4
5
6
7
8
9
10
11
f = function(o) sum(rexp(o, rate = 1/1000))
g = Vectorize(f)
p = function(n){
b= sapply(1:n, function(o) {
N = rpois(365, lambda = 10)
F = cumsum(11000 - g(N)) + 25000
d = I(sum(F>0)==365)
d})
c(mean = mean(b), sd = sd(b))
}

设一个系统由 20 个独立的部件组成，部件 $i$ 失效的概率为 $0.5 + i 50 , i = 1 , 2 , . . . , 20 0.5 + frac{i}{50}, i = 1,2,...,20$ 以 $X$ 表示系统失效部件的总数。请使用简单 Monte Carlo 方法和 importance sampling 的方法估计 $P (X < = 5)$ ，并估计由两种方法得到的概率的估计的标准差。

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 简单 monte carlo 方法
p1 = function(n){P = 0.5 + (1:20)/50
b = sapply(1:n, function(o){
S = rbinom(20, 1, prob = P)
X = sum(S)
d = I(X<=5)
d})
c(mean = mean(b), sd = sd(b))}
# Importance sampling 方法
p2 = function(n){P = 0.5 + (1:20)/50
smp = sapply(1:n, function(o){
S = rbinom(20, 1, prob = 0.2)
a = (P^S)*((1-P)^(1-S))
b = (0.2^S)*(0.8^(1-S))
X = prod(a/b)
d = X*I(sum(S)<=5)})
c(mean = mean(smp), sd = sd(smp))}

设随机变量X的密度函数为 $e^{-sqrt{x}}sin^2x, x > 0$
我们希望使用 importance sampling 的方法估计上述随机变量的期望。令函数g(x)为 $frac{1}{2Pi(1+frac{x^2}{4})}$ , $x \in R$
请完成下列任务：
1)使用长为 M = 100, 1000, 100000 的 $g (x)$ 的随机数序列，估计密度 $f (x)$ 的期望。对上面的每次估计，重复100次，计算估计的标准差 $σ$ 并比较三个标准差的大小；
2)求使得估计的标准差 $σ$ <0.01的M的大小。

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
f <- function(x) exp(- x^(1/2))*(sin(x))^2
g <- function(x)
(1/(2*pi))*(1/(1 + x^2/4))*2
ratio <- function(x) f(x)/g(x)
# 首先看看如何产生服从柯西分布的正半轴的密度函数的随机数
h <- function(n) abs(rcauchy(n, scale = 2))
a = h(100000)
hist(a[a<=40], freq = FALSE, breaks = seq(0, 40, 0.05), border = "darkblue")
curve(g, from = 0, to = 40, add = TRUE, col = "red")
# 下面用 importance sampling 方法生成服从 f(x) 的随机数
Smp <- function(m,n){d = sapply(1:m, function(o)
{ s <- abs(rcauchy(n, scale = 2))
a <- s*ratio(s)
mean(a)})
c(mean = mean(d), sd = sd(d))}
Smp(100, 100)
Smp(100, 1000)
Smp(100, 100000)
# 检查上面结果是否正确
w = function(x) x*f(x)
integrate(w, lower = 0, upper = 1000)
# 注意到均值的标准差是和根号 n 成反比。
a = NULL
for(i in seq(1000, 10000, 1000)) a = c(a, Smp(100, i)["sd"])
# 得 n = 7000