NB. statistical distributions

NB.*betarand        v random numbers in a beta distribution
NB.
NB.*binomialdist    v discrete values for binomial distribution
NB.*binomialprob    v probability of success in binomial distribution
NB.*binomialrand    v random numbers 0 and 1 in binomial distribution
NB.
NB.*cauchyrand      v random numbers in a cauchy distribution
NB.
NB.*discreterand    v random numbers in a discrete distribution
NB.
NB.*exponentialrand v random numbers in an exponential distribution
NB.
NB.*gammarand       v random numbers in a gamma distribution
NB.
NB.*normalprob      v probability of success in normal distribution
NB.*normalrand      v random numbers in a standard normal distribution,
NB.
NB.*poissondist     v discrete values for poisson distribution
NB.*poissonprob     v probability of success in poisson distribution
NB.*poissonrand     v random numbers in a poisson distribution

rand01=: ?@$ 0:

NB. =========================================================
NB. betarand
NB. random numbers in a beta distribution
NB. y. has 3 elements p,q,n
NB.   0 = power parameter
NB.   1 = power parameter
NB.   2 = number of trials
betarand=: 3 : 0
'p q n'=. y.
if. (1>p) *. 1>q do.
  b=. n#1
  r=. n#0
  whilst. 1 e. b do.
    m=. +/b
    x=. (rand01 m)^%p
    y=. x+(rand01 m)^%q    
    t=. 1>:y
    z=. (t#x)%t#y
    i=. t#b#i.#b
    b=. 0 i } b
    r=. (z+i{r) i } r
  end.
else.
  s%(gammarand q,n)+s=. gammarand p,n
end.
)

NB. =========================================================
NB. binomial distribution
NB. y. has 2 elements p,n:
NB.  0  =  probability of success in one trial
NB.  1  =  number of trials
NB.
NB. E=n*p, var=n*p*-.p
NB.
NB. e.g.  binomialdist 0.25 5  ...distribution of results from 0 to 5
binomialdist=: 3 : 0
'p n'=. y.
if. b=. p<:0.5 do. p=. -. p end.
q=. -.p
r=. }.i.>:n
r=. (q^n)**/\1,(p%q)*(>:n-r)%r
if. b do. |. r end.
)

NB. =========================================================
NB. probability of successes in binomial distribution
NB. y. has 3 or 4 elements:
NB.   0  =  probability of success in one trial
NB.   1  =  number of trials
NB.   2  =  minimum number of successes
NB.  {3} =  maximum number of successes
NB.
NB. e.g.  binomialprob 0.25 100 30 40    ...probability that there are
NB. between 30 and 40 successes in 100 trials, where probability = 0.25
binomialprob=: 3 : 0
'p n s'=. 3{.y.
m=. >:3{y.,n
+/ s}.m{.binomialdist p,n
)

NB. =========================================================
NB. random numbers 0 and 1 in a binomial distribution
NB.
NB. y. has 2 elements:
NB.  0  =  probability of success in one trial
NB.  1  =  number of trials
NB.
binomialrand=: 3 : 0
'p n'=. y.
r=. 2147483647
s=. <:p*r
s>?n#<:r
)

NB. =========================================================
NB. random numbers in a cauchy distribution
NB. with F(x)=0.5+(arctan x)%o.1
NB.
NB. y.  =  number of trials
cauchyrand=: 3 : 0
3 o. o. _0.5+rand01 y.
)

NB. =========================================================
NB. random numbers in a discrete distribution
NB.
NB. y. has two elements
NB.   0 =  2-row matrix: 0 =  discrete values
NB.                      1 =  corresponding probabilities
NB.   1 =  number of trials
discreterand=: 3 : 0
'm n'=. y.
'v p'=. m
f=. 0,+/\p%+/p
l=. #f
r=. /:f,rand01 n
s=. +/\l>r
r=. s r } r
v {~ <:l}.r
)

NB. =========================================================
NB. random numbers in an exponential distribution
NB. with mean=1. F(x)=1-^-x
NB.
NB. y.  =  number of trials
exponentialrand=: 3 : 0
-^.rand01 y.
)

NB. =========================================================
NB. gammarand
NB. random numbers in a gamma distribution
NB. y. has 2 elements p,n
NB.   0 = power parameter
NB.   1 = number of trials
NB.
NB. if p=1 this is the exponential distribution
gammarand=: 3 : 0
'p n'=. y.
r=. n#0
k=. p-i=. <.p
if. k do.
  r=. betarand k,(-.k),n
  r=. r * -^.rand01 n
end.
if. i do.
  r-^.*/rand01 i,n
end.
)

NB. =========================================================
NB. probabilities of normal distribution
NB.
NB. y. has 3 or 4 elements:
NB.   0  =  mean of distribution
NB.   1  =  standard deviation
NB.   2  =  minimum result
NB.  {3} =  maximum result
NB.
NB. e.g.  normalprob 0 1 2    ...probability that result will be
NB. greater than 2, where mean = 0, standard deviation 1
NB.
NB. Result is rounded to 7 decimal places
NB. Algorithm has accuracy ~ 7.5*1e_8
normalprob=: 3 : 0
'm s'=. 2{.y.
x=. s %~ m -~ 2}.y.
sgn=. *x
x=. |x
b=. 0.31938153 _0.356563782 1.781477937 _1.821255978 1.330274429
p=. 0.2316419
r=. ((2p1)^_0.5)*^_0.5**:x
r=. r*b+/ .**/\5#(1,#x)$%>:p*x
r=. -/|(sgn=_1)-r
(<.0.5+r*c)%c=. 1e7
)

NB. =========================================================
NB. random numbers in a standard normal distribution,
NB. with  mean=0, standard deviation=1.
NB. y.  = number of trials.

NB. NB. Polar method (J502 and earlier)
NB. normalrand=: 3 : 0
NB. r=. i.0
NB. while. y. ~: #r do.
NB.   j=. >.0.65*y.-#r
NB.   t=. (2,j)$<:+: rand01 +:j
NB.   b=. 1>s=. +/*:t
NB.   s=. 0.5^~(_2*^.s)%s=. b#s
NB.   r=. r,,(b#"1 t)*(2,#s)$s
NB.   r=. (y.<.#r){.r
NB. end.
NB. )

normalrand=: 3 : 0
(2 o. +: o. rand01 y.) * %: - +: ^. rand01 y.
)

NB. =========================================================
NB. discrete values for poisson distribution
NB.
NB. y. has 2 elements:
NB.   0 = mean of distribution
NB.   1 = maximum value to show
NB.
NB. e.g.  poissondist 2 10
NB.       = list of probabilities of values from 0 to 10
NB.         in poisson distribution of mean 2
poissondist=: 3 : 0
'm n'=. y.
(^-m)**/\1,m%}.i.>:n
)

NB. =========================================================
NB. probability of successes in poisson distribution
NB.
NB. y. has 2 or 3 elements:
NB.  0   =  mean of distribution
NB.  1   =  minimum number of successes
NB.  {2} =  maximum number of successes
NB.
NB. e.g.
NB.   poissonprob 5 4 7    ...probability that result will be
NB.         between 4 and 7 inclusive, where mean result is 5.
NB.   poissonprob 5 6 6    ...probability that result will be
NB.         exactly 6, where mean result is 5.
poissonprob=: 3 : 0
m=. {.y.
s=. }.y.
b=. 1=#s
'x y'=. _2{.0,s-b
if. 0>y do. 1 return. end.
|b-(^-m)*+/x}.*/\1,m%}.i.>:y
)

NB. =========================================================
NB. random numbers in a poisson distribution
NB.
NB. y. has 2 elements:
NB.   0  = mean of distribution (=variance)
NB.   1  = number of trials
poissonrand=: 3 : 0
'm n'=. y.
roll=. -@^.@rand01
r=. b=. m >: t=. roll n
i=. i.n
while. #i=. b#i do.
  b=. m >: t=. (b#t) + roll #i
  r=. (b + i{r) i } r
end.
r
)
