Расчёт вероятности поглощения в состоянии

Поглощающее(существенное) состояние цепи Маркова - состояние с вероятностью перехода в самого себя [math]p_{ii}=1[/math]. Составим матрицу G, элементы которой [math]g_{ij}[/math] равны вероятности того, что, выйдя из i, попадём в поглощающее состояние j.

Теорема:

[math] G = N \cdot R [/math]

Доказательство:

[math]\triangleright[/math]

Пусть этот переход будет осуществлён за r шагов: i → [math]i_{1}[/math] → [math]i_{2}[/math] → ... → [math]i_{r-1}[/math] → j, где все [math]i, i_{1}, ... i_{r-1}[/math] являются несущественными. Тогда рассмотрим сумму [math]\sum\limits_{\forall(i_{1} ... i_{r-1})} {p_{i, i_{1}} \cdot p_{i_{1}, i_{2}} \cdot ... \cdot p_{i_{r-1}, j}} = Q^{r-1} \cdot R[/math], где Q - матрица переходов между несущественными состояниями, R - из несущественного в существенное. Матрица G определяется их суммированием по всем длинам пути из i в j: [math]G = \sum\limits_{r = 1}^{\infty}{Q^{r-1} \cdot R} = (I + Q + Q^{2} + Q^{3} + ...) \cdot R = NR[/math], т.к. [math](I + Q + Q^2 + ...) \cdot (I - Q) = I - Q + Q - Q^{2} + ... = I[/math], а фундаментальная матрица марковской цепи [math]N = (I - Q)^{-1}[/math]

[math]\triangleleft[/math]

Псевдокод

Пусть [math]n[/math] - количество состояний Марковской цепи, [math]m[/math] - количество переходов. Состояния пронумерованы от [math]0[/math] до [math]n - 1[/math], переходы от [math]0[/math] до [math]m - 1[/math]. Входные данные хранятся в массиве [math]input[/math] где [math]i[/math]-ая строка характеризует [math]i[/math]-ый переход таким образом: [math]input[i][2][/math] - вероятность перехода из состояния [math]input[i][0][/math] в состояние [math]input[i][1][/math]. Создадим массив [math]absorbing[/math] типа boolean, где [math]i[/math]-ое true обозначает что [math]i[/math]-ое состояние является поглощающим и наоборот. Обнаружим поглощающие состояния по такому признаку: если состояние поглощающее то с вероятностью 1 оно переходит само в себя. Также посчитаем количество поглощающих состояний [math]abs[/math].

for i=0 to n - 1
   if input[i][0] == input[i][1] and input[i][2] == 1
      absorbing[input[i][0]] = true;
      abs++;

Найдем число несущественных состояний [math]nonabs = n - abs[/math]. Теперь нужно заполнить матрицы [math]Q[/math] (переходов между несущественными состояниями) и [math]R[/math] (переходов из несущественных состояний в поглощающие). Для этого создадим сначала массив [math]position[/math] где [math]i[/math]-ый элемент указывает под каким номером будет находиться [math]i[/math]-ое состояние среди существенных если оно существенное или несущественных в обратном случае, и заполним эти массивы.

count_q = 0;
count_r = 0;
for i = 0 to n - 1
   if abs[i]
      position[i] = count_r;
      count_r++;
   else 
      position[i] = count_q;
      count_q++;
for i = 0 to m - 1
   if absorbing[input[i][1]]
      if absabsorbing[input[i][0]]
         R[position[input[i][0]]][position[input[i][1]]] = input[i][2];
   else
      Q[position[input[i][0]]][position[input[i][1]]] = input[i][2];

Найдем Матрицу [math]E = I - Q[/math] и создадим единичную матрицу [math]N[/math].

for i = 0 to nonabs - 1
   N[i][i] = 1;
   E[i][i] = 1;
   for j = 0 to nonabs - 1
      E[i][j] -= Q[i][j];  

Теперь приведем матрицу [math]E[/math] к единичной методом Гаусса - Жордана, применяя те же преобразования к матрице [math]N[/math].

for i = 0 to nonabs - 1
   if E[i][i] != 1
      mul = E[i][i];
      for j = 0 to nonabs - 1
         E[i][j] /= mul;
         N[i][j] /= mul;
   for row = 0 to nonabs - 1
      if i != row
         mul = E[row][i];
         for j = 0 to nonabs - 1
            E[row][j] -= mul * E[i][j];
            N[row][j] -= mul * N[i][j];

В результате [math]N = E^{-1}[/math] т.е. [math]N[/math] - фундаментальная матрица Марковской цепи. Найдем матрицу [math]G = N * R[/math].

for i = 0 to nonabs - 1
   for j = 0 to absorbing - 1
      G[i][j] = 0;
      for k = 0 to nonabs - 1
         G[i][j] += N[i][k] * R[k][j];

Выведем ответ: в [math]i[/math]-ой строке вероятность поглощения в [math]i[/math]-ом состоянии. Естественно для несущественного состояния это 0, в ином случае [math]p=(($$\sum_{k=1}^n G[k][j]$$)+1)/n[/math] где [math]j[/math] - номер соответствующий состоянию в [math]G[/math]. Прибавлять 1 нужно т.к. вероятность поглотиться в [math]i[/math]-ом поглощающем состоянии, оказавшись изначально в нем же равна 1.

for i = 0 to n - 1
   prob = 0;
   if absorbing[i]
      for j = 0 to nonabs - 1
         prob += G[j][position[i]];
      prob++;
      prob /= n;
   println(prob);

Литература

• Википедия - Цепи Маркова
• Кемени Дж., Снелл Дж. "Конечные цепи Маркова".