洛谷 P3317 [SDOI2014]重建【矩阵树定理】【概率期望】

作者: wjyyy 分类: 概率期望,矩阵树定理,解题报告,高斯消元 发布时间: 2019-06-30 18:27

点击量:29

不是那么简单的矩阵树。

题目描述

T 国有 $N$ 个城市,用若干双向道路连接。一对城市之间至多存在一条道路。

在一次洪水之后,一些道路受损无法通行。虽然已经有人开始调查道路的损毁情况,但直到现在几乎没有消息传回。

幸运的是,此前 T 国政府调查过每条道路的强度,现在他们希望只利用这些信息估计灾情。具体地,给定每条道路在洪水后仍能通行的概率,请计算仍能通行的道路恰有 $N-1$ 条,且能联通所有城市的概率。

输入输出格式

输入格式:

输入的第一行包含整数 $N$。

接下来 $N$ 行,每行 $N$ 个实数,第 $i+1$ 行,$j$ 列的数 $G_{i,j}$ 表示城市 $i$ 与 $j$ 之间仍有道路联通的概率。

输入保证 $G_{i,j}=G_{j,i}$,且 $G_{i,i}=0$;$G_{i,j}$ 至多包含两位小数。

输出格式:

输出一个任意位数的实数表示答案。

你的答案与标准答案相对误差不超过 $10^{-4}$ 即视为正确。

输入输出样例

输入样例#1:

3
0 0.5 0.5
0.5 0 0.5
0.5 0.5 0

输出样例#1:

0.375

说明

$1<N\le 50$

数据保证答案非零时,答案不小于 $10^{-4}$。

题解:

注意到这道题的边是有概率存在。

矩阵树定理可以支持重边,因此它所做的贡献会乘上这条边出现的次数。这一点很容易从矩阵行列式的角度来理解。

不过由于这里是实数,不容易直接联想。我们在其他题目中总是通过取模来减小精度误差,但是我们也可以认为取模后的数变成了另外一个数域。

既然模意义下有意义,那么就可以认为原数代表的东西有意义。

因此我们把实数放在边权上,就可以求出某 $n-1$ 条边同时存在的概率。但是我们也要保证剩下的 $m-n+1$ 条边不存在。因此一棵树出现的概率是 $\prod_{i\in\text{出现的边}}p_i\prod_{i\in\text{没有出现的边}}(1-p_i)$。

前面的 $\prod p_i$ 很容易求出,而后面的并不是非常显然。

考虑到我们知道的信息是有关前面“$i\in\text{出现的边}$”的信息,我们可以通过用全集 $-$ 出现过的边 $=$ 没有出现的边来解决这个问题。

那么可以让前面变成
$$
\prod_{i\in\text{出现的边}}\frac{p_i}{1-p_i}\cdot\prod_{i\in U}(1-p_i)
$$
其中 $U$ 是边的全集。

这样的话,把边权改为 $\frac{p_i}{1-p_i}$ 做矩阵树定理,最后再乘上 $\prod_{i\in U}(1-p_i)$ 就可以了。

然而会出现 $p_i=1$ 的情况,这样的话我们可以把边的两端缩成一个点。不过本题是在实数运算下的,因此可以有一定的误差,这样认为 $p_i=1-\varepsilon$ 就可以了。

发现这样做的话数值差异会比较大,最大可能达到 $\varepsilon^n$。为了进一步减小误差,我们用到高斯消元中的一个技巧,当进行到第 $i$ 行时,可以找第 $i+1$ 到第 $n$ 行中第 $i$ 列元素的绝对值最大的一行与第 $i$ 行交换。

注意矩阵树定理要把基尔霍夫矩阵去掉一行一列。

时间复杂度 $O(n^3)$。

Code:

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define eps 1e-7
double f[55][55];
int main()
{
    int n;
    scanf("%d",&n);
    double P=1;
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n;++j)
        {
            if(i==j)
            {
                scanf("%*lf");
                continue;
            }
            scanf("%lf",&f[i][j]);

            if(fabs(f[i][j]-1)<eps)
                f[i][j]-=eps;
            if(fabs(f[i][j])<eps)
                f[i][j]+=eps;
            if(i<j)
                P*=(1-f[i][j]);
            f[i][j]/=(1-f[i][j]);
            f[i][i]+=f[i][j];
            f[i][j]=-f[i][j];
        }
    double ans=1;
    for(int i=1;i<n;++i)
    {
        for(int j=i+1;j<n;++j)
            if(fabs(f[j][i])>fabs(f[i][i]))
            {
                std::swap(f[j],f[i]);
                ans=-ans;
            }
        double t=f[i][i];
        ans*=t;
        for(int j=i;j<n;++j)
            f[i][j]/=t;
        for(int j=i+1;j<n;++j)
        {
            double r=f[j][i];
            for(int k=i;k<n;++k)
                f[j][k]-=r*f[i][k];
        }
    }
    printf("%.8lf\n",ans*P);
    return 0;
}

说点什么

avatar
  Subscribe  
提醒
/* */