概率生成函數學習筆記

前言

概率生成函數好像是個很厲害的東西啊……如果有擲骰（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\neq 0 \]

其中\(F^{(k)}(x)\)表示\(F(x)\)的\(k\)階導

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

其中\(DX\)表示隨機變量\(X\)的方差

最后一個證明一下好了因為我也不太會方差這個東西……

方差的定義為$$DX=E(X-EX)^2$$，即括號里那個東西的平方的期望

根據期望的線性我們可以推倒

\[\begin{aligned} DX &=E(X-EX)^2\\ &=E(X^2-2XEX+(EX)^2)\\ &=E(X^2)-2EX\times 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({1\over m}x\right)^L=\sum_{i=1}^La_iF(x)\left({1\over m}x\right)^{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]='\n';
}
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\n",res);
	}
	return 0;
}

洛谷P3706 [SDOI2017]硬幣游戲

題意：給定\(n\)個長度為\(m\)的\(01\)序列\(A_i\)，\(m\)個序列互不相同，有一個初始為空的序列\(B\)，每次等概率生成\(0\)或\(1\)並將其放入序列\(B\)，若這一過程中\(A_i\)最先作為\(B\)的子串出現則稱\(i\)獲勝，請對\(\forall i\in [1,m]\)求出\(i\)獲勝的概率

因為證明基本和上面差不多，下面我就不給證明直接放柿子了

定義\(P(A_i)=\prod_{i\in 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}^n\sum_{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}^n\sum_{k=1}^ma_{i,j,k}F_j(1)P({1\over 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\n",mp[i][n+2]);
	return 0;
}