地震后的幻想乡

题意

一幅图的每条边的权值为互不相关的$[0,1]$内的连续随机变量。求图的最小生成树的最大边的权值的期望,结果保留6位小数。

概率密度函数

在研究连续型随机变量的过程中,我们一般用概率密度函数$\rho_a(x)$描述一个随机变量的取值。

$\rho_a(x)$在$x$处的取值越大,说明随机变量$a$有更大的概率取$x$。但是,由于我们研究的$a$是一个连续型随机变量,虽说$a$可能为$x$,但它的值为$x$的概率实际上是$0$,这样也就不太好比较$a$为某个值的概率的大小。因此为了方便,我们用一个更直观的方式表示$a$取值的分布。

$$
P(a\in[l,r]) = \int_l^r\rho_a(x)dx
$$

这也恰好符合”$\rho_a(x)$越大,$a$越容易取得$x$”这一特性。同时,根据这个定义,我们知道,$\int_{-\infty}^{+\infty}\rho_a(x)dx = 1$。

问题规约

求最小生成树的最大边的权值的期望并不是一件容易的事。但是如下的一个问题实际上比较容易解决:原图是否存在一个最大边不大于$x$的生成树。

这个问题相当于:每条边有$x$的概率存在,有$1-x$的概率不存在,求原图联通的概率。利用状态压缩dp可以在$O(n3^n)$的时间复杂度内解决。

设$F[S]$为原图内$S$代表的点集联通的概率。那么我们选定该点集内的一个点为“观察点”,我们尝试观察出这个点集中与这个点联通的点有哪些。如果并不是所有点都与观察点联通,那么这个点集就是不联通的。恰好,每个不连通的情况都不重复地对应着一个观察点所在的联通块。所以我们可以断言,点集$S$联通的概率,等于1减去观察点所在的每一种联通块存在的概率(这个联通块不等于S)。

由于我们只需要求全图联通的概率,我们将观察点设为1号节点即可,这样已经可以补充不漏地计算到所有情况,即:
$$
\forall S\ (1\in S)\ \ : \ \ F[S] = 1-\sum_{T\subsetneqq S, 1\in T}F[T]G(S-T,T)
$$
其中$G(A,B)$表示点集$A,B$之间的所有边都断开的概率。

那么怎么把原问题规约到这个问题上来呢?如果我们枚举最小生成树的最大边的权值$x$,那么每条边出现的概率就是$x$,进而通过上面的状态转移方程求出的每一个状态都会是一个关于$x$的多项式。而多项式$F[S](x)$就表示“S的最小生成树最大边不超过$x$”的概率。

那么,我们又要求$F[s](x)$的期望,怎么办呢?可以这样:

先对$F(x)$求导,$F'(x)$表示最小生成树最大边等于$x$的概率密度函数。为什么呢?反过来考虑:设$f(x)$表示最小生成树最大边的概率密度函数,则$F=\int_0^af(x)dx$表示的是最小生成树最大边不超过$x$的概率。由于$F(0)=0$,所以求导再定积分的结果就是原函数,这样证明是没有毛病的。

则$\int_0^1xF'(x)$为最小生成树最大边期望。为什么呢?当我们在计算期望时,当函数取值为$a$时,它对期望的贡献也为$a$。因此当最大边为$x$时,它的贡献就是$x$,而最大边为$x$的概率又是通过$F'(x)$描述的。所以将两个多项式卷积,意味着点值的直接相乘,这样求出的就是期望。

所以总的来说,我们只需要将状态压缩DP的$F[S]$设为一个多项式即可,原先的乘法、加法、减法全部用多项式的运算代替。最后利用上述方法求出期望就行了。

实现上的问题

在运算的过程中,自变量$x$的系数都是整数,因此我们可以用long long存储。但是在最终求答案时,我们需要将多项式积分,这时会出现小数。因此需要将long long转化为浮点数。在转化过程中会丢失精度,再将转化后的数除以一个较小的整数,这时会严重影响结果。因此我们需要一个技巧避免精度问题:

对于较大的$a$和较小的$b$,如下的方法可以避免精度问题:
$$
\frac{a}{b}=\lfloor\frac{a}{b}\rfloor + \frac{a\mod b}{b}
$$

代码

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <assert.h>
#define mov(x) (1<<(x))

using namespace std;

typedef long long ll;

struct POLY
{
    ll a[50];
    int d;

    POLY () {memset(a, 0, sizeof(a)); d = 0;}
    POLY deriv ()
    {
        POLY c; c.d = d-1;
        for(int i=0; i<c.d; i++) c.a[i] = a[i+1]*(i+1);
        return c;
    }
    POLY operator + (const POLY& t)
    {
        POLY c; c.d = max(d, t.d);
        for(int i=0; i<c.d; i++) c.a[i] = a[i]+t.a[i];
        return c;
    }
    POLY operator - (const POLY& t)
    {
        POLY c; c.d = max(d, t.d);
        for(int i=0; i<c.d; i++) c.a[i] = a[i]-t.a[i];
        return c;
    }
    POLY operator * (const POLY& t)
    {
        POLY c; c.d = d+t.d-1;
        for(int i=0; i<d; i++)
            for(int j=0; j<t.d; j++)
                c.a[i+j] += a[i]*t.a[j];
        return c;
    }
    double f(double x)
    {
        double y = 0, xi = 1;
        for(int i=0; i<d; i++) y += xi*a[i], xi*=x;
        return y;
    }
    void print()
    {
        for(int i=0; i<d; i++) printf("%lld ", a[i]);
        putchar('\n');
    }
};

POLY F[1024], MI[100];
int N, M, E[12], SE[1024][1024], BN[1024];

void dp()
{
    F[1].d = 1;
    F[1].a[0] = 1;
    for(int u=2; u<mov(N); u++)
    {
        if(!(u & 1)) continue;
        for(int s=(u-1)&u; s; s=(s-1)&u)
            if(s & 1)
                F[u] = F[u] + F[s]*MI[SE[s][u^s]];
        for(int i=0; i<F[u].d; i++) F[u].a[i] *= -1;
        F[u].a[0] ++;
    }
}

void calc()
{
    POLY pos = F[mov(N)-1];
    POLY rho = pos.deriv();
    POLY mul; mul.d = 2, mul.a[0] = 0, mul.a[1] = 1;
    POLY ans = rho*mul;
    long double prt = 0;
    for(int i=0; i<ans.d; i++) prt += (long double)(ans.a[i]%(i+1))/(long double)(i+1);
    while(prt < 0) prt++;
    while(prt >= 1) prt--;
    printf("%.6LF\n", prt);
}

void init()
{
    MI[0].d = 1;
    MI[0].a[0] = 1;
    MI[1].d = 2;
    MI[1].a[0] = 1, MI[1].a[1] = -1;
    for(int i=2; i<=M; i++) MI[i] = MI[i-1]*MI[1];
    for(int i=1; i<mov(N); i++) BN[i] = BN[i>>1] + (i&1);
    for(int s=1; s<mov(N); s++)
    {
        int u = (~s)&(mov(N)-1);
        for(int t=u; t; t=(t-1)&u)
            for(int i=0; i<N; i++)
                if(s & mov(i))
                    SE[s][t] += BN[E[i]&t];
    }
}

void input()
{
    scanf("%d%d", &N, &M);
    for(int i=0; i<M; i++)
    {
        int a, b;
        scanf("%d%d", &a, &b);
        E[a-1] |= mov(b-1);
        E[b-1] |= mov(a-1);
    }
}

int main()
{
    input();
    init();
    dp();
    calc();
    return 0;
}
分类: 文章

发表评论

电子邮件地址不会被公开。 必填项已用*标注

你是机器人吗? =。= *