6

我正在尝试为 Haskell 中的假设行星生成随机质量。我想通过对双峰分布（理想情况下是两个正态分布的叠加：一个对应于小行星，一个对应于气态巨行星）进行采样来产生这些质量。我查看了统计包，它提供了quantile可以将一个均匀分布Double变成一个Double对多个分布的函数。但似乎没有任何支持编写发行版。

这种特殊情况可以通过选择一个分布或另一个分布来预先采样，但我想用一个分布来做，特别是因为我以后可能需要调整整体分布。最终，我可能会用来自天空调查的真实数据代替正态分布。

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

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

4

2 回答 2

6

在简单混合分布的情况下，您可以通过您首先提到的“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 例程。 这是对 Monte Carlo 和 MCMC 的一个很好的介绍。

于 2012-05-31T08:27:00.300 回答
3

嗯。我熟悉的最好方法是调整MonadRandom 包以获得“概率单子”，从http://en.wikipedia.org/wiki/Normal_distribution#Generating_values_from_normal_distribution借用一些工具：

getRandomStrictlyBetween :: (Ord a, Random a, RandomGen m) => 
    (a, a) -> a
getRandomStrictlyBetween (lo, hi) = do
  x <- getRandomR (lo, hi)
  -- x is uniformly randomly chosen from the *closed* interval
  if lo < x && x < hi then return x else getRandomStrictlyBetween (lo, hi)

normalValue :: MonadRandom m => m Double
normalValue = do
  u <- getRandomStrictlyBetween (0, 1)
  v <- getRandomStrictlyBetween (0, 2 * pi)
  return (sqrt (-2 * log u) * cos v) -- according to Wikipedia

然后您可以得出或多或少的任意分布；例如，要获得y具有概率p和z概率的随机变量的分布(1 - p)，您只需编写

do alpha <- getRandom -- double chosen from [0, 1)
   if alpha < p then y else z

其中双峰分布似乎是一种特殊情况。要从这些分布evalRandIO distribution中采样，只需在IOmonad 中采样。

于 2012-05-31T08:01:26.463 回答