[JZOJ100004]【NOI2017模拟.4.1】 Dice

题目大意

给你一个骰子，点数1~6，每个点数都有抛出的概率P_i。你要抛n轮，每轮中，如果抛出了跟上一轮一样的点数，要重新抛。
现在求两个值，n轮点数和ans的期望E[ans],方差的期望V[ans]。
方差定义：对于一个点数和为x的局面，它的方差为(x−E1)2
n<=100000，保证六个面概率加起来严格为1，精度误差10−4。

分析

首先看看E[ans]怎么算。
设Xi为第i轮点数期望，则E[ans]=E(∑Xi)=∑Xi=∑∑6j=1P(第i轮甩出j)∗j
考虑怎么求P(~)。
P(第i轮甩出j)=∑6k=1P(i−1轮抛出k)∗P(i−1轮抛出k后，i轮抛出j)
注意到：
1，P(第1轮甩出j)，就是Pj
2，最后面的那个P是可以提前算的，设g[x][y]为上次抛了x，这次抛y的概率。因为上次抛了x，这次不能抛x，设方程g[x][y]=Py+Px∗g[x][y]，可得g[x][y]=Py1−Px
那么我们可以递推出所有的P
设h[i][j]表示第i轮扔出j的概率，则
h[i][j]=∑6k=1h[i−1][k]∗g[k][j]
E[ans]=所有h[i][j]∗j
就算出了E[ans]。

接下来要要算V[ans]了。
搞一波事情：
V[ans]
=E[(ans−E[ans])2]
=E[ans2+E[ans]2−2∗ans∗E[ans]]
=E[ans2]+E[ans]2−2∗E[ans]∗E[ans]]（不相关的变量期望可以分开算）
=E[ans2]−E[ans]2
发现后面那一项是已经算过的了，只用算前一项。
E[ans2]=E[(∑Xi)2]
受上面求E[ans]的启发，我们试试继续递推。
设f[i][j]表示第i轮扔出j，前i轮所有点数的期望和，f2[i][j]表示期望和的平方。
显然f[i][j]=∑6k=1f[i−1][k]∗g[k][j]+j∗h[i][j]
而f2要稍微注意一下。先从期望式子看看吧。
设当前我们第i轮扔出了j。
f2[i][j]=E((∑i−1k−1Xk+j)2)
设sigma的东西为S。
=E((S+j)2)
=E(S2+j2+2∗S∗j)
=E(S2)+E(j2)+2∗E(S∗j)
看出来第一个是前面的f2[i-1]的东西，中间的是j*j*h[i][j]。后面的那个比较容易推错，我们单独看看。
是：E(S∗j)=E(S)∗E(j)吗？不对，因为S和j实际上是不独立的，因为第i−1轮的状态能够决定j的概率。
那怎么办呢？实际上，只要确定i-1的状态，E(s*j)里面的j就可以拿出来算，因为确定了i-1的状态之后，j在这个局面下是独立的了。也就是说，你把第i-1轮的甩出数当做常数了，那么一个随机变量跟常数有关，就等于独立了嘛。
当确定上一轮是k的时候，我们可以抽离j出来变成E(j)了，E(j)=j*g[k][j]，当然此时的E(S)也不是原来没有确定k的E(S)了。

从刚刚的分析我们可以得出要确定上一轮的值才能求出f2[i][j]，递推式也可以写出来了。
f2[i][j]=f2[i−1][k]∗g[k][j]+j2∗h[i][j]+2∗f[i−1][k]∗j∗g[k][j]
这样我们搞出了V[ans].
其实考场上的时候我被后面的哪个E(S*j)搞了好久，就是因为期望分析的时候乱拆括号，然后递推式就打错了。
从这里我们可以看出，递推式一定要以期望的式子为基础，然后期望的式子也可以由递推式启发来继续推下去。
哦写这篇东西的时候发现了一个东西：两个不是互相独立的变量，虽然不能拆开，但是你确定了一个变量之后，就相当于变成了常数，然后就可以和另一个变量独立开来，然后期望式就可以拆开了。分析如何确定一些变量使得期望式可以用递推式算出来，这应该是做题的一般套路吧。

其实问题没完，这东西卡精度。
为什么呢？我们最后的答案是(∑if2[n][i])−E[ans]2，这里是最后两个很大的数相减，想一下，E[ans]起码有n这么大，那么平方后有10位了，相减很难保证到10−4。怎么办？
回到之前方差的式子
V[ans]=E[(ans−E[ans])2]
其实我们不必这么快拆括号，仔细看看，发现E[ans]是一个常数了，式子相当于
V[ans]=E[(ans−const)2]
在一开始，我们就ans=∑iXi，那么
V[ans]=
E[((∑ni=1Xi)−const)2]
E[(∑ni=1(Xi−constn))2]
也就是说，对于每个点数i=1~6，原来每次投出后的贡献是i，现在只需要把贡献改成i−constn，再递推f2，最后就不用减去E[ans]2了。
这样就做完了。

代码

#include<cstdio>#include<algorithm>#include<cstring>#include<cmath>#include<bitset>using namespace std;typedef long long ll;typedef double db;#define fo(i,j,k) for(i=j;i<=k;i++)#define fd(i,j,k) for(i=j;i>=k;i--)const int N=100005;int i,j,k,n;double d[N],e[N],g[7][7],p[7],f[N][7],f2[N][7],h[N][7],incr,eans,vans,p2[N],pos,check[N];int main(){    freopen("dice.in","r",stdin);    freopen("dice.out","w",stdout);    fo(i,1,6) scanf("%lf",p+i);    fo(i,1,6)    {        incr=1/(1-p[i]);        fo(j,1,6)            if (i!=j)                g[i][j]=p[j]*incr;    }    scanf("%d",&n);    fo(i,1,6)    {        d[i]=db(i);        e[i]=d[i]*d[i];    }    fo(i,1,6)     {        f[1][i]=p[i]*d[i],f2[1][i]=p[i]*e[i],h[1][i]=p[i];    }    fo(i,1,n-1)    {        fo(j,1,6)            fo(k,1,6)                if (j!=k)                {                    pos=h[i][j]*g[j][k];                    f[i+1][k]=f[i+1][k]+g[j][k]*f[i][j]+d[k]*pos;                    h[i+1][k]=h[i+1][k]+pos;                }    }    fo(i,1,6)         eans=eans+f[n][i];    fo(i,1,6)    {        d[i]=db(i)-eans/db(n);        e[i]=d[i]*d[i];    }    fo(i,1,n)        fo(j,1,6) f[i][j]=0;    fo(i,1,6)         f[1][i]=p[i]*d[i],f2[1][i]=p[i]*e[i];    fo(i,1,n-1)    {        fo(j,1,6)            fo(k,1,6)                if (j!=k)                {                    pos=h[i][j]*g[j][k];                    f[i+1][k]=f[i+1][k]+g[j][k]*f[i][j]+d[k]*pos;                    f2[i+1][k]=f2[i+1][k]+e[k]*pos+(2*d[k]*f[i][j]+f2[i][j])*g[j][k];//why???                }    }    fo(i,1,6)         vans=vans+f2[n][i];    printf("%.10lf\n",eans);    printf("%.10lf\n",vans);}