Bzoj3270:博物馆:概率与期望,高斯消元

本文介绍了一个博物馆问题的求解方法,通过定义状态转移矩阵并利用高斯消元法求解稳态概率。具体地,文章定义了状态转移概率,并考虑了四个主要转移情况,最终通过矩阵运算得出每个状态的长期停留概率。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

题目链接:博物馆

我们用id[i][j]代表一个人到了i另一个人在j的状态

假设id[i][j]代表的状态可以一步走到id[x][y]代表的状态,那么id[x][y]一步也可以走到id[i][j]

所以状态之间的转移形成了一个环,这时候要用高斯消元来解决一下了

设a[id[x][y]][id[i][j]]为从状态id[i][j]转移到id[x][y]的概率来列个方程组

假设现在在id[i][j]这个状态上,接下来有这么几种情况

1:两个人同时停留在原点,a[id[i][j]][id[i][j]]=p[i]*p[j];

2:第一个人走到了x(前提是i->x有边存在),a[id[x][j]][id[i][j]]+=(1-p[i])/du[i]*p[j];

3:第二个人走到了y,a[id[i][y]][id[i][j]]+=p[i]*(1-p[j])/du[j];

4:第一个人走到了x,第二个人走到了y,a[id[x][y]][id[i][j]]+=(1-p[i])/du[i]*(1-p[j])/du[j];

这样我们列出了n*n个方程,形如a[id[x][y]][id[i][j]]=ax+by+cz+...

然后移项即可,只有一个特例即初始点,他的概率要+1,所以移项后常数项为-1,其余的都为0

消元即可

#include<cmath>
#include<queue>
#include<vector>
#include<cstdio>
#include<cstdlib>
#include<iostream>
#include<algorithm>
#define pb push_back
using namespace std;
const int maxn=21;
const int maxm=410;
int n,m,id[maxn][maxn],cnt=0,x,y,du[maxn];
double a[maxm][maxm],equ[maxm],p[maxn];
vector<int>nxt[maxn];

void guass(){
	for (int i=1;i<=cnt;++i){
		int tmp=i;
		for (int j=i+1;j<=cnt;++j)
		    if (abs(a[j][i])>abs(a[tmp][i])) tmp=j;
		if (tmp!=i)
		for (int j=1;j<=cnt+1;++j) swap(a[i][j],a[tmp][j]);
		for (int j=i+1;j<=cnt;++j){
			double N=a[j][i]/a[i][i];
			for (int k=i;k<=cnt+1;++k) a[j][k]-=a[i][k]*N;
		}
	}
	equ[cnt]=a[cnt][cnt+1]/a[cnt][cnt];
	for (int i=cnt-1;i>=1;--i){
		double N=a[i][cnt+1];
		for (int j=i+1;j<=cnt;++j) N-=a[i][j]*equ[j];
		equ[i]=N/a[i][i];
	}
	for (int i=1;i<=n;++i) printf("%.6lf ",equ[id[i][i]]);
}

int main(){
	scanf("%d%d%d%d",&n,&m,&x,&y);
	for (int i=1;i<=m;++i){
		int u,v; scanf("%d%d",&u,&v);
		nxt[u].pb(v); nxt[v].pb(u);
		du[u]++; du[v]++;
	}
	for (int i=1;i<=n;++i) scanf("%lf",&p[i]);
	for (int i=1;i<=n;++i)
		for (int j=1;j<=n;++j) {id[i][j]=++cnt;if(i!=j)a[cnt][cnt]=p[i]*p[j];}
	for (int i=1;i<=n;++i)
	    for (int j=1;j<=n;++j)
	        if (i!=j){
				for (int k=0;k<nxt[i].size();++k){
					int tx=nxt[i][k];
					a[id[tx][j]][id[i][j]]+=(1-p[i])/du[i]*p[j];
				}
				for (int k=0;k<nxt[j].size();++k){
					int ty=nxt[j][k];
					a[id[i][ty]][id[i][j]]+=p[i]*(1-p[j])/du[j];
				}
				for (int k=0;k<nxt[i].size();++k)
				    for (int q=0;q<nxt[j].size();++q){
						int tx=nxt[i][k],ty=nxt[j][q];
						a[id[tx][ty]][id[i][j]]+=(1-p[i])*(1-p[j])/du[i]/du[j];
				    }
	        }
	for (int i=1;i<=n;++i)
	    for (int j=1;j<=n;++j) a[id[i][j]][id[i][j]]-=1.0;
	a[id[x][y]][cnt+1]=-1;
	guass();
}


评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值