binomial-random ( n p -- x )

Vocabulary

random

Inputs

n	an object
p	an object

Outputs

x

an object

Definition

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

IN: random

:: binomial-random ( n p -- x )
    n assert-non-negative drop p 0.0 1.0 between? [
        {
            { [ p 0.0 = ] [ 0 ] }
            { [ p 1.0 = ] [ n ] }
            { [ n 1 = ] [ random-unit p < 1 0 ? ] }
            { [ p 0.5 > ] [ n dup 1.0 p - binomial-random - ] }
            {
                [ n p * 10.0 < ]
                [
                    1.0 p - log :> c 0 0 random-generator get [
                        (random-unit) log c /i 1 + + dup n <=
                        dup ~quotation~ when
                    ] curry loop drop
                ]
            }
            [
                n p * 10.0 >= p 0.5 <= and t assert=
                random-generator get :> rnd 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?
                    [ ~array~ 0&& ~quotation~ ~quotation~ if ]
                    [ t ] if
                ] loop
            ]
        } cond
    ] [ "p must be in the range 0.0 <= p <= 1.0" throw ] if ;