6

Haskellの仮想惑星のランダムな質量を生成しようとしています。バイモーダル分布（理想的には、2つの正規分布の重ね合わせ：1つは小さな惑星に対応し、もう1つは巨大ガスに対応）をサンプリングすることによって、これらの質量を生成したいと思います。一様分布を多数の分布に変換できる関数を提供する統計パッケージを見てきました。しかし、ディストリビューションの作成はサポートされていないようです。quantileDoubleDouble

この特定のケースは、事前にサンプリングするためにいずれかのディストリビューションを選択することでハッキングされる可能性がありますが、特に後でディストリビューション全体を微調整する必要がある場合があるため、単一のディストリビューションでそれを実行したいと思います。最終的には、正規分布を空の調査からの実際のデータに置き換える可能性があります。

棄却サンプリングを自分で実装することを検討しています。これは、任意の分布をかなり簡単に処理できますが、かなり非効率的であり、ソリューションがライブラリとしてすでに存在する場合は、実装するのは確かに良い考えではありません。

構成された、または明示的に指定された分布からのサンプリングをサポートするHaskellライブラリはありますか？または、棄却サンプリングの既存のHaskell実装ですか？あるいは、2つの正規分布の合計のCDFの逆数の明示的な式はありますか？

4

2 に答える 2

6

分布の単純な混合の場合、最初に言及した「ハック」を介して効率的なサンプラーを取得できます。

この特定のケースは、事前にサンプリングするためにいずれかのディストリビューションを選択することでハッキングされる可能性がありますが、特に後でディストリビューション全体を微調整する必要がある場合があるため、単一のディストリビューションでそれを実行したいと思います。

これは実際にはギブスサンプリングの場合であり、統計で非常に一般的です。それは非常に柔軟性があり、使用している混合物の数を知っていれば、おそらく打ち負かすのは難しいでしょう。アンサンブル全体から1つの個別の分布を選択してサンプリングし、その条件付き分布からサンプリングします。すすぎ、繰り返します。

これは、ガウス混合ギブスサンプラー用の単純で最適化されていない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.

実行例を次に示します。2つのガウス分布の1次元混合から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の良い紹介です。

于 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でサンプリングするだけです。IO

于 2012-05-31T08:01:26.463 に答える