标题

p103 8(2)
p103 10
p144 5
p144 7
p134 13

p103 8(2)

用模拟的方法近似计算下列积分，并和已知的精确答案进行比较：

∫ 0 ∞ x ( 1 + x 3 ) − 4 d x \int^{\infty }_{0}x(1+x^3)^{-4}dx ∫0∞x(1+x3)−4dx

解答：

f1 = function(n,a,b,g){x = runif(n)sum((b-a)*g(a+(b-a)*x))/n
}
f = function(x) x*(1+x^3)^-4
f1(99999,0,Inf,f)

得到结果

[1]NaN

。因此我们需要近似计算积分。将Inf改为100试试：

[1] 0.2055264

我们使用integrate函数获得精确答案：

integrate(f,0,Inf)

0.2089975 with absolute error < 1.7e-05

误差不大。

p103 10

对(0,1)上的均匀随机变量U1、U2…，定义

N = m i n { n : ∑ i = 1 n U i > 1 } N=min \{n:\sum^{n}_{i=1} Ui>1\} N=min{n:∑i=1nUi>1}.

即，N等于使其和超过1的随机数的个数。

解答：

这题比较简单，可以直接写：

f=function(x){EN = 0sumN = 0for(i in 1:x){N = 0i = 0while(N <= 1){rand = runif(1)N = N + randi = i + 1}sumN = sumN + i}EN = sumN/x;EN
}

我猜理论值是2.72。

p144 5

在用随机投点法估计 p i pi pi时，做多少次模拟才能获得一个长度小于0.1的区间，使之以95%的可靠性包括 p i pi pi？

解答：
我们先写一个函数使用随机投点法估计pi。

f=function(x){k=9;u1=runif(x);u2=runif(x)X=2*u1-1;Y=2*u2-1for(i in 1:x){if(X[i]^2+Y[i]^2<=1) k=k+1}4*k/x
}

我们很容易可以写出测试程序：

f=function(x){k=9;u1=runif(x);u2=runif(x)X=2*u1-1;Y=2*u2-1;res=0;num=1for(i in 1:x){if(X[i]^2+Y[i]^2<=1) k=k+1}res=4*k/xwhile(abs(res-pi)>=0.1||abs((1-res/pi))>0.05){num = num + 1k=9;u1=runif(x);u2=runif(x)X=2*u1-1;Y=2*u2-1;res=0for(i in 1:x){if(X[i]^2+Y[i]^2<=1) k=k+1}res=4*k/x}res
}

很明显，这一次所有的值都在题目范围内，接下来我们只需要输出模拟次数就可以了（把res改成num）：

逐步更改x值，我们可以让它逐渐迅速而准确地得到pi值。（以上为改x值时对应的模拟次数）

p144 7

为了估计 θ \theta θ，产生了20个相互独立的具有均值 θ \theta θ的随机数，其值如下：

解答：

我们观察这组数字，得到其平均值为122.85，平均值与最值差为29.85。

我们不妨使用代码：

round(rnorm(20,122.85,29.85))

可以大致模拟。如果要确定置信系数为0.99长度为1的最终估计量，我们这组数据仍然是不够的。我们很容易写出：

f=function(x,adv,sd){sumN=sum(round(rnorm(x,adv,sd)))advN=sumN/xtryN=1while(abs(advN-adv)>1||abs(advN/adv)<0.99){tryN=tryN+1sumN=sum(round(rnorm(x,adv,sd)))advN=sumN/x}tryN
}

我们不妨假设 θ \theta θ为120，标准差为30。

可以看出一个产生1个，由于达不到要求，会额外生成20～100不定的模拟次数才能最终达到要求。

如果一次产生10个，大约需要产生70～80个随机数。

p134 13

一个意外伤亡保险公司有1000个客户，每个客户独立地在下个月以概率0.05索赔.假设索赔量是独立的具有均值$ 800指数随机变量.用模拟方法估计这些索赔量的和超过$ 50000 的概率。

解答：

作1000次模拟，设X为下个月需要索赔的客户数，易知X服从参数为(1000,0.05)的二项分布；设Y表示每个客户需要索赔的数额，由题意知，Y服从参数为1/800的指数随机变量。所求的概率即：

k=0
for(i in 1:1000){X=rbinom(1,1000,0.05)Y=rexp(X,1/800)if(sum(Y)>50000){k=k+1}
}
k/1000

[1] 0.116

R语言：作业八（用函数进行近似计算）相关推荐

R语言使用download.file函数下载网络文件到本地（Download File from the Internet）
R语言使用download.file函数下载网络文件到本地(Download File from the Internet) 目录 R语言使用download.file函数下载网络文件到本地(Down ...
R语言stringr包str_dup函数字符串多次复制实战
R语言stringr包str_dup函数字符串多次复制实战目录 R语言stringr包str_dup函数字符串多次复制实战 #导入stringr包 #仿真数据
R语言stringr包str_count函数计算字符串匹配个数实战
R语言stringr包str_count函数计算字符串匹配个数实战目录 R语言stringr包str_count函数计算字符串匹配个数实战 #导入stringr包 #仿真数据
R语言ggpubr包ggsummarystats函数可视化分组条形图(自定义分组颜色、添加抖动数据点jitter、误差条)并在X轴标签下方添加分组对应的统计值（样本数N、中位数、四分位数的间距iqr)
R语言ggpubr包ggsummarystats函数可视化分组条形图(自定义分组颜色.添加抖动数据点jitter.误差条error bar)并在X轴标签下方添加分组对应的统计值(样本数N.中位数med ...
R语言使用t.test函数计算两组独立数据的t检验（Independent t-test）
R语言使用t.test函数计算两组独立数据的t检验(Independent t-test) 目录 R语言使用t.test函数计算两组独立数据的t检验(Independent t-test) #仿真数据
R语言plyr包round_any函数将向量数据近似到任意精度实战
R语言plyr包round_any函数向量将数据近似到任意精度实战目录 R语言plyr包round_any函数向量将数据近似到任意精度实战 #导入plyr包 #仿真数据
R语言ggplot2使用geom_line函数geom_point函数可视化哑铃图、并对哑铃图进行排序（reorder dumbbell plot）
R语言ggplot2使用geom_line函数geom_point函数可视化哑铃图.并对哑铃图进行排序(reorder dumbbell plot) 目录
R语言stringr包str_detect函数检测字符串中模式存在与否实战
R语言stringr包str_detect函数检测字符串中模式存在与否实战目录 R语言stringr包str_detect函数检测字符串中模式存在与否实战 #导入stringr包
R语言使用upper.tri函数、lower.tri函数、diag函数改变matrix矩阵上三角形、下三角形、对角线的数值
R语言使用upper.tri函数.lower.tri函数.diag函数改变matrix矩阵上三角形.下三角形.对角线的数值目录
R语言使用read.delim函数读取带分隔符的文本文件
R语言使用read.delim函数读取带分隔符的文本文件目录 R语言使用read.delim函数读取带分隔符的文本文件 #read.delim函数语法

R语言：作业八（用函数进行近似计算）

标题