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

Поглощающее(существенное) состояние цепи Маркова - состояние с вероятностью перехода в самого себя [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]_[math]num[/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_num++;

Найдем число несущественных состояний [math]nonabs=n-abs[/math]_[math]num[/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
   N[i][i]=1;'
   E[i][i]=1;
   for j=0 to nonabs
     E[i][j]-=Q[i][j];  

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

 for i = 0 to nonabs
    if E[i][i] != 1
       mul = E[i][i];
       for j = 0 to nonabs
         E[i][j] /= mul;
         N[i][j] /= mul;
     for row = 0 to nonabs
         if i != row
            mul = E[row][i];
            for j = 0 to nonabs
                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
     for j = 0 to absorbing
         G[i][j] = 0;
         for k = 0 to nonabs
             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] где j - номер соответствующий состоянию в [math]G[/math]. Прибавлять 1 нужно т.к. вероятность поглотиться в [math]i[/math]-ом поглощающем состоянии, оказавшись изначально в нем же равна 1.

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

Литература

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