关于编写计算马尔可夫联合分布脚本的两个问题（python）

Question

我是python的新手，最近正在做一些项目来计算马尔可夫过程的联合分布。

An example of a stochastic kernel is the one used in a recent study by Hamilton (2005), who investigates a nonlinear statistical model of the business cycle based on US unemployment data. As part of his calculation he estimates the kernel

pH :=   0.971 0.029 0
        0.145 0.778 0.077
        0     0.508 0.492

Here S = {x1, x2, x3} = {NG, MR, SR}, where NG corresponds to normal growth, MR to mild recession, and SR to severe recession. For example, the probability of transitioning from severe recession to mild recession in one period is 0.508. The length of the period is one month.

基于上述马尔可夫过程的习题是

With regards to Hamilton’s kernel pH, and using the same initial condition ψ = (0.2, 0.2, 0.6) , compute the probability that the economy starts and remains in recession through periods 0, 1, 2 (i.e., that xt = NG fort = 0, 1, 2).

我的脚本是这样的

import numpy as np
## In this case, X should be a matrix rather than vector
## and we compute w.r.t P rather than merely its element [i][j]
path = []
def path_prob2 (p, psi , x2):  # X a sequence giving the path
     prob = psi                # initial distribution is an row vector
     for t in range(x2.shape[1] -1): # .shape[1] grasp # of columns
        prob = np.dot(prob , p)       # prob[t]: marginal distribution at period t
        ression = np.dot(prob, x2[:,t])
     path.append(ression)
     return path,prob



p = ((0.971, 0.029, 0    ),
      (0.145, 0.778, 0.077),
      (0    , 0.508, 0.492)) 
# p must to be a 2-D numpy array     
p = np.array(p)      

psi = (0.2, 0.2, 0.6)  
psi = np.array(psi)  
x2 = ((0,0,0),
      (1,1,1),
      (1,1,1))
x2 = np.array(x2)      
path_prob2(p,psi,x2)

在执行过程中，出现了两个问题。第一个是，在第一轮循环中，我不需要初始分布psi来postmultiply交易矩阵p，所以"remaining in recession"的概率应该是0.2+0.6 = 0.8，但我不知道如何编写 if 语句。 第二个是，正如您可能注意到的，我使用一个名为 path 的列表来收集每个时期 "remaining in recession" 的概率。最后我需要将列表中的每个元素一个一个地相乘，我没有找到一种方法来实现这样的任务，比如 path[0]*path[1]*path[2] （np.multiply 到目前为止只能接受两个参数我所知）。如果确实存在这种方法，请给我一些线索。另一个问题是请给我任何您认为可以使代码更高效的建议。谢谢。

Answer 1

如果我对你的理解正确，这应该可行（我喜欢对某些 steps/outcome 进行一些手动计算），请注意我没有使用 if/else 语句而是从第二列开始迭代：

import numpy as np

# In this case, X should be a matrix rather than vector
# and we compute w.r.t P rather than merely its element [i][j]
path = []


def path_prob2(p, psi, x2):  # X a sequence giving the path
    path.append(np.dot(psi, x2[:, 0]))  # first step
    prob = psi  # initial distribution is an row vector
    for t in range(1, x2.shape[1]):  # .shape[1] grasp # of columns
        prob = np.dot(prob, p)  # prob[t]: marginal distribution at period t
        path.append(np.prod(np.dot(prob, t)))
    return path, prob


# p must to be a 2-D numpy array
p = np.array([[0.971, 0.029, 0],
              [0.145, 0.778, 0.077],
              [0, 0.508, 0.492]])

psi = np.array([0.2, 0.2, 0.6])
x2 = np.array([[0, 0, 0],
               [1, 1, 1],
               [1, 1, 1]])

print path_prob2(p, psi, x2)

对于你的第二个问题，我相信 Numpy.prod 会给你一个 list/array.

的所有元素之间的乘法

您可以这样使用产品：

>>> np.prod([15,20,31])
9300

关于编写计算马尔可夫联合分布脚本的两个问题（python）

Two problems on writing a script to compute markov joint distribution (in python)

python

numpy

markov-chains

mcmc

markov