机器算法验证 - 重要性采样以评估 R 中的积分 - 吾爱随笔录

重要性采样以评估 R 中的积分

机器算法验证 r 采样蒙特卡洛重要性抽样数值积分

2022-04-07 01:52:32

我也在这里问过这个问题。但是，我的理论理解可能有问题，因此我在这里问，因为它更相关。请不要在没有先看的情况下diss。

我已经提到了这些链接—— 这里和这里。后一个链接是我的问题，但事实证明，答案并不完全正确。

因此，我正在创建一个问题，并尝试回答这个问题。

我有以下积分要使用重要性采样 MC 积分解决此问题，需要选择与函数大致相同的重要性 pdf阴谋 $\int_{0.01}^{1} x^{- 0.5} d x$ $\int_{0.01}^1 x^{-0.5}\,\text{d}x$

我解决相同问题的 R 代码是这样的：

#function 1 - importance sampling
w <- function(x) dunif(x,0.01,1)/dbeta(x,0.7,1)
f <- function(x) x^(-0.5)
X <- rbeta(1000,0.7,1)
Y <- w(X)*f(X)
c(mean(Y),var(Y))

真正的整数值- 1.8

使用上面的重要性采样代码- 1.82（我的重要性 PDF 是Beta(0.7,1)

这很好，所以我假设代码是正确的。

现在我有这个积分 $\int_{0.3}^{8} [1 + \sinh (2 x) \log (x)]^{- 1} d x$ $\int_{0.3}^8 [1+\sinh(2x)\log(x)]^{-1}\text{d}x$

我的代码是：

#function 2 
w <- function(x) dunif(x,0.01,1)/dnorm(x,0.5,0.25)
f <- function(x) (1+sinh(2*x)*log(x))^(-1)
X <- rnorm(1000,0.5,0.25)
Y <- w(X)*f(X)
Y <- Y[!is.na(Y)]
c(mean(Y),var(Y))

真积分值 ~0.7014

执行上述代码的值〜3.25（我的重要性PDF是正常的（0.5，sd = 0.25）

我究竟做错了什么？

1) 取函数为 f(x)。

2) 从间隔之间截断的重要性 PDF g(x) 生成样本。

3) 得到积分的均值(f(x)/g(x))。

1个回答

运行为第二个函数提供的代码时，我得到

> c(mean(Y),var(Y))
[1] 3.2981238 0.5203621
> integrate(f,0.01,1)
3.19264 with absolute error < 1.1e-06

这意味着积分的真实值接近 3.2，而不是 0.70。

如果f要从 0.3 积分到 8，那么重要性函数必须是

> w <- function(x) dunif(x,0.3,8)/dnorm(x,0.5,0.25)

我改变了问题中 NA 的处理方式

> f <- function(x) (x>0)*(1+sinh(2*x)*log(abs(x)))^(-1)

这导致

> Y <- w(X)*f(X)
> integrate(f,0.3,8)
2.77512 with absolute error < 2.4e-05
> integrate(f,0.3,8)$val/(8-.3)
[1] 0.3604052

即使正态重要性函数 N(0.5,0.5) 可能过于集中，即它的方差太小而无法以合理的概率覆盖 (0.3,8)，这也显示出良好的拟合。如果改变法线的比例，

> X <- rnorm(1e5,0.5,2.25)
> Y <- w(X)*f(X)
> c(mean(Y),var(Y))
[1] 0.3664046 1.0230197
> integrate(f,0.3,8)$va/(8-.3)
[1] 0.3604052

这表明相当合适。当使用像这样的较小标准偏差时，我发现拟合较差 $\sigma=0.1$

> c(mean(Y),var(Y))
[1]   0.3370849 484.6246147

即使与相比： $\sigma=1$

> c(mean(Y),var(Y))
[1] 0.3606815 0.3232931

为了理解这些变化，我对一系列的值进行了实验，从到，进行了次正常模拟，并得到了下图（第一个轴上有对数刻度） $\sigma$ $0.1$ $4$ $10^6$

这说明了足够大所带来的改进 $\sigma$ 的。和一个最佳的周围 $\sigma=0.74$ .

这是相关的R代码

w <- function(x) dunif(x,0.3,8)/dnorm(x,0.5,sigma) 
sigs=varz=meanz=seq(.1,4,le=50)
for (i in 1:50){
sigma=sigs[i];X=0.5+sigma*X0;Y=w(X)*f(X)
varz[i]=var(Y);meanz[i]=mean(Y)}

其它你可能感兴趣的问题

上一篇如果测试的水平降低，测试的力量会增加吗？下一篇带有“单向”噪声的回归