用Python计算马尔可夫链的平稳分布

Question

我已经在谷歌foobar问题上工作了几天,并且只有一次测试通过了,而且我在这一点上非常困难.如果您有任何想法,请告诉我!我使用描述的方法在这里,我有一个工作例如多达上repl.it 这里.这是问题规范:

世界末日燃料

为LAMBCHOP的反应堆堆芯制造燃料是一个棘手的过程,因为涉及到奇特的物质.它从原始矿石开始,然后在加工过程中,开始在形式之间随机变化,最终达到稳定的形式.样品最终可能存在多种稳定形式,并非所有形式都可用作燃料.

Lambda指挥官已经责成您通过预测给定矿石样品的最终状态来帮助科学家提高燃料产生效率.您仔细研究了矿石可以采取的不同结构以及它经历的转变.似乎随机,每个结构转换的概率是固定的.也就是说,每次矿石处于1状态时,它具有进入下一状态(可能是相同状态)的相同概率.您已在矩阵中记录观察到的过渡.实验室中的其他人假设矿石可以变成更奇特的形式,但你还没有看到它们全部.

编写一个函数answer(m),它接受一个非负整数数组的数组,表示该状态进入下一个状态的次数,并为每个终端状态返回一个int数组,给出每个终端状态的确切概率,表示为每个州的分子,然后是最后和最简单形式的所有分母.矩阵最多为10乘10.确保无论矿石处于哪种状态,都存在从该状态到终端状态的路径.也就是说,处理总是最终以稳定状态结束.矿石在状态0开始.只要分数被定期简化,分母将在计算期间符合带符号的32位整数.

*例如,考虑矩阵m:

[
  [0,1,0,0,0,1],  # s0, the initial state, goes to s1 and s5 with equal probability
  [4,0,0,3,2,0],  # s1 can become s0, s3, or s4, but with different probabilities
  [0,0,0,0,0,0],  # s2 is terminal, and unreachable (never observed in practice)
  [0,0,0,0,0,0],  # s3 is terminal
  [0,0,0,0,0,0],  # s4 is terminal
  [0,0,0,0,0,0],  # s5 is terminal
]

因此,我们可以考虑到终端状态的不同路径,例如:

s0 -> s1 -> s3
s0 -> s1 -> s0 -> s1 -> s0 -> s1 -> s4
s0 -> s1 -> s0 -> s5

追踪每个的概率,我们发现s2有概率0 s3有概率3/14 s4有概率1/7 s5有概率9/14所以,把它放在一起,并做一个共同的分母,给出一个答案的形式[s2.numerator,s3.numerator,s4.numerator,s5.numerator,denominator]是 [0, 3, 2, 9, 14].*

测试用例

Inputs:
    (int) m = [[0, 2, 1, 0, 0], [0, 0, 0, 3, 4], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]]
Output:
    (int list) [7, 6, 8, 21]

Inputs:
    (int) m = [[0, 1, 0, 0, 0, 1], [4, 0, 0, 3, 2, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]]
Output:
    (int list) [0, 3, 2, 9, 14]

到目前为止,这是我的代码.

from __future__ import division
from itertools import compress
from itertools import starmap
from operator import mul
import fractions


def convertMatrix(transMatrix):
  probMatrix = []
  for i in range(len(transMatrix)):
    row = transMatrix[i]
    newRow = []
    rowSum = sum(transMatrix[i])
    if all([v == 0 for v in transMatrix[i]]):
      for j in transMatrix[i]:
        newRow.append(0)
      newRow[i] = 1
      probMatrix.append(newRow)
    else:
      for j in transMatrix[i]:
        if j == 0:
          newRow.append(0)
        else:
          newRow.append(j/rowSum)
      probMatrix.append(newRow)
  return probMatrix



def answer(m):
  # convert matrix numbers into probabilities
  probMatrix = convertMatrix(m)

  # find terminal states
  terminalStateFilter = []
  for row in range(len(m)):
    if all(x == 0 for x in m[row]):
      terminalStateFilter.append(True)
    else:
      terminalStateFilter.append(False)

  # multiply matrix by probability vector
  oldFirstRow = probMatrix[0]
  probVector = None
  for i in range(3000):
    probVector = [sum(starmap(mul, zip(oldFirstRow, col))) for col in zip(*probMatrix)]
    oldFirstRow = probVector

  # generate numerators
  numerators = []
  for i in probVector:
    numerator = fractions.Fraction(i).limit_denominator().numerator
    numerators.append(numerator)

  # generate denominators
  denominators = []
  for i in probVector:
    denominator = fractions.Fraction(i).limit_denominator().denominator
    denominators.append(denominator)

  # calculate factors to multiply numerators by
  factors = [max(denominators)/x for x in denominators]
  # multiply numerators by factors
  numeratorsTimesFactors = [a*b for a,b in zip(numerators, factors)]
  # filter numerators by terminal state booleans
  terminalStateNumerators = list(compress(numeratorsTimesFactors, terminalStateFilter))

  # append numerators and denominator to answer
  answerlist = []
  for i in terminalStateNumerators:
    answerlist.append(i)
  answerlist.append(max(denominators))

  return list(map(int, answerlist))