Рассмотрим реализацию цепи Маркова на C++ с использованием STL. Сначала некоторые сведения из теории.
Нас будут интересовать дискретные цепи Маркова. Для понимания достаточно вводного изложения из Википедии. Будем рассматривать последовательность случайных величин $X_{i}$, которые принимают значения в некотором пространстве состояний $S=\{s_{1}...s_{n}\}$. Важное свойство нашего процесса - значение процесса на шаге $X_{i}$ зависит только от состояния процесса на предыдущем шаге $X_{i-1}$, или более формально - $P(X_{n+1}=s_{j}|X_{0}=s_{k}..X_{n}=s_{i})=P(X_{n+1}=s_{j}|X_{n}=s_{i})$.
Так как мы рассматриваем процессы с конечным числом состояний, то переходные вероятности между состояниями можем задать в виде матрицы размерности $N*N$, где $N=|S|$. Далее, если переходные вероятности для состояний постоянны и не зависит от номера текущего шага, то данная матрица называется однородной. Задав вектор распределения начальных состояний, скажем, $\mu^{(0)}$, мы в принципе обладаем исчерпывающем описанием процесса.
Для компьютерной симуляции нам также потребуются две важные функции: функция инициации и функция обновления. Определяются они довольно просто:
Нас будут интересовать дискретные цепи Маркова. Для понимания достаточно вводного изложения из Википедии. Будем рассматривать последовательность случайных величин $X_{i}$, которые принимают значения в некотором пространстве состояний $S=\{s_{1}...s_{n}\}$. Важное свойство нашего процесса - значение процесса на шаге $X_{i}$ зависит только от состояния процесса на предыдущем шаге $X_{i-1}$, или более формально - $P(X_{n+1}=s_{j}|X_{0}=s_{k}..X_{n}=s_{i})=P(X_{n+1}=s_{j}|X_{n}=s_{i})$.
Так как мы рассматриваем процессы с конечным числом состояний, то переходные вероятности между состояниями можем задать в виде матрицы размерности $N*N$, где $N=|S|$. Далее, если переходные вероятности для состояний постоянны и не зависит от номера текущего шага, то данная матрица называется однородной. Задав вектор распределения начальных состояний, скажем, $\mu^{(0)}$, мы в принципе обладаем исчерпывающем описанием процесса.
Для компьютерной симуляции нам также потребуются две важные функции: функция инициации и функция обновления. Определяются они довольно просто:
1) Функция инициации $\psi :[0,1]-S$. Зададим ее в виде отображения, которое ставит в соответствию некоторому числу из промежутка [0,1] номер состояния в виде $s_{i}$ если некоторое сгенерированное нами число $x$ принадлежит интервалу $(\sum_{j=0}^{i-1} \mu^{(0)}(s_{j});\sum_{j=0}^{i} \mu^{(0)}(s_{j}))$. Очевидно, что число $x$ принадлежит промежутку от нуля до единицы.
2) Функция обновления $\phi :s\times[0,1]$. Зададим ее в виде отображения, которое ставит в соответствию некоторой паре чисел $x, s_{i}$ номер следующего состояния в виде $s_{j}$ если некоторое сгенерированное нами число $x$ принадлежит интервалу $(\sum_{l=0}^{j-1} p_{ij};\sum_{l=0}^{j} p_{ij})$.
Проще говоря, на первом шаге мы подставляем в функцию инициации некоторое случайное число от нуля до единицы и получаем начальное состояние (состояния я нумерую начиная с нуля), а затем состояние на предыдущем шаге и новое случайное число итеративно подставляются в функцию обновления. На C++ это может выглядеть так:
2) Функция обновления $\phi :s\times[0,1]$. Зададим ее в виде отображения, которое ставит в соответствию некоторой паре чисел $x, s_{i}$ номер следующего состояния в виде $s_{j}$ если некоторое сгенерированное нами число $x$ принадлежит интервалу $(\sum_{l=0}^{j-1} p_{ij};\sum_{l=0}^{j} p_{ij})$.
Проще говоря, на первом шаге мы подставляем в функцию инициации некоторое случайное число от нуля до единицы и получаем начальное состояние (состояния я нумерую начиная с нуля), а затем состояние на предыдущем шаге и новое случайное число итеративно подставляются в функцию обновления. На C++ это может выглядеть так:
#include "stdafx.h"
#include <iostream>
#include <vector>
#include <assert.h>
#include "MyMatrix.h" //тут используем пользовательский класс матрицы
#include <ctime>
using std::vector;
vector<double> genCumProb(vector<double>& initState);
size_t initiation(double x, vector<double>& initState);
size_t update(int s,double x, MyMatrix& m); //функция обновления
double unifRand()
{
return std::rand() / double(RAND_MAX);
}
double unifRand(double a, double b)
{
return (b-a)*unifRand() + a;
}
void seed()
{
srand(time(0));
}
template <typename T> void print(vector<T>& initState){
vector<T>::iterator it;
for (it=initState.begin();it<initState.end();it++){
std::cout<<*it;
std::cout<<std::endl;
}
}
int _tmain(int argc, _TCHAR* argv[])
{
seed();
vector<double> mu;
mu.push_back(0.1);
mu.push_back(0.2);
mu.push_back(0.1);
mu.push_back(0.2);
mu.push_back(0.1);
mu.push_back(0.05);
mu.push_back(0.1);
mu.push_back(0.05);
mu.push_back(0.05);
mu.push_back(0.05);
std::cout<<"Initial state probabilities: "<<std::endl;
std::cout<<std::endl;
print(mu);
std::cout<<std::endl;
vector<double> cP=genCumProb(mu);
MyMatrix m=MyMatrix(10,10);
m.Mfill();
for (size_t i=0;i<10;i++){
for (size_t j=0;j<10;j++)
m.at(i,j)=0.1;}
std::cout<<std::endl;
std::cout<<"Transition probablities matrix: "<<std::endl;
std::cout<<std::endl;
std::cout<<m<<std::endl;
std::cout<<std::endl;
vector<int> statesHistory;
std::cout<<"Initial chain state: "<<std::endl;
double x=unifRand(0,1);
statesHistory.push_back(initiation(x,cP));
std::cout<<statesHistory.at(0)<<std::endl;
std::cout<<"With initial random value of "<<x<<std::endl;
int s=statesHistory.at(0);
for (int j=0;j<1000;j++){
statesHistory.push_back(update(s,unifRand(0,1),m));
s=statesHistory.back();
std::cout<<s<<" ";
}
system("pause");
return 0;
}
size_t initiation(double x, vector<double>& cumProb){
size_t ind=0;
if (x<cumProb.at(0)) return 0;
for (size_t i=0;i<cumProb.size()-1;i++){
if (x>=cumProb.at(i)&&x<cumProb.at(i+1)) ind=i+1;
}
return ind;
}
vector<double> genCumProb(vector<double>& initState){
assert(initState.size()!=0);
vector<double> cumProb(initState);
double temp=0;
for (vector<double>::iterator it=cumProb.begin()+1;it<cumProb.end();it++)
{temp=*(it-1);
(*it)+=temp;}
return cumProb;
};
size_t update(int s,double x, MyMatrix& m){
vector<double> tempCum;
double temp=0;
tempCum.push_back(m.at(s,0));
for (size_t i=1;i<m.columns();i++){
temp=tempCum.at(i-1);
tempCum.push_back(m.at(s,i)+temp);
temp+=m.at(s,i);
}
size_t ind=0;
if (x<tempCum.at(0)) return 0;
for (size_t i=0;i<tempCum.size()-1;i++){
if (x>=tempCum.at(i)&&x<tempCum.at(i+1)) ind=i+1;
}
return ind;
}
При необходимости легко заменить <double> на шаблон и использовать соответствующим образом (для первоначальных целей нужен был именно <double>). Если цепь является неоднородной, то матрица переходных вероятностей на каждом шаге меняется - единственное условие: она должна быть стохастической.
