The randomized algorithm does not behave as expected

Question

I use an approximate calculation algorithm, where we:

Maintain X counter using the log (log n) bit

• Initialize X to 0

• When an item arrives, increase X by 1 with the probability of (& frac12;) ^X

• When the stream is finished, print 2 ^X & minus; 1, so that E [2 ^X ] = n + 1

My implementation is as follows:

import System.Random type Prob = Double type Tosses = Int -- * for sake of simplicity we assume 0 <= p <= 1 tos :: Prob -> StdGen -> (Bool,StdGen) tos ps = (q <= 100*p, s') where (q,s') = randomR (1,100) s toses :: Prob -> Tosses -> StdGen -> [(Bool,StdGen)] toses _ 0 _ = [] toses pns = let t@ (b,s') = tos ps in t : toses p (pred n) s' toses' :: Prob -> Tosses -> StdGen -> [Bool] toses' pn = fmap fst . toses pn 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) sn' = if h then succ n else n main :: IO Int main = do s <- newStdGen return $ morris s [1..10000] 

The problem is that my X is always incorrect for any |stream| > 2 |stream| > 2 , and it seems that for all StdGen and |stream| > 1000 |stream| > 1000 , X = 7

I tested the same algorithm in Matlab and it works there, so I assume that it

• problem with my random number generator or

• rise 1/2 to large n in Double

Please suggest a way forward?

Answer 1

The problem is actually very simple: with randomR (1,100) you exclude values ​​within the first percent, so you have full cutoff at high 1/2 degrees (they all lie in this small interval). In fact, the general: ranges should start from scratch , and not on one ^{& dagger;} unless there is a specific reason.

But why use the range of 100 at all? I would just do it

 tos :: Prob -> StdGen -> (Bool,StdGen) tos ps = (q <= p, s') where (q,s') = randomR (0,1) s 

---

^{& dagger;} _{I know that Matlab is mistaken everywhere. Just one of the many terrible things about this language.}

---

Not related to your problem: as Qi noted, such code looks much nicer if you use a suitable random monad instead of manually StdGen around 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