标签: mcmc

我有一个模型,如下图所示:

我有几个人口(在这张照片中索引1 ... 5).人口参数(A和B,但可以有更多)确定每个人的潜在变量的分布L[i].潜在变量以概率方式L[i]确定观察X[i].在大多数节点没有直接连接它们的边缘的意义上,该模型是"稀疏的".

我试图使用PyMC来推断人口参数,以及每个人的潜在变量.(一个相关的问题,其中详细描述了我的具体情况,是在这里.)我的问题是:我应该使用自适应大都市,而不是另一种方法,如果是这样,有没有"猫腻",以正确分组随机变量？

如果我理解正确的自适应都市报采样(和我可能不会...),该算法提出了一种未知的新值(A,B以及所有L[i]通过考虑这些变量是如何在目前为止在运行构造的后验分布相关) .如果A且B是负相关的,那么增加的提案A将倾向于减少B,反之亦然,以增加提案被接受的机会.

问题是,在这个模型中,每个L[i]都是由A和确定的基础人口分布的独立抽取B.因此,虽然他们将被视为在后方相关联,但这些相关性实际上是由于A并且是B单独的,因此它们以某种方式"混淆".所以当我调用这个函数时

M.use_step_method(pymc.AdaptiveMetropolis, stochastics)

是否所有人都应该L[i]在随机指标列表中？或者我应该多次调用use_step_method,每次stochastics=[A, B, L[i]]只调用一个L[i]？我的想法是,对于不同的随机指标组多次调用该函数将构成问题并使PyMC更容易通过告诉它只关注重要的相关性.它是否正确？

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

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

如何在pymc中代表这些？它如何从模型对象计算这些数量？

我想知道是否有人能给我一个关于这种方法的高级描述,或者指出我能找到它的地方.

我是matlab的死硬用户,主要是因为这是我第一次学到的东西而且我没有遇到过切换差异很大的问题.我来自数值优化/线性代数,我已经在数百万个自由度中进行了优化和特征值计算.最近,我进入了随机性的境界,我最初的印象是我不得不改变.但是,在优化代码并小心地将种子初始化为随机数生成器之后,我能够在大致相同的时间内完成与同时代相同的蒙特卡罗任务.我的理解是基础级别'if'语句等在matlab中明显变慢.但是,如果在每个循环中有可以进行矢量化的重要计算,我' 我不相信C会更好.而且,无论如何,matlab似乎做得很好,因为我的上限不亚于任何其他专业(在很多情况下,似乎更多).我有一种感觉,我会从这里的亲C人那里得到很多回复,他们很久以前就已经把matlab写成了一些琐碎的玩具语言.我是一名专业研究人员,并认为matlab对于最高级别的计算密集型数学编程具有竞争力.我错了吗 ？我是否需要考虑更改为较低级别的语言,例如C/FORTRAN？为什么或者为什么不 ？有没有像我这样的人？非常感谢!干杯 我有一种感觉,我会从这里的亲C人那里得到很多回复,他们很久以前就已经把matlab写成了一些琐碎的玩具语言.我是一名专业研究人员,并认为matlab对于最高级别的计算密集型数学编程具有竞争力.我错了吗 ？我是否需要考虑更改为较低级别的语言,例如C/FORTRAN？为什么或者为什么不 ？有没有像我这样的人？非常感谢!干杯 我有一种感觉,我会从这里的亲C人那里得到很多回复,他们很久以前就已经把matlab写成了一些琐碎的玩具语言.我是一名专业研究人员,并认为matlab对于最高级别的计算密集型数学编程具有竞争力.我错了吗 ？我是否需要考虑更改为较低级别的语言,例如C/FORTRAN？为什么或者为什么不 ？有没有像我这样的人？非常感谢!干杯

首先,我要说我缺乏科学数学或统计学的经验 - 所以这可能是一个众所周知的问题,但我不知道从哪里开始.

我有一个函数f(x1, x2, ..., xn),我需要猜测x'并找到最高值f.该函数具有以下属性:

• 总数或参数通常在40到60左右,因此蛮力方法是不可能的.

• 每个x的可能值范围为0.01到2.99

• 该函数是稳定的,这意味着较高的f值意味着参数的猜测更好,反之亦然.

到目前为止,我在python中实现了一个非常基本的方法.它最初将所有参数设置为1,随机猜测新值并检查f是否高于之前.如果没有,请回滚到先前的值.在具有10,000次迭代的循环中,这似乎以某种方式起作用,但结果可能远非完美.

关于如何改进搜索最佳参数的任何建议将不胜感激.当谷歌搜索这个问题时,linke MCMC出现了,但这似乎是一种非常先进的方法,我需要花费大量时间来理解这种方法.基本提示或概念对我的帮助不仅仅是详细阐述方法和算法.

是否可以在python中传递带有一些参数的模块？

所有我的意思是参数是模块中存在一个未在该模块中初始化的变量,我仍然在该模块中使用该变量.简而言之,我希望行为类似于函数,但与函数不同,我希望模块的变量在调用代码中公开.

例如

a.py

#lists like data, count, prob_distribution are constructed from training_pool (not initialized in this file)
x = pymc.Uniform('x', lower = 0, upper = 1)
rv = [ Multinomial("rv"+str(i), count[i], prob_distribution[i], value = data[i], observed=True) for i in xrange(0, len(count)) ]

b.py

import a  #I want some way tr pass value of training_pool
m = pymc.MCMC(a)

我希望a.py中的所有随机变量都暴露给MCMC.我对我手头的问题采取了更好的方法,但我也想知道在python中是否可以传递模块的参数

我正在寻求关于如何使用Metropolis-Hastings算法将C或C++代码合并到我的R代码中来加速MCMC程序的建议.我使用MCMC方法来模拟可能性,给定各种协变量,个人将被第三方(法官)分配给社会地位等级中的特定等级:每个法官(大约80个,跨越4个村庄)被要求根据对每个人的社会地位的评估,对一组个人(大约80个,跨越4个村庄)进行排名.因此,对于每个法官,我有一个等级向量,对应于他们对每个人在等级中的位置的判断.

为了模拟这种我认为,分配的行列时,法官对个人效用,一些潜在措施的相对价值根据他们的决定ü.鉴于此,然后可以假定行列,向量[R ,由给定的法官产生是未观测到的矢量的函数,ü,描述个人的效用被排名,其中与个人第k的最高值ü将被分配第k级.我建模ü,,,使用感兴趣的协变量作为多变量正态分布变量,然后确定所观察到的行列的可能性给定的分布ü由模型生成的.

除了估计最多5个协变量的影响之外,我还估计了描述法官和项目之间差异的超参数.因此,对于链的每次迭代,我估计多变量法向密度大约为8-10倍.因此,5000次迭代可能需要长达14个小时.显然,我需要运行它超过5000次运行,所以我需要一种方法来大大加快这个过程.鉴于此,我的问题如下:

(i)我是否有权假设通过在C或C++中运行一些(如果不是全部)链条来获得最佳速度增益？

(ii)假设问题1的答案是肯定的,我该怎么做呢？例如,有没有办法让我保留所有R函数,只需在C或C++中循环:即我可以从C调用R函数然后循环吗？

(iii)我想我真正想知道的是如何最好地将C或C++代码合并到我的程序中.

在pymc3中如何配置截断的正常先验？在pymc2中,它非常简单(下图),但在pymc3中似乎不再有截断的正态分布.

Pymc2:

TruncatedNormal('gamma_own_%i_' % i, mu=go, tau=v_gamma_inv, value=0, a=-np.inf, b=0)

Pymc3 :?

我正在使用psclR中的软件包并尝试使其生成可测试/可重现的结果.我已经看到了在底层的C代码,它看起来好像GetRNGstate()和PutRNGstate()被称为在正确的地方,但它似乎是不可能从MCMC模型重复输出.

我已经simulationResult从SoDA包中打包了函数,因此我可以验证R端的每个模拟R的启动状态.

library(pscl)
library(SoDA)
run1 <- simulationResult(
  ideal(s109, 
    normalize=TRUE,
    maxiter = 500,
    thin = 10,
    burnin = 0),
  seed = 42)

run2 <- simulationResult(
  ideal(s109, 
    normalize=TRUE,
    maxiter = 500,
    thin = 10,
    burnin = 0),
  seed = 42)

我们可以验证起始状态至少在R方面是相同的:

all.equal(run1@firstState, run2@firstState)

但输出是不同的:

all.equal(run1@result$xbar, run2@result$xbar)

我可以增加迭代次数,但如果RNG状态得到传播则这并不重要.我错过了一些非常简单的事吗？谢谢.

编辑:我还应该注意all.equal(run1@lastState, run2@lastState)(每次运行的结束状态)应该是相同的但它们最终会有所不同.我的猜测是,被C称为R RNG功能外应急的一些源被撞击数的倍那些RNG函数被调用.好奇.

EDIT2

我还应该在OS X 10.8.4上使用pscl 1.04.4添加我的R 3.0.1.

我是PyMC的新手,并尝试使用最大后验估计来拟合我的非齐次泊松过程和分段恒定速率函数.

我的过程描述了一天中的一些事件.因此,我将一天分成24小时,这意味着,我的速率函数中有24个常数(分段常数).

结合以下思路:

• 黑客的概率规划和贝叶斯方法:第1章
• SO问题:使用可变速率拟合带顶点的泊松过程

我提出了以下一段代码,这是不满意的(结果明智,我确定这是错的):

import numpy as np
import pymc

eventCounter = np.zeros(24) # will be filled with real counts before going on
alpha = 1.0 / eventCounter.mean()

a0 = pymc.Exponential('a0', alpha)
a1 = pymc.Exponential('a1', alpha)
a2 = pymc.Exponential('a2', alpha)
a3 = pymc.Exponential('a3', alpha)
a4 = pymc.Exponential('a4', alpha)
a5 = pymc.Exponential('a5', alpha)
a6 = pymc.Exponential('a6', alpha)
a7 = pymc.Exponential('a7', alpha)
a8 = pymc.Exponential('a8', alpha)
a9 = pymc.Exponential('a9', alpha)
a10 = pymc.Exponential('a10', alpha) …

我正在尝试使用MCMCglmmR中的包创建模型

数据结构如下,其中二元,焦点,其他都是随机效应,预测1-2是预测变量,响应1-5是捕获不同子类型的观察行为#的结果变量:

 dyad focal other r    present  village  resp1 resp2 resp3 resp4 resp5 
 1    10101 14302 0.5  3        1        0     0     4     0     5
 2    10405 11301 0.0  5        0        0     0     1     0     1
 …

所以只有一个结果(教学)的模型如下:

 prior_overdisp_i <- list(R=list(V=diag(2),nu=0.08,fix=2), 
 G=list(G1=list(V=1,nu=0.08), G2=list(V=1,nu=0.08), G3=list(V=1,nu=0.08), G4=list(V=1,nu=0.08)))

 m1 <- MCMCglmm(teaching ~ trait-1 + at.level(trait,1):r + at.level(trait,1):present, 
 random= ~idh(at.level(trait,1)):focal + idh(at.level(trait,1)):other + 
 idh(at.level(trait,1)):X + idh(at.level(trait,1)):village, 
 rcov=~idh(trait):units, family = "zipoisson", prior=prior_overdisp_i, 
 data = data, nitt = nitt.1, thin = 50, burnin = 15000, …