ma er ke fu.cpp

该程序完成马尔可夫分析

#include<iostream.h>
#include<math.h>
#include<iomanip.h>
int n;
double pr[10],p[10][10],pnew[10][10],ptemp[10][10],
       table[100][10],extim[10][10],a[10][10];
double min=1e-5;
void mul()     //矩阵乘法
{
	for(int i=0;i<=n;i++)
		for(int j=0;j<=n;j++)
		{
			pnew[i][j]=0;
			for(int k=0;k<=n;k++)			
				pnew[i][j]+=ptemp[i][k]*p[k][j];
		}
}

void state(int k)  //计算状态分布
{
	for(int i=0;i<=n;i++)
	{
		table[k][i]=0;
		for(int j=0;j<=n;j++)
			table[k][i]+=pr[j]*ptemp[i][j];
	}
}

double hope()      //计算期望值，以确定是否达到稳态
{
	double dis=0;
	int np=n+1;
	for(int i=0;i<=n;i++)
		for(int j=0;j<=n;j++)
			dis+=(pnew[i][j]-ptemp[i][j])*(pnew[i][j]-ptemp[i][j]);
	return sqrt(dis/(np*np));
}

void print()      //输出稳态概率
{
    
	for(int i=0;i<=n;i++)
	{
		for(int j=0;j<=n;j++)
			cout<<setw(8)<<ptemp[j][i]<<"   ";
		cout<<endl;
	}
}

void cass()       //高斯消元法
{
	int i,j,k;
	for(k=0;k<n-1;k++) 
		for(i=k+1;i<n;i++)  
		{   
			double client = a[i][k]/a[k][k];   
			for(j=k+1;j<n;j++)    
				a[i][j]=a[i][j]-client*a[k][j];   
			a[i][n]=a[j-1][n]-client*a[k][n];   
		} 
		a[n-1][n]=a[n-1][n]/a[n-1][n-1]; 
		for(i=n-2;i>=0;i--) 
		{ 
			double temp=0;
			for (j=i+1;j<n;j++)   
				temp+=a[i][j] * a[j][n];  
			a[i][n]=(a[i][n]-temp)/a[i][i];
		}
}

void main()
{
	int np;
	cout<<"输入状态数:"<<endl;
	cin>>np;
	n=np-1;
	cout<<"输入初始各状态的概率分布:"<<endl;
	for(int i=0;i<=n;i++)
		cin>>pr[i];
	cout<<"输入 1 步平稳转移概率:"<<endl;
	for(int j=0;j<=n;j++)
		for(i=0;i<=n;i++)
			cin>>p[i][j];
	for(i=0;i<=n;i++)
		for(j=0;j<=n;j++)
			ptemp[i][j]=p[i][j];
	for(i=0;i<=n;i++)
		table[1][i]=pr[i];
	for(int k=2;k<=100;k++)
	{
		mul();
		if(hope()<=min)
			break;
	    for(i=0;i<=n;i++)		
			for(j=0;j<=n;j++)			
				ptemp[i][j]=pnew[i][j];
		state(k);
		if(k%5==0)
		{		
			cout<<k<<" 步平稳转移概率: "<<endl;		
			print();        
			cout<<endl;
		}
	}
	cout<<"各步的概率分布:"<<endl;
	for(i=1;i<k;i++)
	{
		cout<<i<<": ";
		for(j=0;j<=n;j++)
			cout<<setw(8)<<table[i][j]<<"   ";
		cout<<endl<<endl;
	}
	cout<<"稳态概率:"<<endl;
	for(i=0;i<=n;i++)
		cout<<table[k-1][i]<<"    ";
	cout<<endl;
	for(i=0;i<=n;i++)
		extim[i][i]=1/table[k-1][i];
    for(i=0;i<=n;i++)
	{
		for(j=0;j<=n;j++)
			for(k=0;k<=n;k++)
				a[j][k]=p[k][j];		
		for(k=0;k<=n;k++)
			for(j=i;j<n;j++)
				a[k][j]=a[k][j+1];
		for(k=i;k<n;k++)
			for(j=0;j<=n;j++)
				a[k][j]=a[k+1][j];
		for(k=0;k<n;k++)
			for(j=0;j<n;j++)
				a[k][j]*=-1;
		for(k=0;k<n;k++)
			a[k][k]+=1;
		for(k=0;k<n;k++)
			a[k][n]=1;
		cass();
		j=0;
		for(k=0;k<=n;k++)
			if(k!=i)
			{
				extim[k][i]=a[j][n];
				j++;
			}
	}
	cout<<"期望首次到达时间和再现时间:"<<endl;
	for(k=0;k<=n;k++)
	{
		for(j=0;j<=n;j++)
			cout<<extim[k][j]<<"    ";
		cout<<endl;
	}
}