F#奇妙游（22）：Monte Carlo方法的F#实现

一个小问题的求解

问题

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

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

物理的视角

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

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

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

数学的视角

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

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

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

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

$$\begin{array}{rl}x& =\underset{i=1,2}{\min }{r}_i\\ y& =\left|r_1-r_2\right|\\ z& =1-\underset{i=1,2}{\max }{r}_i\end{array}$$

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

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

蒙特卡洛模拟

蒙特卡洛模拟的思路是，随机生成 $r_1,r_2$ ，计算 $x,y,z$ ，然后计算 min ⁡ { x , y , z } \min \{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$ 。

更广义的伪随机采样算法

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

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

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

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

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

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

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

多重积分的计算

问题定义和算法

$$I=\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$ 的增加，计算结果会越来越趋近真值。

$$I={\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.reduce (*) // fixed 2023/08/24 13:18

// 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")

可以很好的看到10重积分依然能做到很高的收敛精度。

这个算法也算是求解多重积分的一个很高效的方法了，相对于分步积分，这个方法的收敛速度更快，而且可以很容易的并行化。

结论

1. F#编写数值算法代码非常方便，通过对数组进行函数式的操作，避免采用C和C++的for循环，算法的逻辑更加清晰，代码也更加简洁。
2. Monte Carlo方法可以用于求解很多问题，比如求解积分、求解概率、求解物理量等等。其算法简单粗暴容易实现。

修正

• 2023/08/24 修正文章的函数volume，更正最后的收敛图。