概率生成函数-掷骰子问题详解

概率生成函数-掷骰子问题详解

做SDOI2007，第5题啃不动了来看论文（杨懋龙）

概率生成函数

[F(x) = sum_{i=0}^{inf}Pr(X=i)x^i ]

性质

[F(1)=sum_{i=0}^{inf}Pr(X=i) =1 ]

​ 相当于所有情况的概率总和

[F'(x)=sum_{i=0}^{inf}i imes Pr(X=i)x^{i-1} ]

[E(x)=F'(1)=sum_{i=1}^{inf}iPr(X=i) ]

​ (常数项被求导搞掉了，不过求期望可以不管)

还有些不知道咋用的结论

[Eleft(X^{underline k} ight)=F^{(k)}(1), k eq 0 ]

​ 就是求导后前面的系数成下降幂的形式了，再乘上概率

还有个随机变量的方差

因为$D(x)=E(X^2)-(E(X))^2 $

上面的证明：

[egin{aligned}n imes D(x) &=sum_{i=1}^{n}left(x_{i}-E(x) ight)^{2}=sum_{i=1}^{n}left(x_{i}^{2}-2 x_{i} E(x)+E(x)^{2} ight) \&=sum_{i=1}^{n} x_{i}^{2}-2 E(x) sum_{i=1}^{n} x_{i}+E(x)^{2} imes n \=& n Eleft(x^{2} ight)-2 E(x) imes(n E(x))+n E(x)^{2} \&=n Eleft(x^{2} ight)-n E(x)end{aligned} ]

所以(D(x)=E(x*(x-1))+E(x)-(E(x))^2=F''(1)+F'(1)-(F'(1))^2)

掷骰子题型

大概是一直持续一个游戏，直到满足某个条件获胜终止，求一个东西的期望（一般是持续时间）

例题：CTST2006歌唱王国

题意：给一个长L值域为m的序列A。每个时间掷一个1∼m的公平骰子并将这个数字加入到初始为空的序列B的末尾，当A是B子串时，停止，求期望停止时间。((n,mleq 1e5))

分析

令 (f_i) 为结束时随机序列长度为 (i) 的概率，其概率生成函数为 (F(x)) 。令辅助数列 (g_i) 为随机序列长度达到i且还未结束的概率，其普通生成函数为 (G(x)) 。

考虑下一个骰子，有 (F(x)+G(x)=1+G(x)*x)

然后考虑两种情况，第一种是在每一个未结束的序列上直接接上长为L的序列A，但是这样会有很多多出来的情况，比如A是(12312)，此时的B是(123'12312')，引号部分是接上的，可见只接前两个数字就已经匹配，相当于被重复统计，此时把(k+i)项结束的概率的项数加上(L-i)，来对应第一种情形，并乘上相应概率即可。

不难发现此时一定接上的是一个border，即‘12’，设(a_i=[A[1...i]是border])，那么我们有

[G(x) cdotleft(frac{1}{m} x ight)^{L}=sum_{i=1}^{L} a_{i} cdot F(x) cdotleft(frac{1}{m} x ight)^{L-i} ]

对应两种情况

由上文，(E(x)=F'(1))，考虑求(F'(1))

将第一个式子求导，(egin{aligned} F^{prime}(x)+G^{prime}(x) &=G^{prime}(x) cdot x+G(x) \ F^{prime}(1) &=G(1) end{aligned})，于是考虑求(G(1))

再将x=1带入第二个式子，注意(F(1)=1)，有

[egin{aligned} G(1) cdotleft(frac{1}{m} ight)^{L} &=sum_{i=1}^{L} a_{i} cdot F(1) cdotleft(frac{1}{m} ight)^{L-i} \ G(1) &=sum_{i=1}^{L} a_{i} cdot m^{i} end{aligned} ]

这个就能随便算了

code（hash有点丑，能过）

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
#include<ctime>
#define ll long long
#include<cstdlib>
#include<queue>
const int N = 1e5+1;
using namespace std;
inline int read(){
    char ch=getchar();int x=0;int pos=1;
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') pos=0;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';
    return pos?x:-x;
}
int a[N],s[N],n,t,m;
ll h[N],b[N];
const ll mod = 10000,sd=131,md1=998244353;
ll ksm(ll a,ll b){
	ll res=1;
	while(b){
		if(b&1) res=res*a%md1;
		a=a*a%md1;b>>=1;
	}
	return res;
}
ll ginv(ll a){
	return ksm(a,md1-2);
}
void get_hash(){
	memset(h,0,sizeof(h));ll res=1;
	for(int i=1;i<=n;i++){
		res=res*sd%md1;
		h[i]=(h[i-1]+a[i]*res%md1)%md1;
	}
}
ll ha(int l,int r){
	return ((h[r]-h[l-1]+md1)%md1)*b[l-1]%md1;
}
int main(){
	m=read();t=read();
	b[0]=1;ll si=ginv(sd);
	for(int i=1;i<N;i++) b[i]=b[i-1]*si%md1;
	while(t--){
		n=read();for(int i=1;i<=n;i++) a[i]=read(); 
		get_hash();
		ll ans=0,res=1;
		for(int i=1;i<=n;i++){
			res=res*m%mod;
			if(ha(1,i)==ha(n-i+1,n)) ans=(ans+res)%mod;
		}
		printf("%d%d%d%d
",ans/1000,ans/100%10,ans/10%10,ans%10);
	}
	return 0;
}

SDOI2017

本来写了十多行结果没保存，就不写了

跟上面差不多，对于每个字符串分别列出式子就能一共列出n+1个式子，高斯消元即可

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
#include<ctime>
#define ll long long
#include<cstdlib>
#include<queue>
const int N = 321;
using namespace std;
inline int read(){
    char ch=getchar();int x=0;int pos=1;
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') pos=0;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';
    return pos?x:-x;
}
int a[N],n,t,m;
char s[N];
ll h[N][N],b[N];
double f[N][N],_2[N],ans[N];
const ll mod =998244353,sd=131;
ll ksm(ll a,ll b){
	ll res=1;
	while(b){
		if(b&1) res=res*a%mod;
		a=a*a%mod;b>>=1;
	}
	return res;
}
ll ginv(ll a){
	return ksm(a,mod-2);
}
void get_hash(int t){
	ll res=1;
	for(int i=1;i<=m;i++){
		res=res*sd%mod;
		h[t][i]=(h[t][i-1]+((s[i]=='H')+1)*res)%mod;
	}
}
ll ha(int t,int l,int r){
	return ((h[t][r]-h[t][l-1]+mod)%mod)*b[l-1]%mod;
}
int main(){
	n=read(),m=read();
	b[0]=1;ll tx=ginv(sd);_2[0]=1;
	for(int i=1;i<=m;i++) b[i]=b[i-1]*tx%mod,_2[i]=_2[i-1]+_2[i-1];
	for(int i=1;i<=n;i++){
		scanf("%s",s+1);
		get_hash(i);
	}
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			for(int k=1;k<=m;k++){
				if(ha(i,1,k)==ha(j,m-k+1,m)) (f[i][j]+=_2[k]);
			}
		}
		f[i][n+1]=-1; 
	}
	for(int i=1;i<=n;i++) f[n+1][i]=1;f[n+1][n+2]=1;
	for(int i=1;i<=n+1;i++){
		for(int j=n+2;j>=i;j--) f[i][j]/=f[i][i];
		for(int j=i+1;j<=n+1;j++)
			for(int k=n+2;k>=i;k--)
				f[j][k]-=f[i][k]*f[j][i];
	}
	ans[n+1]=f[n+1][n+2];
	for(int i=n;i>=1;i--){
		ans[i]=f[i][n+2];
		for(int j=i+1;j<=n+1;j++) ans[i]-=ans[j]*f[i][j];
	}
	for(int i=1;i<=n;i++) printf("%.10f
",ans[i]);
	return 0;
}