在给定建议函数的情况下使用重要性采样进行蒙特卡洛积分

2023-12-11

给定拉普拉斯分布提案:

g(x) = 1/2*e^(-|x|)

和样本量n = 1000,我想进行蒙特卡罗(MC)积分来估计 θ:

enter image description here

通过重要性采样。最终,我想在到达那里后计算 R 中 MC 估计的平均值和标准差。


编辑(在回答以下问题后迟到)

这是我到目前为止的 R 代码:

library(VGAM)
n = 1000
x = rexp(n,0.5)
hx = mean(2*exp(-sqrt(x))*(sin(x))^2)
gx = rlaplace(n, location = 0, scale = 1)

enter image description here

现在我们可以编写一个简单的 R 函数来从拉普拉斯分布中进行采样:

## `n` is sample size
rlaplace <- function (n) {
  u <- runif(n, 0, 1)
  ifelse(u < 0.5, log(2 * u), -log(2* (1 - u)))
  }

还编写一个拉普拉斯分布密度函数:

g <- function (x) ifelse(x < 0, 0.5 * exp(x), 0.5 * exp(-x))

现在,你的被积函数是:

f <- function (x) {
  ifelse(x > 0, exp(-sqrt(x) - 0.5 * x) * sin(x) ^ 2, 0)
  }

现在我们使用 1000 个样本来估计积分(set.seed为了重现性):

set.seed(0)
x <- rlaplace(1000)
mean(f(x) / g(x))
# [1] 0.2648853

还可以与使用求积的数值积分进行比较:

integrate(f, lower = 0, upper = Inf)
# 0.2617744 with absolute error < 1.6e-05
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

在给定建议函数的情况下使用重要性采样进行蒙特卡洛积分 的相关文章

  • R闪亮的html小部件之间的交互

    我正在开发一个 R 闪亮应用程序 它使用多个 html 小部件 特别是网络D3 d3热图 and 和弦诊断 这些小部件单独工作正常 但是 在同一页面中使用它们会留下一个空格处他们应该在哪里 这是显示错误的可重现代码 在 UI 中注释绘图线
  • Shiny 中的模态对话框:可以调整宽度但不能调整高度

    在我的 Shiny 应用程序中 我有几个来自闪亮BS 包的模式窗口 我可以像这样调整这些模式窗口的宽度 tags head tags style HTML modal lg width 1200px abs 1 background col
  • 从 R 主题模型中的 DocumentTermMatrix 中删除空文档?

    我正在使用 R 中的 topicmodels 包进行主题建模 我正在创建一个 Corpus 对象 进行一些基本的预处理 然后创建一个 DocumentTermMatrix corpus lt Corpus VectorSource vec
  • 在 for 循环中绘制的多个 ggplot2 绘图的网格

    作为一个新的 ggplot2 用户 我对可能性的数量感到有点迷失 并且很难在网上找到我认为简单问题的简单答案 我想在同一张纸上显示 ggplot2 的多个图 但知道这些图来自 for 循环 以下示例无法编译 仅用于说明 for i in c
  • 如何从 data.frame 中选择行和列的子集

    我有这个 d d Age gt 2 它返回 Age 超过 2 的所有行 但我只想返回几列中的值 例如 d X 和 d Y 而不是全部 无论如何我可以做到这一点吗 Thanks d d Age gt 2 c X Y
  • R 中大型稀疏矩阵的聚类分析

    我有一个包含 250000 笔交易 行 和 2183 项 列 的交易数据集 我想将其转换为稀疏矩阵 然后对其进行分层聚类 我尝试了包 sparcl 但它似乎不适用于稀疏矩阵 关于如何解决这个问题有什么建议吗 或者我可以使用任何其他包对稀疏矩
  • 有效地将环境从内部功能转移到全局环境

    我有一个在其中创建环境的函数 我希望将该环境分配给全局环境 目前我通过将环境分配给来做到这一点globalenv 作为最后一步 如下 funfun lt function inc 1 dataEnv lt new env dataEnv d
  • glm() 模型的交叉验证

    我正在尝试对我之前在 R 中构建的一些 glm 模型进行 10 倍交叉验证 我对cv glm 函数在boot包 尽管我已经阅读了很多帮助文件 当我提供以下公式时 library boot cv glm data glmfit K 10 这里
  • R 中 write.table 文件名中的变量

    请帮助我解决一个幼稚的问题 已经用谷歌搜索 并尝试了很多变体 但失败了 如何使用 R 中 write table 的文件名中的变量保存文件 脚本循环遍历 dir 中的文件 应用一些函数 然后将结果保存到具有相同名称但附加结尾的文件中 谢谢
  • mclapply 用户时间大于已用时间

    我正在尝试使用mclapply的功能parallel封装在R 该函数通过计算对数似然距离将值分配给序列矩阵 这是一个 CPU 密集型操作 所结果的system time价值观令人困惑 gt system time mclapply work
  • 按绝对值排序

    有谁知道如何按绝对值对 R 中的向量进行排序 所以 2 3 1 gt 1 2 3 etc 如果我在 python 中这样做 我会创建一对每个值及其符号 按绝对值对对列表进行排序 然后重新应用符号 但我对 R 很陌生 所以不知道如何执行此操作
  • 在捕食者-被捕食者系统的生态建模中正确使用 deSolve

    我有一个捕食者 被捕食者模型 其中包含指定的参数和初始值 我在这里用两种方法求解微分方程 1 使用 for 循环 2 使用 deSolve 包 我相信 for 循环是正确的 并且应该给出如下图所示的输出 For loop attempt r
  • R中有字典功能吗

    有没有办法在 R 中创建一个 字典 使其具有对 一些效果 x dictionary c Hi Why water c 1 5 4 x Why 5 我问这个是因为我实际上正在寻找两个分类变量函数 所以如果 x dictionary c a b
  • 我可以调整scale_color_brewer的下限吗?

    我已经订购了我想使用 color Brewer 的分类数据 但我很难看到非常低的值 有没有办法去掉这些较低的值或设置范围的下限 ggplot data frame x 1 6 y 10 15 w letters 1 6 aes x y co
  • 使用亚毫秒日期时间从字符->POSIXct->字符准确转换

    我的文件中有一个字符日期时间列 我加载文件 到data table 并执行需要将列转换为的操作POSIXct 然后我需要写POSIXct值返回文件 但日期时间不会相同 因为打印不正确 这个打印 格式问题是众所周知的 并且已经被讨论过多次 我
  • 如何在 R 中查找平衡面板数据(又名,如何查找面板中的哪些条目在给定窗口内完整)

    我有来自 Compustat 的大量数据 我向其中添加了一些手工收集的数据 认真地从一堆旧书中手工收集 但我不想手工收集整个面板 只想随机选择一个子集 为了找到更大的集合 我从中随机选择 我想从 Compustat 的平衡面板开始 我看到p
  • 创建序列组合

    我正在尝试解决以下问题 考虑 5 个简单序列 0 100 100 0 rep 0 101 rep 50 101 rep 100 101 我需要 3 个数字变量的集合 它们的所有组合都具有上述序列 由于有 5 个序列和 3 个变量 因此可以有
  • 在 R Shiny 中显示/隐藏整个框元素

    我目前正在尝试找到一种方法来隐藏 显示 R Shiny 中的整个 box 元素 以及里面的所有内容 我想创建一个可能的按钮 它允许用户展开特定框 然后使用相同 甚至不同 的按钮隐藏它 我不想使用条件面板 因为我的应用程序非常大并且会产生一些
  • 在 ggplot 中过滤管道 df

    我正在使用 dplyr 管道来清理我的 df 然后直接输入到 ggplot 中 但是 我只想一次只绘制一组 因此我需要过滤到该组 问题是 我希望比例保持不变 就好像所有群体都存在一样 是否可以在 ggplot 命令中进一步过滤管道 df 例
  • 使用插入符和方法 = gamLoess 进行训练时 R 崩溃

    当我运行下面的代码时 R 崩溃了 如果我在训练调用中注释掉tuneGrid行 就不会发生崩溃 我已经用另一个数据集尝试过此操作 但仍然使 R 崩溃 崩溃消息是 R 会话中止 R遇到致命错误 会话被终止 开始新会话 代码是 library s

随机推荐