2012-01-09 81 views
5

我需要生成一個隨機整數的無限流,數字在範圍[1..n]。然而,預先給出每個數字p_i的概率,因此分佈不均勻。生成給定概率的隨機整數

在Haskell中是否有庫函數?

回答

3

只是對dflemstr的答案擴大,您可以使用Control.Monad.Random這樣創造的加權值的無窮列表:

import Control.Monad.Random 
import System.Random 

weightedList :: RandomGen g => g -> [(a, Rational)] -> [a] 
weightedList gen weights = evalRand m gen 
    where m = sequence . repeat . fromList $ weights 

而且使用這樣的:

> let values = weightedList (mkStdGen 123) [(1, 2), (2, 5), (3, 10)] 
> take 20 values 
[2,1,3,2,1,2,2,3,3,3,3,3,3,2,3,3,2,2,2,3] 

這並不需要IO monad,但您需要提供用於流的隨機數生成器。

5

Control.Monad.Random提供在fromList:: MonadRandom m => [(a, Rational)] -> m a

形式此功能可以在IO單子與使用它:

import Control.Monad.Random 
-- ... 
someNums <- evalRandIO . sequence . repeat . fromList $ [(1, 0.3), (2, 0.2), (3, 0.5)] 
print $ take 200 someNums 

有運行Rand單子的其他方式,你可以在包看。權重不必加起來1

編輯:Rand比我想象的明顯懶惰,所以replicateM n可以通過sequence . repeat取代,如@shang建議。

11

正如人們指出Control.Monad.Random中有一個函數,但它的複雜度相當差。下面是我今天上午寫的一些巧合的代碼。它使用美麗的Alias算法。

module Data.Random.Distribution.NonUniform(randomsDist) where 
import Data.Array 
import Data.List 
import System.Random 

genTable :: (Num a, Ord a) => [a] -> (Array Int a, Array Int Int) 
genTable ps = 
    let n = length ps 
     n' = fromIntegral n 
     (small, large) = partition ((< 1) . snd) $ zip [0..] $ map (n' *) ps 
     loop ((l, pl):ls) ((g, pg):gs) probs aliases = 
      let prob = (l,pl) 
       alias = (l,g) 
       pg' = (pg + pl) - 1 
       gpg = (g, pg') 
      in if pg' < 1 then loop (gpg:ls) gs (prob:probs) (alias:aliases) 
          else loop ls (gpg:gs) (prob:probs) (alias:aliases) 
     loop ls gs probs aliases = loop' (ls ++ gs) probs aliases 
     loop' [] probs aliases = (array (0,n-1) probs, array (0,n-1) aliases) 
     loop' ((g,_):gs) probs aliases = loop' gs ((g,1):probs) ((g, -1):aliases) 
    in loop small large [] [] 

-- | Generate an infinite list of random values with the given distribution. 
-- The probabilities are scaled so they do not have to add up to one. 
-- 
-- Uses Vose's alias method for generating the values. 
-- For /n/ values this has O(/n/) setup complexity and O(1) complexity for each 
-- generated item. 
randomsDist :: (RandomGen g, Random r, Fractional r, Ord r) 
      => g       -- | random number generator 
      -> [(a, r)]     -- | list of values with the probabilities 
      -> [a] 
randomsDist g xps = 
    let (xs, ps) = unzip xps 
     n = length xps 
     axs = listArray (0, n-1) xs 
     s = sum ps 
     (probs, aliases) = genTable $ map (/ s) ps 
     (g', g'') = split g 
     is = randomRs (0, n-1) g' 
     rs = randoms g'' 
     ks = zipWith (\ i r -> if r <= probs!i then i else aliases!i) is rs 
    in map (axs!) ks 
+0

我只是想這個「總時間55.59s」對執行此:http://idontgetoutmuch.wordpress.com/2014/08/26/haskell-vectors-and-sampling-from-a-categorical分佈/「總時間11.09s」在兩種情況下采樣2 * 10^7個樣本。也許這不是一個公平的比較,因爲使用System.Random和其他System.Random.MWC。 – idontgetoutmuch 2014-08-28 06:57:35

+0

是的,我會假設生成隨機數將主宰我的代碼。它還需要專業化,這可能會在-O2時自動發生。 – augustss 2014-08-28 07:46:13

+0

使用不同的隨機數發生器,我得到了「總時間20.31s」,但仍然不如以前好。我還沒有嘗試過專業化。內存使用情況也不好。我希望兩個表中的每個條目都有4 + 8個字節,因此應該是2 * 12 * 10^7個字節,因此小於1G。我看到大約5G。我可能很天真。我還沒有讀完Devroye和Vose。誰會想到你可以用隨機數字獲得如此多的樂趣。 – idontgetoutmuch 2014-08-29 08:33:40