概率生成函数学习笔记
简称 PGF。
哪个天才把多项式出到概率上的啊
性质
设一个多项式 \(F(x)=\sum P(X=i)x^i\)
然后会发现这玩意有一堆性质，然后就被出成题了。。。
- \(F(1)=1\) 。这个显然，所有情况的概率和就是 \(1\) ，\(\sum P(X=i)=1\)
- \(E(X)=F'(1)\) 。
\(E(X)=\sum i\times P(X=i)\) 带系数，怎么搞？求个导，就有 \(i\) 了。
- \(E(X^{\underline k})=F^{(k)}(1)\) （\(k\) 阶导）。
上面那个式子多求几次导即可。
或许有人和我一样不知道 \(E(X^{\underline{k}})\) 是啥意思。。。\(E(X^{\underline{k}})=\sum P(X=i)i^{\underline{k}}\)
- \(E(X^2)=F''(x)+F'(x)\)
\(E(X^2)=E(X^{\underline{2}}+X^{\underline{1}})=E(X^{\underline{2}})+E(X^{\underline{1}})=F''(x)+F'(x)\)
- 方差 \(V(X)=F''(1)+F'(1)-F'(1)^2\)
\(V(x)=E(X^2)-E(X)^2=F''(1)+F'(1)-F'(1)^2\)
- 补充一下 \(V(X)=E(X^2)-E(X)^2\) 是怎么推出来的，我自己想了好久，最后是在别人的提醒下翻了《具体数学》第八章明白的。
主要是上面那步理解了好久，因为 \(E(X)\) 是个常数，所以 \(E(E(X)^2)=E(X)^2\) 被直接提了出来。\(E(X\ E(X))=E(X)E(E(X))=E(X)^2\) ，这个用期望的线性性+乘法分配律即证。
所以，方差等于“随机变量平方的均值减去其均值的平方”
然后就可以做题了。。。
例1
好像是个套路，但是真的nb
令 \([x^n]F(x)\) 为第 \(n\) 轮恰好唱出一个名字的概率，\([x^n]G(x)\) 表示 \(n\) 轮之后还没有唱出名字的概率。
那么
上一轮如果还得接着唱，那么再唱一轮只可能结束或者没结束。后面那个 \(1\) 就是你唱了 \(0\) 轮必然唱不出名字，然而因为 \(G(x)\) 乘了 \(x\) 常数项没了要手动补上。怎么感觉是废话，但是这个式子很有用
首先，根据上面性质那部分推导，我们要求的是 \(F'(1)\) 。
上面那个式子两边求导
直接把 \(x=1\) 带进去
\(F'(1)=G(1)\) 。
所以只要求 \(G(1)\) 就行了。
还有一波骚操作可以导出 \(F(x),G(x)\) 的关系。
考虑如何表示出一个必然结束的情况。
用 \(G(x)\) 来写就是 \(G(x)(\dfrac{x}{n})^{m}\) ，也就是暴力在后面接上一个名字。
但是有可能会存在提前结束的情况。
非常神奇的是用 \(F(x)\) 能写出一个有相同组合意义的式子，所以并不用管提前结束这玩意。
显然，\(G(x)\) 那个式子有可能提前结束当且仅当现在这个串的末尾是原串的 \(\rm Border\) 。
事实上这东西并不严格是 \(\rm Border\) ，因为 \(\rm Border\) 不能等于自身，但是这题可以。能懂就好了qwq
令 \(A_i\) 表示 \([S[1,\cdots,i]=S[m-i+1,m]]\) 也就是是否为 \(\rm Border\) 。
后面那个幂次是补上剩下的操作。其实可以理解成枚举在哪个位子提前结束。
所以现在我们得到的是
带入 \(x=1\)
再用性质里头那个 \(F(1)=1\) ，两边同时乘 \(n^m\) 得到
感觉好牛逼啊，第一次把具体值往生成函数里代诶。
关于怎么求 \(A_i\) ，记住一句曾经在我们班广泛流传的一句话：\(\rm Border\) 的 \(\rm Border\) 是 \(\rm Border\) 。很有用的（
const int N=100005;
#define mod 10000
void fmod(int&x){x-=mod,x+=x>>31&mod;}
int n,m,p[N],a[N],ans;
bool A[N];
void print(int n){
int d[5];d[0]=0;
for(;n;n/=10)d[++d[0]]=n%10;
while(d[0]<4)d[++d[0]]=0;
while(d[0])printf("%c",'0'+d[d[0]--]);
puts("");
}
signed main(){
n=read();
for(int T=read();T;--T){
m=read(),ans=0;
rep(i,1,m)a[i]=read(),A[i]=0;
p[1]=0;
int j=0;
for(int i=2;i<=m;++i){
while(j&&a[j+1]!=a[i])j=p[j];
if(a[j+1]==a[i])++j;
p[i]=j;
}
j=m;
while(j)A[j]=1,j=p[j];
for(int i=1,j=n;i<=m;++i,j=1ll*j*n%mod)if(A[i])fmod(ans+=j);
print(ans);
}
}
例2
设 \([x^n]G(x)\) 表示 \(n\) 轮之后仍未结束的概率。
考虑答案怎么表示。不能像上面只开一个了。但是注意到 \(n\) 很小，于是可以开 \(n\) 个。
\([x^n]F_{i}(x)\) 表示第 \(i\) 个人在第 \(n\) 轮获胜的概率。
和上一题一样，考虑搞一个必然结束的式子出来：
随便在后面接一个串就好了，不妨接上第 \(i\) 个。
对于 \(F_i(x)\) 可以列一个和上面组合意义相同的式子。
其中 \(A_{i,j,k}=[S_i[1,\cdots,k]=S_j[m-k+1,\cdots,m]]\) 可以哈希 \(O(n^3)\) 搞。
还是在枚举在哪个地方出现了提前结束的情况。
所以对于每一个 \(i\in[1,n]\) 都有
还有一个式子是
就是第 \(i - 1\) 轮没有结束，那么第 \(i\) 轮只有没结束和结束两种状态。
我们要求的是所有的 \(F_i(1)\) 。
把 \(x=1\) 带入第二个式子
然后把第一个式子两边带入 \(x=1\) 。
\(n+1\) 个未知数 \(n+1\) 个方程直接高消就好了。
但是 \(2^{300}\) 次方是啥玩意啊。然后我去看了看题解区有什么高妙做法
woc这玩意double居然能存，存完还能高消，怎么这么nb啊，而且我还不会卡诶（
woc为啥题解都能跑，我就被卡精度了啊/kk
woc我为啥 \(2^{300}\) 开的是 int 啊，手残手残，我是傻逼。。。然后就过了qwq
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp(x,y) make_pair(x,y)
#define pb(x) push_back(x)
#define sz(v) (int)v.size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
const int N=305;
#define mod1 998244353
#define mod2 1000000007
#define base1 31
#define base2 37
int n,m,H1[N][N],H2[N][N],pw1[N],pw2[N];
db a[N][N],ans[N],pw[N];
int hash1(int id,int l,int r){
return (H1[id][r]-1ll*H1[id][l-1]*pw1[r-l+1]%mod1+mod1)%mod1;
}
int hash2(int id,int l,int r){
return (H2[id][r]-1ll*H2[id][l-1]*pw2[r-l+1]%mod2+mod2)%mod2;
}
void Gauss(int n){
rep(i,1,n){
int mx=i;
rep(j,i,n)if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
if(mx!=i)rep(j,1,n+1)swap(a[mx][j],a[i][j]);
rep(j,1,n){
if(i==j)continue;
db div=a[j][i]/a[i][i];
rep(k,i,n+1)a[j][k]-=a[i][k]*div;
}
}
rep(i,1,n)ans[i]=a[i][n+1]/a[i][i];
}
signed main(){
n=read(),m=read();
pw[0]=pw1[0]=pw2[0]=1;
rep(i,1,m)pw1[i]=1ll*pw1[i-1]*base1%mod1,pw2[i]=1ll*pw2[i-1]*base2%mod2,pw[i]=2.*pw[i-1];
rep(i,1,n){
static char str[N];scanf("%s",str+1);
rep(j,1,m)H1[i][j]=(1ll*H1[i][j-1]*base1+(str[j]=='H')+1)%mod1,H2[i][j]=(1ll*H2[i][j-1]*base2+(str[j]=='H')+1)%mod2;
}
rep(i,1,n){
db res=0;
rep(j,1,n)rep(k,1,m)
if(hash1(i,1,k)==hash1(j,m-k+1,m)&&hash2(i,1,k)==hash2(j,m-k+1,m))a[i][j]+=pw[k];
a[i][n+1]=-1,a[i][n+2]=0;
}
rep(i,1,n)a[n+1][i]=1;a[n+1][n+1]=0,a[n+1][n+2]=1;
Gauss(n+1);
for(int i=1;i<=n;++i)printf("%.10lf\n",ans[i]);
return 0;
}
预告：可能还会有道例题，等我会做了再放上来。。。
upd：例题没了，因为我不会做。
upd 2021.5.18：真的又做到概率生成函数题了，搬运题解（
例3
设 \(n\) 个生成函数 \(F_i(x)\)，\([x^k]F_i(x)\) 表示第 \(i\) 个人在第 \(k\) 轮获胜的概率。
再设一个生成函数 \(G(x)\)，\([x^k]G(x)\) 表示游戏在第 \(k\) 轮仍未结束的概率。
接下去可以得到若干等式：
即在第 \(i - 1\) 轮没有结束的概率等于在第 \(i\) 轮结束的概率加上没结束的概率，或许写成 \(f_i+g_i=g_{i-1}\) 更好理解。
然后对于每一个 \(1\le i\le n\) 都可以列出下面的等式
其中 \(P_{i,l, r}\) 表示出现 \(S_i[l, r]\) 的概率。
等式左边的含义是给没结束的字符串接上一个字符串让他结束，但是也会存在撞上另一个人的串导致提前结束的情况。
等式右边就是所有可以结束的状况的和。一个串 \(j\) 已经结束之后又接上了 \(S_i[1,l - k]\) 所以除掉 \(P_{i,l - k + 1, l}\)。
可以发现两边组合意义相同。
因为个人习惯，把 \(P\) 的定义改成前缀和，减小一些常数，即 \(P_{i,k}\) 表示出现 \(S_i[1,k]\) 的概率，化简一下可以得到：
然后注意到每个人获胜的概率就是 \(F_i(1)\)，那么把 \(x=1\) 带进去。
\(n+1\) 个方程 \(n+1\) 个未知数直接高斯消元，复杂度 \(O(n^3)\)。
注意到这个方程在所有串都不可能得到的情况下会无解要特判。然后还有一堆奇怪的边界要注意，看个人写法了，主要是除以 \(0\) 的情况要小心。
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp(x,y) make_pair(x,y)
#define pb(x) push_back(x)
#define sz(v) (int)v.size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
const int N = 13;
const int base = 11;
const int mod = 20050103;// tommy0103/qq
//忽然发现不用取模，但是上面那行不想删了（
int n, l, m, s[N][N];
db p[N], P[N][N], a[N][N], ans[N];
LL H[N][N], pw[N];
inline LL getHash(int i, int l, int r) {
return H[i][r] - H[i][l - 1] * pw[r - l + 1];
}
inline void Gauss(int n) {
rep(i, 1, n) {
int mx = i;
rep(j, i + 1, n) if(fabs(a[j][i]) > fabs(a[mx][i])) mx = j;
if(mx != i) { rep(j, 1, n + 1) swap(a[mx][j], a[i][j]); }
rep(j, 1, n) {
if(i == j) continue;
db d = a[j][i] / a[i][i];
rep(k, i, n + 1) a[j][k] -= a[i][k] * d;
}
}
rep(i, 1, n) ans[i] = a[i][n + 1] / a[i][i];
}
signed main() {
n = read(), l = read(), m = read();
rep(i, 1, m) p[i] = read(), p[i] /= read();
pw[0] = 1;
rep(i, 1, l) pw[i] = pw[i - 1] * base;
rep(i, 1, n) {
static char str[N];
scanf("%s", str + 1), P[i][0] = 1;
rep(j, 1, l)
s[i][j] = str[j] - 'A' + 1,
H[i][j] = H[i][j - 1] * base + s[i][j],
P[i][j] = P[i][j - 1] * p[s[i][j]];
}
bool flg = 1;
rep(i, 1, n) {
if(P[i][l] < 1e-11) {
a[i][i] = 1;
continue;
}
flg = 0;
rep(j, 1, n) if(P[j][l] > 1e-11) {
rep(k, 1, l) if(getHash(i, 1, k) == getHash(j, l - k + 1, l))
a[i][j] += 1. / P[i][k];
}
a[i][n + 1] = -1, a[i][n + 2] = 0;
}
if(flg) {
rep(i, 1, n) puts("0.00");
return 0;
}
rep(i, 1, n) a[n + 1][i] = 1;
a[n + 1][n + 1] = 0, a[n + 1][n + 2] = 1;
Gauss(n + 1);
rep(i, 1, n) printf("%.2lf\n", ans[i] < 1e-5 ? 0. : ans[i]);
return 0;
}
