一个简单的马尔可夫链
假设我们想要估计一个系统的参数,这样我们就可以在给定时间步 t 的状态的情况下预测时间步 t+1 的系统状态。 PyMC 应该能够轻松处理这个问题。
让我们的玩具系统由一维世界中的移动物体组成。状态是对象的位置。我们要估计潜在变量,即物体的速度。下一个状态取决于前一个状态和潜在变量速度。
# define the system and the data
true_vel = .2
true_pos = 0
true_positions = [.2 * step for step in range(100)]
我们假设我们的观察中有一些噪音(但这在这里无关紧要)。
问题是:我如何为下一个状态对当前状态的依赖关系建模。我可以为转换函数提供一个参数 idx 以访问时间 t 的位置,然后预测时间 t+1 的位置。
vel = pymc.Normal("pos", 0, 1/(.5**2))
idx = pymc.DiscreteUniform("idx", 0, 100, value=range(100), observed=True)
@pm.deterministic
def transition(positions=true_positions, vel=vel, idx=idx):
return positions[idx] + vel
# observation with gaussian noise
obs = pymc.Normal("obs", mu=transition, tau=1/(.5**2))
但是,索引似乎是一个不适合索引的数组。可能有更好的方法来访问以前的状态。
最佳答案
最简单的方法是生成一个列表,让PyMC把它当作一个Container来处理。有一个相关example on the PyMC wiki .这是相关的片段:
# Lognormal distribution of P's
Pmean0 = 0.
P_0 = Lognormal('P_0', mu=Pmean0, tau=isigma2, trace=False, value=P_inits[0])
P = [P_0]
# Recursive step
for i in range(1,nyears):
Pmean = Lambda("Pmean", lambda P=P[i-1], k=k, r=r: log(max(P+r*P*(1-P)-k*catch[i-1],0.01)))
Pi = Lognormal('P_%i'%i, mu=Pmean, tau=isigma2, value=P_inits[i], trace=False)
P.append(Pi)
请注意当前对数正态分布的均值是上一个对数分布的函数吗?不优雅,使用 list.append
和所有,但你可以改用列表推导。
关于python - PyMC:马尔可夫系统中的参数估计,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/20222971/