haskell - 如何从 Haskell 中的复杂或复合分布中采样？

Question

我正在尝试为 Haskell 中假设的行星生成随机质量。我想通过采样双峰分布(理想情况下是两种正态分布的叠加:一种对应于小行星，一种对应于气态巨行星)来产生这些质量。我看过statistics package ，它提供了quantile函数，可以将均匀分布的Double转换为多个分布上的Double。但似乎不支持编写发行版。

这种特殊情况可以通过预先选择一个发行版或另一个发行版进行采样来解决，但我想使用单个发行版来完成此操作，特别是因为我稍后可能需要调整整体发行版。最终我可能会用天空调查的真实数据替换正态分布。

我正在考虑实现 rejection sampling我自己，它可以相当简单地处理任意分布，但似乎效率相当低，如果解决方案已经作为库存在，那么实现它肯定不是一个好主意。

是否有 Haskell 库支持从组合或明确指定的分布中采样？或者现有的 Haskell 拒绝抽样实现？或者，是否有两个正态分布之和的 CDF 的倒数的明确公式？

Answer 1

在简单混合发行版的情况下，您可以通过您首先提到的“黑客”获得高效的采样器:

This particular case could be hacked around by picking one distribution or the other to sample before-hand, but I'd like to do it with a single distribution, especially since I might need to tweak the overall distribution later.

这实际上是统计中非常普遍的吉布斯抽样的情况。它非常灵活，如果您知道所使用的混合物的数量，它可能很难被击败。从整个集合中选择一个单独的分布进行采样，然后从该条件分布中进行采样。冲洗并重复。

这是一个简单的、未优化的 Haskell 实现，用于混合高斯吉布斯采样器。这是非常基本的，但你明白了:

import System.Random
import Control.Monad.State

type ModeList = [(Double, Double)]                 -- A list of mean/stdev pairs, for each mode.

-- Generate a Gaussian (0, 1) variate.
boxMuller :: StdGen -> (Double, StdGen)
boxMuller gen = (sqrt (-2 * log u1) * cos (2 * pi * u2), gen'')
    where (u1, gen')  = randomR (0, 1) gen 
          (u2, gen'') = randomR (0, 1) gen'

sampler :: ModeList -> State StdGen Double
sampler modeInfo = do
    gen <- get
    let n           = length modeInfo
        (z0, g0)    = boxMuller gen
        (c,  g1)    = randomR (0, n - 1) g0        -- Sample from the components.
        (cmu, csig) = modeInfo !! c                
    put g1
    return $ cmu + csig * z0                       -- Sample from the conditional distribution.

下面是一个运行示例:从两个高斯的一维混合中采样 100 次。众数为 x = -3 和 x = 2.5，每个混合分量都有自己单独的方差。您可以在此处添加任意数量的模式。

main = do
let gen      = mkStdGen 42
    modeInfo = [(2.5, 1.0), (-3, 1.5)]
    samples     = (`evalState` gen) . replicateM 100 $ sampler modeInfo
print samples

这是这 100 个样本的平滑密度图(使用 R 和 ggplot2):

更通用的算法是拒绝或重要性采样器，在更复杂的分布的情况下，您可能需要手动滚动适当的 MCMC 例程。 Here是对蒙特卡洛和 MCMC 的很好的介绍。