如何从Haskell中的复杂或复合分布中进行采样？

Question

我正在尝试为Haskell中的假想行星生成随机质量.我想通过采样双模态分布(理想情况下是两个正态分布的叠加:一个对应于小行星,一个对应于气体巨行)来产生这些质量.我查看了统计软件包,它提供了quantile函数,可以将统一分布Double转换Double为多个分布.但似乎没有任何支持撰写发行版.

这个特殊情况可以通过选择一个分布或另一个分类来进行预测,但是我想用一个分发来做,特别是因为我可能需要稍后调整整体分布.最终,我可能会用天空测量中的真实数据替换正态分布.

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

是否有一个Haskell库支持从组合或显式指定的分发中进行采样？或者现有的Haskell实现拒绝采样？或者,是否存在两个正态分布之和的CDF逆的显式公式？

Answer 1

在简单混合分布的情况下,您可以通过首次提到的"hack"获得高效的采样器:

这个特殊情况可以通过选择一个分布或另一个分类来进行预测,但是我想用一个分发来做,特别是因为我可能需要稍后调整整体分布.

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

这是一个简单的,未经优化的Haskell实现,用于高斯混合的Gibbs采样器.这是非常基本的,但你明白了:

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例程. 以下是蒙特卡洛和MCMC的精彩介绍.