среда, 14 ноября 2012 г.

Цепь Маркова

    Рассмотрим реализацию цепи Маркова на 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)}$, мы в принципе обладаем исчерпывающем описанием процесса.

    Для компьютерной симуляции нам также потребуются две важные функции: функция инициации и функция обновления. Определяются они довольно просто:

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++ это может выглядеть так:

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

1 комментарий: