我正在实施一个近似计数算法,其中我们:
使用log(log n)位维护计数器X.
将X初始化为0
当物品到达时,增加X的概率(1/2)由1 X
当流结束时,输出2 X - 1,使E [2 X ] = n + 1
我的实现如下:
import System.Random type Prob = Double type Tosses = Int -- * for sake of simplicity we assume 0 <= p <= 1 tos :: Prob -> StdGen -> (Bool,StdGen) tos p s = (q <= 100*p, s') where (q,s') = randomR (1,100) s toses :: Prob -> Tosses -> StdGen -> [(Bool,StdGen)] toses _ 0 _ = [] toses p n s = let t@(b,s') = tos p s in t : toses p (pred n) s' toses' :: Prob -> Tosses -> StdGen -> [Bool] toses' p n = fmap fst . toses p n morris :: StdGen -> [a] -> Int morris s xs = go s xs 0 where go _ [] n = n go s (_:xs) n = go s' xs n' where (h,s') = tos (0.5^n) s n' = if h then succ n else n main :: IO Int main = do s <- newStdGen return $ morris s [1..10000]
问题是,我的X总是不正确的任何|stream| > 2
,而且好像所有的StdGen
和|stream| > 1000
,X = 7
我在Matlab中测试了相同的算法,它在那里工作,所以我认为它也是
我的随机数生成器的问题,或
提高1/2到一个大的n in Double
请建议前进的道路?
这个问题实际上非常简单:randomR (1,100)
你排除了第一个百分比内的值,所以你在1/2的高倍率下有一个完全的截止值(所有这些都在那么小的间隔内).实际上是一般的事情:范围应该从零开始,而不是在一个†,除非有特定的原因.
但为什么甚至首先使用100的范围?我只是做到了
tos :: Prob -> StdGen -> (Bool,StdGen) tos p s = (q <= p, s') where (q,s') = randomR (0,1) s
† 我知道,Matlab在整个地方都弄错了.只是关于那种语言的许多可怕的事情之一.
与你的问题无关:正如chi所说,如果你使用一个合适的随机monad,这种代码看起来好多了,而不是手动传递StdGen
s.
import Data.Random import Data.Random.Source.Std type Prob = Double tos :: Prob -> RVar Bool tos p = do q <- uniform 0 1 return $ q <= p morris :: [a] -> RVar Int morris xs = go xs 0 where go [] n = return n go (_:xs) n = do h <- tos (0.5^n) go xs $ if h then succ n else n morrisTest :: Int -> IO Int morrisTest n = do runRVar (morris [1..n]) StdRandom