问题标签 [mcmc]

我想知道是否有人知道如何创建 mcmc 采样并计算由包制作的线性混合模型的 HPDlme间隔nlme。我熟悉使用lme4包执行此操作，但我需要来自nlme.

我有一个关于在 JAGS 和 BUGS 中运行模型的详细信息的快速问题。

n.burnin=5000假设我使用和运行n.iter=5000模型thin=2。这是否意味着该程序将：

1. 运行 5,000 次迭代，并丢弃结果；接着
2. 再运行 10,000 次迭代，只保留每秒的结果？

如果我将这些模拟保存为一个CODA对象，是全部保存了 10,000 个，还是只保存了稀疏的 5,000 个？我只是想了解使用哪组迭代来制作 ACF 图？

也许这里对 mcmc 对象有更多经验的人可以帮助我。

问题：我有一个存储了 20 多个 mcmc 对象的列表。我需要为每个 mcmc 对象分配变量名。

我有另一个列表，其中包含存储的每个 mcmc 的所有变量名称（作为数据框中的列）。

我可以使用 coda 包中的“varnames”函数单独执行此操作，如下所示：

作为“投票”我的数据框列表和“后验”我的 mcmc 对象列表..

但是，我不想一个接一个地做这件事，我想一次做所有的事情。我试过下面的代码......

但我明白了Error in *tmp*[[x]] : Recursive indexing failed at level 2。我已经尝试了该行的一些变体，但我很难理解如何正确索引它，或者如何做我想做的事情。

我知道这是一个特定的问题，但也许这里有人可以给我一个提示或什么。

在此先感谢您的帮助。抱歉，我无法提供一些数据，但很难获得可行的样本。

问候，费德里科

I am currently working on my final year project for my mathematics degree which is based on giving an overview of the Metropolis-Hastings algorithm and some numerical examples. So far I have got some great results by using my proposal distribution as a Gaussian, and sampling from a few other distributions, however I am trying to go one step further by using a different proposal distribution.

So far I have got this code (I am using Matlab), however with limited resources online about using different proposals it is hard to tell if I am close at all, as in reality I am not too sure how to attempt this, (especially as this gives no useful data output so far).

It would be fantastic if someone could lend a hand if they know or forward me to some easily accessible information (I understand that I am not just asking coding advice, but Mathematics as well).

So, I want to sample from a Gaussian using a proposal distribution of a Laplace, this my code so far:

If anyone can tell me if the above is completely wrong/going about it in the wrong way, or there are just a few mistakes (I am very wary about my generation of 'y' in for the for loop being completely wrong) that would be fantastic.

Thanks, Tom

代码如下：

我不知道为什么，当我替换为时程序会正确theta[j] <- b*x[j-1]，theta[j] <- 1*x[j-1]但我已经定义了b ~ dunif(-1, 1). 确实，我需要theta[j] <- a - b*x[j-1]在最终模型中进行设置，当我尝试添加a并b进入它时结果是错误的。有谁发现问题出在哪里？

我使用 PyMC 的时间不长，但我很高兴能够以多快的速度实现线性回归（这段代码应该在 IPython 中无需修改即可运行）：

在这个模型中，有 40 个主题（观察）和每个主题的 5 个协变量。由于随机数据，该模型不会收敛，但它的采样没有错误（我的真实数据确实会收敛到准确的结果）。

我遇到问题的模型是这个的扩展。每个受试者实际上有 3 个（或 N 个）观察结果，因此我需要为这些观察结果拟合一条线，然后使用该线的截距作为“数据”或最终回归节点的输入。我认为这是一个经典的层次模型，但如果我以错误的方式思考它，请纠正我。

我的解决方案是设置一系列 40 个线性回归（每个主题一个），然后使用截距参数向量作为最终回归的数据。

该模型在采样步骤失败。错误似乎是在尝试将 offsetArr 转换为 dtype=float64 而不是它的 dtype=object 时。这样做的正确方法是什么？我是否需要另一个确定性节点来帮助将我的 offsetArr 转换为 float64？我需要一种特殊的 pymc.Container 吗？谢谢你的帮助！

我的winbugs代码如下：

运行此代码后，我得到了 alpha 和 beta 的后验分布。现在我想看看P(beta>0)。他们说我可以使用pbeta<- step(beta)（pbeta 被视为虚拟变量： 0 ifbeta=0和 1 if beta>0）。但是当我把它放在模型中时，它给了我一个错误通知。

如果 pymc 实现 Metropolis-Hastings 算法以从感兴趣的参数上的后验密度中提取样本，那么为了决定是否移动到马尔可夫链中的下一个状态，它必须能够评估与后验成比例的东西所有给定参数值的密度。

后验密度与基于观察数据乘以先验密度的似然函数成比例。

这些在 pymc 中是如何表示的？它如何从模型对象中计算出每一个数量？

我想知道是否有人可以给我该方法的高级描述或指出我在哪里可以找到它。

我正在研究 NVIDIA CUDA GPU 上的马尔可夫链蒙特卡洛 (MCMC) 算法实现。CPU MCMC 算法使用高质量的 Mersenne twister 随机数生成器，我想在我编写的 GPU 内核中使用相同的。我一直在寻找 cuRand MT 代码示例。不幸的是，我从未见过任何使用 Mersenne twister 的内核代码示例。标准的 cuRand 库文档为 MTGP（MT for Graphic Processor）提供了一组函数，但不清楚如何使用它们。

CUDA 示例为 MersenneTwisterGP11213.tar.gz 提供了一个示例，但它似乎专门用于请求在 GPU 上快速生成随机数数组、将它们下载到 CPU 内存并在 CPU 上进行的主机代码。还有一篇论文“Massively Parallel RNG using CUDA C, Thrust and C#”。同样，上一节“使用 CUDA C 的 Mersenne Twister 实现”中的作者仅提供了来自“CUDA 示例”的上述主机代码的简化片段。

所以，我的第一个问题是：谁能给我一个使用 cuRand Mersenne twister 的全局或设备函数示例？

我还有一个问题。目前我使用 cuRand 库随机数生成器，但我不知道使用的是什么生成器！让我提供几段我的代码。这是生成器初始化：

在其他内核中，我从均匀分布和正态分布中采样数字。所有blockDim.x*gridDim.x生成器的状态数组都保存在全局内存数组mc->rndst[]中。例如，curand_uniform()用于：

或者，从高斯分布中采样，curand_normal()使用：

谁能告诉我这里使用了哪些 cuRand 生成器（xorwow、lcs、mtgp ...）（实际上，默认情况下）？

我正在使用psclR 中的包并试图让它产生可测试/可重现的结果。我查看了底层 C 代码，它看起来好像在正确GetRNGstate()的PutRNGstate()位置被调用，但似乎不可能重复 MCMC 模型的输出。

simulationResult我已经从SoDA包中打包了函数，这样我就可以在 R 端验证每个仿真 R 的启动状态。

我们可以验证至少在 R 端的起始状态是相同的：

但是输出不同：

我可以增加迭代次数，但如果 R​​NG 状态正在传播，这并不重要。我错过了一些非常简单的东西吗？谢谢。

编辑：我还应该注意all.equal(run1@lastState, run2@lastState)（每次运行的最终状态）应该是相同的，但它们最终会不同。我的猜测是，C 调用的 R RNG 函数之外的一些意外事件来源正在影响调用这些 RNG 函数的次数。好奇的。

编辑2

我还应该在 OS X 10.8.4 上添加我在 R 3.0.1 和 pscl 1.04.4 上。