这个作业是用R语言进行数据建模和蒙特卡洛算法仿真模拟
STA510 Statistical modelling and simulation
Mandatory assignment 2
问题一：
考虑积分
11
-1
11
-1
1D（x，y）dxdy（1）
其中1D（x，y）是定义的指标函数，因此
1D（x，y）=（
如果x则为1
2 + y
2≤1
否则为0。
即积分（1）给出半径为1的（二维）圆盘的面积，我们知道它具有
值π。这个问题的目的是应用蒙特卡罗积分和各种方差
减少技术来逐步找到更好的π近似值。作为粗略的第一次尝试，
考虑估计量
θCMC=
4
ñ
Xn
i = 1
1D（西
，Yi），（2）
其中Xi〜iid均匀（-1，1），Yi〜iid均匀（-1，1），Xis和Yis也独立。
a）争论为什么（2）是积分（1）的蒙特卡洛估计量。
显示1D（Xi
，Yi）具有成功概率为π/ 4的伯努利分布。
b）Pn做哪个分布
i = 1 1D（西
，易）有吗？
使用上述分布，证明（2）具有正确的期望。
找出（2）的确切方差。
c）R对n = 1000实施R中的蒙特卡洛估计（2）。获得许多独立的
π的估计。
用适当的方法检查结果是否近似正态分布（为什么？）
上一点的均值和方差。
d）R找到n = 1000的（2）正确估计π到第二个小数的概率
（即，估算值以“ 3.14”开头）。
e）根据您对蒙特卡洛中被积和对数变量的了解
积分，您是否希望引入常规的对立变量（即Vi =
a + b − Xi = −Xi
，Wi = -Yi）将减少蒙特卡洛方差。为什么或者为什么不？
f）R通过仿真验证您从e）点得出的结论是正确的。
考虑通过定义的一组不同的对立变量
shift（u）=（
-1 + u如果0≤u≤1，
如果− 1≤u <0，则为1 + u。
（3）
可以证明，当U具有均匀（-1,1）分布时，shift（U）具有均匀（-1,1）分布。使用该功能可以在R中有效地实现
shift <-function（u）{return（（（（u + 2.0）%% 2.0）-1.0）}
g）R使用对数变量（Xi
，Vi = shift（Xi））和（Yi
，Wi = shift（Yi））。
相对于原始的蒙特卡洛估计（2），这是否意味着改进？
现在我们考虑使用重要性抽样估计（1），特别是使用零均值
具有相同方差σ的独立普通建议
2
，即Xi〜N（0，σ2
），Yi〜N（0，σ2
），所以
那
f（x，y）= 1
2πσ2
经验值
–
X
2
2σ
2
</ s> </ s> </ s>
经验值
–
ÿ
2
2σ
2
</ s> </ s> </ s>
h）R使用独立的正常建议实施（1）的重要性抽样程序
f（x，y）以上。
使用σ= 0.3、0.624127、1.0测试您的实现。
相对于粗糙的蒙特卡洛估算器（2）的性能进行评论。 （注意，
在这种情况下，σ= 0.624127提供最佳的重要性抽样估算器。
Problem 2:
In this problem, imagine that you work for an insurance company. Further, you are responsible
for setting aside enough capital to, with high probability, be able to fulfill claims originating from
storm damage.
You are modeling the number of storms as non-homogenous Poisson process, with intensity
λ(t) = 297 (1 + cos (2π (t + 1/10))) (1 − exp(−t/10)/2)
10
+ 3/5 (4)
where time t is measured in years since Jan 1st 2020 (i.e. t = 0 correspond to Jan 1st 2020,
t = 1.5 is July 1st 2021 and so on)1
. As can be seen from the specification of λ(t), there is a
strong seasonal effect, and also the winter intensity is increasing due to climate change.
In the following point, you may find the R-function integrate useful.
a)R Let X be the number of storms in the year 2020 (i.e. t ∈ [0, 1)). Which distribution does
X have?
What is the expected number of storms in 2025 (i.e. t ∈ [5, 6))?
What is the expectation and standard deviation of the total number of storms in 2020 and
2021 combined.
b) Find the smallest possible λmax so that λmax ≥ λ(t) for all t ≥ 0.
c)R Write a function that simulates the times storms occur in some time interval [a, b] where
a ≥ 0.0 and b > a.
Check your algorithm by replicating the expectations and standard deviation in point a)
using simulations.
Now assume that the size of total claims2
in million Norwegian kroner, from any given storm has
an exponential distribution with mean c(t), where
c(t) = 10 exp
5t
100
(5)
Here, t is the time the storm occurred (also in years after Jan 1st 2020). It is seen that the size
of the claims increases, which reflect both economic growth and that insurance company expects
to grow in the future.
d)R Based the function written in point c) above, write a function that simulates the total claim
size resulting from all storms in some given time interval [a, b].
Find the mean and standard deviation of the claim size resulting from storms in the year
2020.
Find a confidence 95% interval for the mean.
How much funds must the insurance company set aside to be 97.5% certain to be able to
cover the total claim size of the year 2020?
Suppose you wish to improve the estimate of the mean of the total claim size. A common technique
for reducing variance in Monte Carlo integration, Rao-Blackwelliztion, may be though of as a way
to exploit the law of total expectation3
to reduce the variance. In the present case, this would
amount to exploiting that we know the expectation of sums of exponentially distributed random
variables. Suppose X is the total claim size originating in time span [a, b], and N is the number
of storms occurring in the time span [a, b]. Then X may be written as
X =
(PN
i=1 c(ti)εi
if N > 0
0 if N = 0
where ti
, i = 1, . . . , N are the times the storms occur and the εis iid exponentially distributed
with mean 1. Then the expectation of X may be written as
E(X) = E(E(X|N, t1, . . . , tN )) = E (Y ) (6)
where
E(X|N, t1, . . . , tN ) = Y =
(PN
i=1 c(ti) if N > 0
0 if N = 0
(7)
e)R Use (6,7) to implement another Monte Carlo estimator for E(X).
Compare this Monte Carlo estimator to the one obtained in point d) for the year 2020.
f)R More difficult: Based on the law of total variance V ar(X) = E(V ar(X|N, t1, . . . , tN )) +
V ar(E(X|N, t1, . . . , tn)), propose and implement an improved estimator for the variance of
X. Compare this estimator to the one in point d) (where you simply simulate X) for the
year 2020