时间戳

首页 > 解决方案 > 使用蒙特卡罗方法从密度核计算均值和方差

r - 使用蒙特卡罗方法从密度核计算均值和方差

问题描述

假设概率分布的密度核为https://i.stack.imgur.com/irjRs.png，我可以使用哪些蒙特卡罗方法来估计分布的均值和方差？

标签： rstatisticsmeankdevariance

解决方案

我们可以在这里使用数值方法。首先，我们创建一个函数来表示您的概率密度函数（尽管它尚未缩放，因此它的积分在其整个域上为 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 包于 2021-11-06 创建 (v2.0.0 )

推荐阅读