概率生成函数学习笔记

前言

概率生成函数好像是个很厉害的东西啊……如果有掷骰（tou）子的问题似乎可以直接套板子的说……

本篇文章全部都是抄《浅谈生成函数在掷骰子问题上的应用》（杨懋龙）这篇论文的

定义

我们定义一个形式幂级数(A(x))，称它为离散随机变量(X)的概率生成函数，当且仅当对于(A(x))的每一项(a_i)，都有(a_i=P(X=i))

性质

容易发现以下几个性质

1.$$A(1)=sum_{i=0}^infty P(X=i)=1$$

2.$$A'(1)=sum_{i=0}^infty iP(X=i)1^{i-1}=E(X)$$

进一步推倒可以得出

[E(x^underline{k})=F^{(k)}(1),k eq 0 ]

其中(F^{(k)}(x))表示(F(x))的(k)阶导

[DX=F''(1)+F'(1)-(F'(1))^2 ]

其中(DX)表示随机变量(X)的方差

最后一个证明一下好了因为我也不太会方差这个东西……

方差的定义为$$DX=E(X-EX)^2$$，即括号里那个东西的平方的期望

根据期望的线性我们可以推倒

[egin{aligned} DX &=E(X-EX)^2\ &=E(X^2-2XEX+(EX)^2)\ &=E(X^2)-2EX imes EX+(EX)^2\ &=E(X^2)-(EX)^2 end{aligned} ]

把上面的柿子代入发现成立

关于这个东西怎么用……还是拿几道题来讲一下好了……

例题

洛谷P4548 [CTSC2006]歌唱王国

题意：给定一个长度为(L)的序列(A)。然后每次掷一个标有(1)到(m)的公平骰子并将其上的数字加入到初始为空的序列(B)的末尾，如果序列B中已经出现了给定序列(A)，即(A)是(B)的子串，则停止，

求序列(B)的期望长度。(L ≤ 10^5)

我们定义一个字符串的前缀(S[1,i])为这个字符串的(border)当且仅当(S[1,i]=S[L-i+1,L])，其中(L)为串长

定义(a_i)，当且仅当(S[1,i])是(S)的(border)时(a_i)为(1)，否则(a_i)为(0)

定义答案的概率生成函数为(F(x))，即(f_i)表示掷了(i)次骰子游戏结束的概率，以及(G(x))，(g_i)表示掷了(i)次骰子游戏仍未结束的概率

那么容易发现两个性质

[G(x)+F(x)=1+xG(x) ]

也就是说(g_i=g_{i+1}+f_{i+1})，即如果第(i)次未结束，那么第(i+1)次只有结束或未结束，(+1)是因为常数项

[G(x)left({1over m}x ight)^L=sum_{i=1}^La_iF(x)left({1over m}x ight)^{L-i} ]

即如果我们在一个未结束的串后面加上整个(A)肯定结束，然而还有可能没有加完整个串就已经结束了。通过分析可知，如果我们加了(A)的前(i)个字符之后结束，即有(A[L-i+1,L]=A[1,i])，那么根据定义，(A[1,i])是一个(border)，然后再把剩下的(L-i)个字符加进去就行了

对于(1)式，我们对两边求导，再把(x=1)代入，得

[G'(x)+F'(x)=G(x)+xG'(x) ]

[F'(1)=G(1) ]

即我们所需要求的(E(x)=F'(1)=G(1))

对于(2)式，我们把(1)代入，在两边同乘上(m^L)，得

[G(1)=sum_{i=1}^La_im^iF(1) ]

又因为(F(1)=1)，最终可以化作

[E(x)=sum_{i=1}^La_im^i ]

而对于(a_i)我们是可以(O(L))求出来的，直接哈希或者跑(kmp)就行了（代码里写的是(kmp)）

于是复杂度为(O(L))

//minamoto
#include<bits/stdc++.h>
#define R register
#define fp(i,a,b) for(R int i=(a),I=(b)+1;i<I;++i)
#define fd(i,a,b) for(R int i=(a),I=(b)-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
char buf[1<<21],*p1=buf,*p2=buf;
inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
int read(){
    R int res,f=1;R char ch;
    while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1);
    for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0');
    return res*f;
}
char sr[1<<21],z[20];int C=-1,Z=0;
inline void Ot(){fwrite(sr,1,C+1,stdout),C=-1;}
void print(R int x){
    if(C>1<<20)Ot();if(x<0)sr[++C]='-',x=-x;
    while(z[++Z]=x%10+48,x/=10);
    while(sr[++C]=z[Z],--Z);sr[++C]='
';
}
const int N=1e5+5,P=1e4;
inline int add(R int x,R int y){return x+y>=P?x+y-P:x+y;}
inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}
int bin[N],kmp[N],a[N],n,p,res,pos;
int main(){
//	freopen("testdata.in","r",stdin);
	p=read(),bin[0]=1;fp(i,1,1e5)bin[i]=mul(bin[i-1],p);
	for(int T=read();T;--T){
		n=read();fp(i,1,n)a[i]=read();
		kmp[0]=kmp[1]=0;
		for(R int i=2,j=0;i<=n;++i){
			while(j&&a[j+1]!=a[i])j=kmp[j];
			j+=(a[j+1]==a[i]),kmp[i]=j;
		}
		pos=n,res=0;
		while(pos)res=add(res,bin[pos]),pos=kmp[pos];
		printf("%04d
",res);
	}
	return 0;
}

洛谷P3706 [SDOI2017]硬币游戏

题意：给定(n)个长度为(m)的(01)序列(A_i)，(m)个序列互不相同，有一个初始为空的序列(B)，每次等概率生成(0)或(1)并将其放入序列(B)，若这一过程中(A_i)最先作为(B)的子串出现则称(i)获胜，请对(forall iin [1,m])求出(i)获胜的概率

因为证明基本和上面差不多，下面我就不给证明直接放柿子了

定义(P(A_i)=prod_{iin A_i}P_i)

定义(a_{i,j,k})，当且仅当(A_i[1,k]=A_j[m-k+1,m])时值为(1)否则为(0)，可以用(hash)从而在(O(n^3))的时间内解得

定义(f_{i,j})表示首次出现的序列是(A_i)且随机序列长度为(j)的概率，(F_i(x))为其生成函数，定义辅助序列(g_i)表示随机序列长度为(i)时仍未结束的概率，生成函数为(G(x))

容易得到

[G(x)+sum_{i=1}^nF_i(x)=1+xG(x) ]

[G(x)P(A_i)x^m=sum_{j=1}^nsum_{k=1}^ma_{i,j,k}F_j(x)P(A_i[k+1,m])x^{L_i-k} ]

前一个柿子这里就不用管了，我们考虑后一个柿子，把(x=1)代入，可以得到

[G(1)=sum_{j=1}^nsum_{k=1}^ma_{i,j,k}F_j(1)P({1over A_i[1,k]}) ]

对于每一个(i)都有这么一个方程，我们需要解出(F_i(1))和(G(1))，那么总共有(n)个方程和(n+1)个变量

等会儿好像还是不能解啊……

我们再转过头来看看……(F_i(1))表示第(i)个人获胜的概率……那么似乎有

[sum_{i=1}^nF_i(1)=1 ]

这样就有(n+1)个方程了，高斯削元就是了，时间复杂度(O(n^3))

//minamoto
#include<bits/stdc++.h>
#define R register
#define fp(i,a,b) for(R int i=(a),I=(b)+1;i<I;++i)
#define fd(i,a,b) for(R int i=(a),I=(b)-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
const int N=305,P=1e9+7;const double eps=1e-10;
double mp[N][N],b[N];char s[N];int bin[N],h[N][N],n,m;
inline int Hash(R int i,R int l,R int r){return ((h[i][r]-1ll*h[i][l-1]*bin[r-l+1])%P+P)%P;}
void Gauss(int n){
	fp(i,1,n){
		if(mp[i][i]>-eps&&mp[i][i]<eps){
			fp(j,i+1,n)if(mp[j][i]<-eps||mp[j][i]>eps){
				fp(k,i,n+1)swap(mp[i][k],mp[j][k]);
				break;
			}
		}
		double t=1.0/mp[i][i];fp(j,i,n+1)mp[i][j]*=t;
		fp(j,i+1,n){
			t=mp[j][i];
			fp(k,i,n+1)mp[j][k]-=mp[i][k]*t;
		}
	}
	fd(i,n-1,1)fp(j,i+1,n)mp[i][n+1]-=mp[j][n+1]*mp[i][j];
}
int main(){
//	freopen("testdata.in","r",stdin);
	scanf("%d%d",&n,&m);
	bin[0]=b[0]=1;
	fp(i,1,m)bin[i]=(bin[i-1]<<1)%P,b[i]=b[i-1]*2;
	fp(i,1,n){
		scanf("%s",s+1);
		fp(j,1,m)h[i][j]=((h[i][j-1]<<1)+(s[j]=='H'))%P;
	}
	fp(i,1,n){
		fp(j,1,n)fp(k,1,m)(Hash(i,1,k)==Hash(j,m-k+1,m))?mp[i][j]+=b[k]:0;
		mp[i][n+1]=-1;
	}
	fp(i,1,n)mp[n+1][i]=1;mp[n+1][n+2]=1;
	Gauss(n+1);
	fp(i,1,n)printf("%.8lf
",mp[i][n+2]);
	return 0;
}