Анализ цепи Маркова и стационарное распределение

В этом примере показано, как вывести символьное стационарное распределение тривиальной Цепи Маркова путем вычисления ее собственного разложения.

Стационарное распределение представляет ограничение, независимое от времени, распределение состояний для процесса Маркова как количество увеличения шагов или переходов.

Задайте (положительные) вероятности перехода между состояниями A через F как показано в вышеупомянутом изображении.

syms a b c d e f cCA cCB positive;

Добавьте дальнейшие предположения, ограничивающие вероятности перехода. Это будет полезно в выборе желательных стационарных распределений позже.

assumeAlso([a, b, c, e, f, cCA, cCB] < 1 & d == 1);

Задайте матрицу перехода. Состояния A через F сопоставлены со столбцами и строками 1 через 6. Отметьте значения в каждой сумме строки до одной.

P = sym(zeros(6,6));
P(1,1:2) = [a 1-a];
P(2,1:2) = [1-b b];
P(3,1:4) = [cCA cCB c (1-cCA-cCB-c)];
P(4,4) = d;
P(5,5:6) = [e 1-e];
P(6,5:6) = [1-f f];
P
P = 

(a1-a00001-bb0000cCAcCBc1-cCA-cCB-c00000d000000e1-e00001-ff)

Вычислите все возможные аналитические стационарные распределения состояний Цепи Маркова. Это - проблема извлечения собственных векторов с соответствующими собственными значениями, которые могут быть равны 1 для некоторого значения вероятностей перехода.

[V,D] = eig(P');

Аналитические собственные вектора

V
V = 

(b-1a-10-c-dcCB-bcCA-bcCB+ccCAσ10-1010-c-dcCA-acCA-acCB+ccCBσ101000-c-dc+cCA+cCB-10000011000f-1e-1000-1010001)where  σ1=c+cCA+cCB-1a+b-ac-bc+c2-1

Аналитические собственные значения

diag(D)
ans = 

(11cda+b-1e+f-1)

Найдите собственные значения, которые точно равны 1. Если существует неоднозначность в определении этого условия для какого-либо собственного значения, остановитесь с ошибкой - этот способ, которым мы уверены, что ниже списка индексов надежно, когда этот шаг успешен.

ix = find(isAlways(diag(D) == 1,'Unknown','error'));
diag(D(ix,ix))
ans = 

(11d)

Извлеките аналитические стационарные распределения. Собственные вектора нормированы с 1 нормой или sum(abs(X)) до display.

for k = ix'
    V(:,k) = simplify(V(:,k)/norm(V(:,k)),1);
end
Probability = V(:,ix)
Probability = 

(b-1a-1σ2001σ2000000010f-1σ1e-1001σ10)where  σ1=f-12e-12+1  σ2=b-12a-12+1

Вероятность устойчивого состояния, являющегося A или B в первом собственном векторе случай является функцией вероятностей перехода a и b. Визуализируйте эту зависимость.

fsurf(Probability(1), [0 1 0 1]);
xlabel a
ylabel b
title('Probability of A');

Figure contains an axes object. The axes object with title Probability of A contains an object of type functionsurface.

figure(2);
fsurf(Probability(2), [0 1 0 1]);
xlabel a
ylabel b
title('Probability of B');

Figure contains an axes object. The axes object with title Probability of B contains an object of type functionsurface.

Стационарные распределения подтверждают следующее (Вспомните, утверждает A через F соответствуйте индексам строки 1 через 6 ):

  • C состояния никогда не достигается и поэтому переходный i.e. третья строка является полностью нулем.

  • Остальная часть состояний формирует три группы, {A B D } и {E F }, которые не связываются друг с другом и являются текущими.