Рассмотрим реализацию цепи Маркова на 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>). Если цепь является неоднородной, то матрица переходных вероятностей на каждом шаге меняется - единственное условие: она должна быть стохастической.
Much wow
ОтветитьУдалить