Набор инструментов не требуется:
Мы знаем, что вероятности стационарного состояния являются уникальным решением приведенных ниже уравнений для цепи Маркова с дискретным временем (DTM C).
Вы можете переписать это, чтобы показать, что это распределение должно l ie в пустом пространстве матрицы вероятностей перехода (см. Ниже). Это позволяет использовать MATLAB null().
Нормализация суммы дает единственное решение и удовлетворяет другому уравнению.
P = [.2 .3 .3 0 .2;
.1 .2 .3 .4 0;
0 .5 .4 .1 0;
.2 .2 .2 .3 .1;
.5 .2 .1 .1 .1];
ssp0 = null(eye(size(P))-P.');
SteadyStateProb = ssp0./sum(ssp0) % [0.1263 0.3017 0.2969 0.2224 0.0528]
Обратите внимание, что промежуточный ssp0 имеет отрицательные значения, но окончательный ответ, SteadyStateProb, является верным в том смысле, что все значения являются действительными вероятностями, а все вероятности составляют 1.
Обновление на основе правки ОП : подготовка матрицы
Одна из проблем с матрицей вероятности перехода, размер 84 × 84 TPM в редактировании , что сумма строки больше единицы. Вы можете проверить это с помощью find(sum(TPM,2) > 1), чтобы получить индексы строк, где это происходит. Кроме того, этот «избыток» очень мал: max(sum(TPM,2))-1 показывает, что наибольшая сумма строк превышает 1 на 4.4409e-16, для меня это в два раза больше eps (2.2204e-16).
Используя приведенный ниже очень технологичный код, я нормализовал каждую строку.
% Re-Normalize: Divide the kth row by the sum of that row
TPM2 = TPM;
for k = 1:size(TPM2,1)
TPM2(k,:) = TPM2(k,:)/sum(TPM2(k,:));
end
В результате сумма строк не будет больше единицы; вы можете проверить с помощью find(sum(TPM2,2) > 1), который возвращает пустой двойной вектор-столбец 0 × 1.
Без дальнейшего массирования матрицы мы можем получить распределение в устойчивом состоянии.
ssp0 = null(eye(size(TPM2))-TPM2.');
SteadyStateProb_TPM2 = ssp0./sum(ssp0);
Обратите внимание, что min(SteadyStateProb_TPM2) равно -9.6502e-16 (ноль!), А сумма по существу равна 1 (abs(1-sum(SteadyStateProb_TPM2))).
Я бы сейчас просто навязал неотрицательность.
SteadyStateProb_TPM2(SteadyStateProb_TPM2<0) = 0;
Сравните с итерационным приближением:
Мы знаем, что TPM2^n сходится к SteadyStateProb_TPM2, когда n становится достаточно большим. Итак, давайте проверим это.
Мы умножаем матрицу перехода многократно (здесь мы получаем матрицу перехода с 2500 шагами). Мы можем видеть, что она кажется более или менее сходящейся к цели.
Approx = TPM2^500;
Approx = Approx*Approx*Approx*Approx*Approx;
ind = find(SteadyStateProb_TPM2 > 0.05); % find states with steady state prob > 0.05
SteadyStateProb_TPM2(ind).'
Approx(1:6,ind)
Обратите внимание, что состояния [4 8 12 16 20 24].
ans =
0.0821 0.2052 0.2565 0.2138 0.1336 0.0668
ans =
0.0727 0.1818 0.2273 0.1894 0.1184 0.0592
0.0760 0.1900 0.2374 0.1979 0.1237 0.0618
0.0763 0.1908 0.2385 0.1988 0.1242 0.0621
0.0821 0.2052 0.2565 0.2138 0.1336 0.0668
0.0727 0.1818 0.2273 0.1894 0.1184 0.0592
0.0760 0.1900 0.2374 0.1979 0.1237 0.0618
В зависимости от вашего приложения вы можете считать этот метод итерации достаточно хорошим, но я бы не стал использовать его, когда методы прямого анализа c все еще жизнеспособны. Я полагаю, что матрица 2600 × 2600 выполнима в MATLAB, если вам нужно время, чтобы убедиться, что она хорошо подготовлена.
Обновление re: Комментарий ОП : А как насчет метода собственных значений?
Он работает нормально.
% Using the 5x5 P matrix from above
[V,D] = eigs(double(P.'),1);
Py = -abs(V)/sum(V)
Или, альтернативно,
Pyalt = V./sum(V)
Вы можете проверить это также на 84x84:
[V2,D2] = eigs(double(TPM2'),1);
Py_TPM2 = -abs(V2)/sum(V2);
Py_TPM2alt = V2./sum(V2);
% Evaluate to compare
SteadyStateProb_TPM2(ind).'
Py_TPM2(ind).'
Py_TPM2alt(ind).'
, что дает следующее. Обратите внимание, что состояния [4 8 12 16 20 24].
SteadyStateProb_TPM2(ind).' = 0.0821 0.2052 0.2565 0.2138 0.1336 0.0668
Py_TPM2(ind).' = 0.0821 0.2052 0.2565 0.2138 0.1336 0.0668
Py_TPM2alt(ind).' = 0.0821 0.2052 0.2565 0.2138 0.1336 0.0668