chain.c

解马尔科夫链CHAIN.C

#include <stdio.h>
#include <math.h>


double a[10][10];

void Guass(int n){
	int i,j,k;
	double t;
	for(k=0;k<n-1;k++)
		{t=a[k][k];
		for(j=k;j<n;j++)
			a[k][j]=a[k][j]/t;
		for(i=0;i<n-1;i++)
			if(i!=k)
				{t=a[i][k]/a[k][k];
				for(j=k;j<n;j++)
					a[i][j]=a[i][j]-a[k][j]*t;
				}
		}
	return;
}


void chain(){
	static double p[10][10],pr[10],diff,table[100][10],pnew[10][10],ptemp[10][10],temp[10],exr[10][10];
	int n,i,j,k,s,m,found,inr,inc;
	printf("Welcome to the MARKOV CHAIN ANALYSIS system!\n\n");
	printf("how many states =?        ");
	scanf("%d",&n);
	printf("\nthe steady transmit possibility of step 1 ?\n");
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			scanf("%lf",&p[i][j]);
	printf("\nthe initiate state of step 1 ?\n");
	for(i=0;i<n;i++)
		scanf("%lf",&pr[i]);
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			pnew[i][j]=p[i][j];
	for(i=0;i<n;i++)
		table[0][i]=pr[i];
	printf("\n step 1");
	for(i=0;i<n;i++)
		{printf("\n");
		for(j=0;j<n;j++)
			printf("%f\t",p[i][j]);
		}
	printf("\n");
	for(k=2;k<100;k++)
		{for(j=0;j<n;j++)
			{temp[j]=0;
			for(i=0;i<n;i++)
				temp[j]=temp[j]+pr[i]*pnew[i][j];
			}
		for(i=0;i<n;i++)
			table[k-1][i]=temp[i];
		for(i=0;i<n;i++)
			for(j=0;j<n;j++)
				{ptemp[i][j]=0;
				for(m=0;m<n;m++)
					ptemp[i][j]=ptemp[i][j]+p[i][m]*pnew[m][j];
				}
		for(i=0;i<n;i++)
			for(j=0;j<n;j++)
				pnew[i][j]=ptemp[i][j];
		for(j=0;j<n;j++)
			{for(i=0;i<n-1;i++)
				{
				for(m=i+1;m<n;m++)
					{diff=pnew[i][j]-pnew[m][j];
					if(diff<0)
						diff=-diff;
					if(diff>0.001)
						{found=0;
						break;
						}
					found=1;
					}
				if(diff>0.001) break;
				}
			if(diff>0.0001) break;
			}
		if(found==0)
			{if(k%5==0)
				{printf("\n step %d",k);
				for(i=0;i<n;i++)
					{printf("\n");
					for(j=0;j<n;j++)
						printf("%f\t",pnew[i][j]);
					}
				//getch();
				}
			if(k>=100)
				{printf("\nsteady_state probability have not been detained in 100");
				return;
				}
			}
		else
			{printf("\nstep %d",k);
			for(i=0;i<n;i++)
				{printf("\n");
				for(j=0;j<n;j++)
					printf("%f\t",pnew[i][j]);
				}
			break;
			}
		}
        for(j=0;j<n;j++)
			{temp[j]=0;
			for(i=0;i<n;i++)
				temp[j]=temp[j]+pr[i]*pnew[i][j];
			}
		for(i=0;i<n;i++)
			table[k][i]=temp[i];
		printf("\nThe steady-state probability of being in \n");
		for(j=0;j<n;j++)
			printf("state %d  is %f\n",j,pnew[n-1][j]);
		printf("probability of being in state\n\n");
		for(i=0;i<=k;i++)
			{printf("%d",i);
			for(j=0;j<n;j++)
				printf("\t%f",table[i][j]);
			printf("\n");
			//if(i%10==0)	getch();
			}

		for(s=0;s<n;s++)
			{ inr=0;
			for(j=0;j<n;j++)
				if(j==s)
					continue;
				else
					{inc=0;
					for(i=0;i<n;i++)
						if(i==s)
							continue;
						else
							{
							a[inr][inc]=-p[j][i];
							if(j==i)
								a[inr][inc]=1+a[inr][inc];
							inc++;
							}
					inr++;
					}
			for(i=0;i<n-1;i++)
				a[i][n-1]=1;
			Guass(n);
			i=0;
			for(j=0;j<n;j++)
				if(j!=s)
					exr[j][s]=a[i++][n-1];
				else
					exr[j][s]=1/pnew[n-1][s];
			}
		printf("\nTable of expected first passage times and recurrence times");
		for(i=0;i<n;i++)
			{printf("\n%d",i);
			for(j=0;j<n;j++)
				printf("\t%f",exr[i][j]);
			}
}

void main(){
//	clrscr();
	chain();
//	getch();
}