一个小问题的求解

问题

一根 1m 长的玻璃棒，摔倒地上断成 3 段，最短一段的平均值是多少？

假设玻璃棒一定会摔成三段，且玻璃棒质地均匀，为理想状态。

物理的视角

玻璃棒摔成三段，其物理过程是什么样的？

有一位网友分析了倾斜碰撞的应力分布，提出了断裂位置的可能性分布，并把一次断成三阶分解为同时发生的两次物理上相似的断裂。

其实玻璃棒的断裂更可能是地面的冲击力导致的，而不是倾斜碰撞。这样的话，断裂位置的可能性分布就是均匀的，而不是倾斜的。这里的均匀意思是地面可能与玻璃棒接触（形成单点接触）的位置可能是导致断裂的原因，而不一定是一段倾斜碰撞带来的应力集中。

数学的视角

从数学的角度，忽略物理的过程，可以有不同的假设。

单纯考虑断裂为三段，因为题目中已经假设了必然断裂为三段，类似于剪刀剪断绳索的情况。那么这三段的长度分别为 $x, y, z$ ，且有 $x + y + z = 1$ 。

很明显，三者的长度无法倍当作独立事件来考虑，最终只能定义为两个独立事件，也就是断裂的点（两个）考虑为独立事件。

两个断裂段的位置定义为 $r_1, r_2$ ，则断裂的三段长度为：

$\begin{aligned} x &= \min_{i=1,2}r_i \\ y &= \left|r_1-r_2\right| \\ z &= 1 - \max_{i=1,2}r_i\\ \end{aligned}$

其中 $r_i$ 的概率密度函数为：

$f(r_i) = \begin{cases} 1, & 0 \le r_i \le 1 \\ 0, & \text{otherwise} \end{cases}$

蒙特卡洛模拟

蒙特卡洛模拟的思路是，随机生成 $r_1, r_2$ ，计算 $x, y, z$ ，然后计算 $min \{x,y,z\}$ 的平均值。

open System
let r = Random(42)

let break2 () =
    let a = r.NextDouble()
    let b = r.NextDouble()     
    [min a b; abs (b-a); 1. - max a b]


let breakN n =
    seq {
        for _ in (bigint 1) .. n do
            yield break2 ()
    }

let minAverage n =
    breakN n
    |> Seq.map List.min
    |> Seq.average

for i in 1..7 do
    let n = 10. ** float i |> bigint
    printfn $"%4A{i} %20A{n} %12f{minAverage n}"

得到的结果：

1 10     0.129573
2 100     0.096320
3 1000     0.105467
4 10000     0.110859
5 100000     0.111424
6 1000000     0.111146
7 10000000     0.111107

接近某网友分析的 $1/9$ 。

更广义的伪随机采样算法

上面这个问题，属于一类称为伪随机采样的计算方法。伪随机采样的思路是，随机生成一些数据，然后计算这些数据的某个统计量，这个统计量就是我们想要的结果。典型的问题就是上面的概率问题，也就是估计一阶统计量（平均值）。还有一类就是估计二阶统计量，比如方差。

求解问题的关键有两个部分：

自变量的统计分布和采样生成；
统计量的计算。

对于第一个问题，我们可以使用一些已知的分布，比如均匀分布、正态分布等等。也可以使用一些已知的采样方法，比如拒绝采样、重要性采样等等。

第二个问题一般包括了对物理过程的数学模型的计算，以及对采样数据的统计量计算。

除了对统计量的计算，Monte carlo方法还可以用于确定问题的转化计算，就是可以把一个物理量视为统计量或者统计量的函数。比如前面的 $\pi$ 的计算问题，转化为面积的比值，也就是概率的比值。

一个典型的应用就是求解积分。上面的圆周率计算，实际上也就是求圆的面积，也就是求解积分。多重积分的Monte Carlo计算也是一个很常用的方法。

多重积分的计算

问题定义和算法

$\int\int\int\cdots\int f(\mathbf{x})d\mathbf{x}$

形如上面的积分，当积分重数很大时，通过分步积分的方法很难求解。这时候可以使用Monte Carlo方法来求解。

求解过程也很简单，就是随机生成一些点，然后计算这些点的函数值的平均值，乘以积分区域的体积。

$\mathbf{x} \in \prod_{i=1}^n [a_i, b_i]$

其中， $a_i, b_i$ 伪定义域的上下界。

$\hat{I}_N = \int\int\int\cdots\int f(\mathbf{x})d\mathbf{x} = \prod_{i=1}^n (b_i - a_i) \frac{1}{N}\sum_{j=1}^N f(\mathbf{x}_j)$

这里 $N$ 就是采样数目，随着 $N$ 的增加，计算结果会越来越趋近真值。

$\hat{I}_{N} |_{N\to \infty}$

F#代码示例

下面就用F#来实现一个多重积分的计算程序。

// Monte Carlo Integration
open System

// random float in [lb, ub]
let r = Random(42)
let rf (lb, ub) = r.NextDouble() * (ub - lb) + lb

// random sample in bounds
let sample (bounds: (float*float) array) =
    bounds
    |> Array.map rf

// volumn of bounds
let volumn (bounds: (float*float) array) = 
    bounds
    |> Array.map (fun (lb, ub) -> Math.Abs(ub - lb))
    |> Array.sum

// Monte Carlo Integration
let mc f (bounds: (float*float) array) N =
    [|1..N|]
    |> Array.map (fun _ -> sample bounds)  // sample N points
    |> Array.map f // f(x)
    |> Array.average // 1/N * sum(f(x))
    |> fun x -> x * (volumn bounds)

上面就是多重积分的全部代码，要求目标函数的参数是一个float array，返回值是一个float。积分上下限是一个float*float的元组数组。

下面是测试算法收敛性的代码，计算的采样数目是 $2^N$ ， $N$ 从10到25。

// Test
let f (x: float array) =
    x
    |> Array.map (fun x -> x * x)
    |> Array.sum
let b = Array.init 10 (fun _ -> (-1.0, 1.0))

let Ns = [|10..25|] |> Array.map (fun i -> int(Math.Pow(2, i)))
let Ifs = Ns |> Array.map (mc f b)

利用ScottPlot绘制图像：


// Plot

#r "nuget: ScottPlot, 4.1.0"

open ScottPlot
let plt = Plot(600, 400)
plt.AddScatter(Ns |> Array.map float, Ifs)
plt.Title("Monte Carlo Integration convergence")
plt.XLabel("N")
plt.YLabel("I")
plt.SaveFig("img/mc.png")

在这里插入图片描述