Project Euler 697

题目

Randomly Decaying Sequence

Given a fixed real number $c$, define a random sequence $(Xn){n\ge 0}$ by the following random process:

• $X_0 = c$ (with probability $1$).
• For $n>0$, $Xn = U_n X{n-1}$ where $Un$ is a real number chosen at random between zero and one, uniformly, and independently of all previous choices $(U_m){m<n}$.

If we desire there to be precisely a $25\%$ probability that $X{100}<1$, then this can be arranged by fixing $c$ such that $\log{10} c \approx 46.27$.

Suppose now that $c$ is set to a different value, so that there is precisely a $25\%$ probability that $X_{10\,000\,000}<1$.

Find $\log_{10} c$ and give your answer rounded to two places after the decimal point.

解决方案

令$n=10^7$。原问题是求$c$使得$P {c\cdot\prod_{i=1}^n U_i<1}=0.25$。

令$Vi=-\ln U_i,v=\ln c$。那么问题就变成求$v$使得$P{v-\sum{i=1}^n V_i< 0}=0.25$。

由于$U_i\sim U(0,1)$，因此$V_i$是服从参数为$1$的指数分布，即$V_i\sim \text{Exp}(1)$。

那么，随机变量$Y=\sum_{i=1}^n V_i$服从参数为$(n,1)$的伽马分布，即有$Y\sim \Gamma(n,1)$。

由于目前是求$v$使得$P{Y>v}=0.25$，也就是$\Gamma(n,1)$的$0.75$分位点。根据伽马分布和伽马函数的定义，随机变量$Y$的概率密度函数为

$$f_Y(x)=x^{n-1}\cdot e^{-x}$$

那么$Y$的分布函数恰好为不完全伽马函数$\Gamma(n,x)$。我们使用scipy.special库中的方法gammainc(n,x)来计算$\Gamma(n,x)$的值。为了确定$v$使得$\Gamma(n,v)=0.75$，使用二分法进行完成。

最终，$\dfrac{v}{\ln 10}$为答案。

代码

from scipy.special import gammainc
from math import log

N = 10 ** 7
R = 0.25
R = 1 - R

l, r = 0, 10 ** 20
for _ in range(100):
    m = (l + r) * 0.5
    v = gammainc(N, m)
    if v < R:
        l = m
    else:
        r = m

ans = l / log(10)
print("{:.2f}".format(ans))