使用 Monte Carlo 方法从密度核计算均值和方差

Question

假设概率分布的密度核为，我可以使用哪些 Monte Carlo 方法来估计分布的均值和方差？

Answer 1

我们可以在这里使用数值方法。首先，我们创建一个函数来表示您的概率密度函数（尽管尚未对其进行缩放，因此其整个域的积分为 1）：

pdf <- function(x) x^2 * exp(-x^2/4)

plot(pdf, xlim = c(0, 10))

为了获得将其转换为真正的 pdf 的比例因子，我们可以在其 c(0, Inf).

域上集成此函数

integrate(pdf, 0, Inf)$value
#> [1] 3.544908

所以现在我们可以通过将我们的原始 pdf 除以这个数量来生成一个真正的 pdf:

pdf <- function(x)  x^2 * exp(-x^2/4) / 3.544908

plot(pdf, xlim = c(0, 10))

现在我们有了pdf，我们可以创建一个带有数值积分的cdf：

cdf <- function(x) sapply(x, \(i) integrate(pdf, 0, i)$value)

plot(cdf, xlim = c(0, 10))

cdf 的逆函数是我们需要的，它能够将从 0 和 1 之间的均匀分布中提取的样本转换为从我们的新分布中提取的样本。我们可以使用 uniroot 创建这个反函数来查找 cdf 的输出与 0 到 1 之间的任意数字相匹配的位置：

inverse_cdf <- function(p) 
{
  sapply(p, function(i) {
              uniroot(function(a) {cdf(a) - i}, c(0, 100))$root
  })
}

逆 cdf 看起来像这样：

plot(inverse_cdf, xlim = c(0, 0.99))

我们现在准备好从我们的分布中抽取样本：

set.seed(1) # Makes this draw reproducible

x_sample <- inverse_cdf(runif(1000))

现在我们可以绘制样本的直方图并确保它与 pdf 匹配：

hist(x_sample, freq = FALSE)

plot(function(x) pdf(x), add = TRUE, xlim = c(0, 6))

现在我们有了从 x 中抽取的样本，我们可以使用样本均值和方差作为分布均值和方差的估计值：

mean(x_sample)
#> [1] 2.264438

var(x_sample)
#> [1] 0.9265678

我们可以通过在调用 inverse_cdf(runif(1000)) 时采用越来越大的样本来提高这些估计的准确性，方法是将 1000 增加到更大的数字。

^{由 reprex package (v2.0.0)}

于 2021-11-06 创建

Calculating mean and variance from a density kernel using Monte Carlo methods