binomial-random* ( n p rnd -- elt )

Vocabulary

random

Inputs

n	an object
p	an object
rnd	an object

Outputs

elt

an object

Definition

USING: combinators combinators.short-circuit kernel math
math.functions math.order ;

IN: random

:: binomial-random* ( n p rnd -- elt )
    n assert-non-negative drop {
        {
            [ p 0.0 1.0 between? not ]
            [ "p must be in the range 0.0 <= p <= 1.0" throw ]
        }
        { [ p 0.0 = ] [ 0 ] }
        { [ p 1.0 = ] [ n ] }
        { [ n 1 = ] [ rnd random-unit* p < 1 0 ? ] }
        { [ p 0.5 > ] [ n dup 1.0 p - rnd binomial-random* - ] }
        {
            [ n p * 10.0 < ]
            [
                1.0 p - log :> c 0 0 [
                    rnd random-unit* log c /i 1 + + dup n <= dup
                    [ ~quotation~ 2dip ] when
                ] loop drop
            ]
        }
        [
            n p * 10.0 >= p 0.5 <= and t assert= n p * 1.0 p - *
            sqrt :> spq 1.15 2.53 spq * +
            :> b -0.0873 0.0248 b * + 0.01 p * + :> a n p *
            0.5 + :> c 0.92 4.2 b / - :> vr 2.83 5.1 b / + spq *
            :> alpha p 1.0 p - / log :> lpq n 1 + p *
            floor :> m m 1 + lgamma n m - 1 + lgamma + :> h f [
                rnd random-unit* 0.5 -
                :> u rnd random-unit* :> v 0.5 u abs -
                :> us drop 2.0 a us / * b + u * c +
                floor >integer dup :> k k 0 n between? [
                    { ~quotation~ ~quotation~ } 0&&
                    [ f ] [
                        alpha a us sq / b + / v * log h k 1 +
                        lgamma - n k - 1 + lgamma - k m - lpq *
                        + >
                    ] if
                ] [ t ] if
            ] loop
        ]
    } cond ;