03-markov-chain.sce

// Однородная цепь Маркова с дискретным временем задана графом состояний. 
// На нулевом шаге она всегда находится в первом состоянии.
// 1. Найти предельные вероятности состояний этой цепи Маркова, если они существуют, решив систему линейных алгебраических уравнений.
// 2. Найти распределение вероятностей состояний цепи через достаточно большое количество шагов путём возведения в соответствующую степень матрицы переходных вероятностей. Сравнить его с предельными вероятностями состояний, если они существуют.
// 3. Статистически оценить вероятности состояний через достаточно большое количество шагов методом Монте-Карло, моделируя поведение цепи достаточно большое количество раз. Сравнить их с предельными вероятностями состояний, если они существуют.

k = 1000; // количество шагов
experiments_number = 1000;

p = [0.0, 0.8, 0.2, 0.0; 1.0, 0.0, 0.0, 0.0; 0.0, 0.9, 0.0, 0.1; 0.0, 0.8, 0.2, 0.0]; // матрица переходных вероятностей; ; -- разделяет строки матрицы

[n, m] = size(p); // n -- количество строк, m -- столбцов; [n1, n2, ...] = size(x), функция возвращает в каждом аргументе значени соответствующего измерения матрицы

// Однородная цепь Маркова с дискретным временем задана графом состояний.
// НУЛЕВОМ ШАГЕ ОНА ВСЕГДА НАХОДИСЯ В ПЕРВОМ СОСТОЯНИИ.
// столбец из n элементов, где 
// первый -- 1, 
// затем все остальные -- 0
// zeros(m1, m2) для создания матрицы размерами (m1, m2)
p_0 = [1; zeros(n - 1, 1)]; // распределение вероятностей состояний на нулевом шаге; записано в виде столбца

// 1) найти предельные вероятности состояний этой цепи Маркова, если они существуют, решив систему линейных алгебраических уравнений:
// решение системы линейных алгебраических уравнений вида Ax = b;
// 1) для значений A формируется матрица коэффициентов при неизвестных, каждая строка которой содержит коэффициенты одного уравнения,
// [p' - eye(n, m); ones(1, m)]
// 2) для значений b формируется вектор-столбец из свободных коэффициентов;
// []zeros(n, 1); -1.0] 
    
// матрица A:
// 1) p' -- транспонированная матрица переходных вероятностей; ' -- оператор транспонирования;
// 2) eye(n, m) -- единичная матрица, n строк, m столбцов;
// 3) ones(1, m) -- нижняя строка -- единицы

// матрица B:
// 1) столбец из n + 1 элементов,
// 2) последний элемент -- -1 
    
// функция возвращает найденные значения неизвестных системы:
function calc_lim_probabilities(p)
    [x0, ker] = linsolve([p' - eye(n, m); ones(1, m)], [zeros(n, 1); -1.0]);
    str = msprintf("предельные вероятности состояний данной цепи Маркова:");
    disp(str);
    disp(x0); // предельные вероятности состояний данной цепи Маркова
endfunction

// 2. Найти распределение вероятностей состояний цепи через достаточно большое количество шагов путём возведения в соответствующую степень матрицы переходных вероятностей.
// 1) p_k -- распределение вероятностей состояний цепи через количество шагов k,
// 2) p_0' -- распределение вероятностей состояний на нулевом шаге,
// 3) p -- матрица переходных вероятностей
function calc_probabilities(p, p_0, k)
    p_k = p_0' * (p ^ k);
    str = msprintf("распределение вероятностей состояний цепи через достаточно большое количество шагов путём возведения в соответствующую степень матрицы переходных вероятностей на шаге №%d:", k); // записывает данные в строку
    disp(str);
    disp(p_k'); // выводит переменные
endfunction

function ans=generate_discrete(p) // генерируем, куда будем переходить
    g = rand(); // равномерное распределение
    ans = 0;
    s = 0.0;
    while s < g,
        ans = ans + 1;
        s = s + p(ans);
    end
endfunction

// 3. Статистически оценить вероятности состояний через достаточно большое количество шагов методом Монте-Карло, 
// моделируя поведение цепи достаточно большое количество раз. 
// Сравнить их с предельными вероятностями состояний, если они существуют.
// Суть метода заключается в следующем: 
// 1) процесс моделируется при помощи генератора случайных величин,
// 2) это повторяется много раз,
// 3) а потом на основе полученных случайных данных вычисляются вероятностные характеристики решаемой задачи.
function estimate_probabilities(p, p_0, k, experiments_number)
    p_k = zeros(n, 1);              // счётчик; столбец длины n из нулей
    for i = 1 : experiments_number  // делаем 1000 (experiments_number) экспериментов
        state = 1;
        for j = 1 : k,              // в каждом эксперименте делаем 1000 (k) шагов
            state = generate_discrete(p(state, :)) // берём state-овскую строчку
        end                         // смотрим где остановились (номер состояния)
        
		// в p_k, в котором указано сколько раз мы остановливались в state-овской вершине через k (1000) шагов
        p_k(state) = p_k(state) + 1; // увеличиваем счётчик в состоянии (state) в котором остановились; (state) -- обратимся к элементу с индексом state
    end
    p_k = p_k / experiments_number; // получаем статистическую оценку; делится каждый элемент ВЕКТОРА
    str = msprintf("статистическая оценка вероятностей (через достаточно большое количество шагов) на шаге №%d:", k);
    disp (str);
    disp(p_k);
endfunction

calc_lim_probabilities(p);				// 1. Найти предельные вероятности состояний этой цепи Маркова, если они существуют, решив систему линейных алгебраических уравнений.
calc_probabilities(p, p_0, k);				// 2. Найти распределение вероятностей состояний цепи через достаточно большое количество шагов путём возведения в соответствующую степень матрицы переходных вероятностей. Сравнить его с предельными вероятностями состояний, если они существуют.
estimate_probabilities(p, p_0, k, experiments_number);	// 3. Статистически оценить вероятности состояний через достаточно большое количество шагов методом Монте-Карло, моделируя поведение цепи достаточно большое количество раз. Сравнить их с предельными вероятностями состояний, если они существуют.