The following file contains a logfact function that uses loggamma
It's pretty old so might need revamping!
--------------------------------------------------------------------------------------------
NB. J.E.H.Shaw probability distributions library
NB.---------------------------------------------
NB. Created 5-Feb-1998
NB. Last modified 11-Aug-1999
NB.
NB. <dist>c = cdf of <dist> at y. [optional pars x.]
NB. <dist>f = pdf/pmf of <dist> at y. [optional pars x.]
NB. <dist>q = inverse cdf (quantiles) of <dist> at y. [optional pars x.]
NB. <dist>r = y. pseudorandom numbers from <dist> [optional pars x.]
NB. <dist>s = survivor function (1-cdf) of <dist> at y. [optional pars x.]
NB.
NB. ====================================================================
NB. Main Reference Books:
NB.
NB. Abramowitz, M. & Stegun, I.A. (1970)
NB. "Handbook of Mathematical Functions"
NB. Dover, New York [A&S]
NB.
NB. Fishman, G.S. (1996)
NB. 'Monte Carlo'
NB. Springer, Berlin [F]
NB.
NB. Press, W.H., Flannery, B.P., Teukolsky, S.A. & Vetterling, W.T. (1986)
NB. 'Numerical Recipes: The Art of Scientific Computing',
NB. Cambridge University Press, Cambridge [PFTV]
NB.
NB. ====================================================================
NB. Constants
NB.
RSQRT2PI=: % %: 2p1
HL2PI=: - ^. RSQRT2PI
NB. ====================================================================
NB. Utilities
NB.
NB. merge see "J Phrases" 3B
NB. underravel apply verb x. to ravelled y., then reshape
NB.
merge=: /:@/:@[ { ]
underravel=: adverb def '$@] $ x.@,'
NB. ====================================================================
NB. Standard mathematical functions
NB.
NB. beta beta function Beta (x. , y.)
NB. erf error function
NB. erfc complementary error function
NB. gamma gamma function Gamma (y.)
NB. ibp incomplete beta function (int from 0 to y.)
NB. ibq incomplete beta function (int from y. to 1)
NB. igp incomplete gamma function (int from -\infty to y.)
NB. igq incomplete gamma function (int from y. to +\infty)
NB. logbeta log(beta) OK for large arguments
NB. logcomb ^. x.!y. OK for large arguments
NB. logfact ^. !y. OK for large arguments
NB. loggamma log(gamma) OK for large arguments
NB. loggamma0 log(gamma) for scalar y.
NB.
NB. Subsidiary functions:
NB. exog '(^-y.) * (y.^x.) % gamma(x.)' OK for large arguments
NB. gt100 1 iff abs(y.) > 100
NB. ibcf continued fraction part of Incomplete Beta [PFTV sec.6.3]
NB. ibcfn ~ number of terms in continued fraction to evaluate ibcf
NB. igq1 continued fraction part of Gamma(x.,y.) [PFTV sec.6.2]
NB. igqn ~ number of terms in continued fraction to evaluate igq
NB. lgbig log(gamma) for large arguments
NB. x1mxob '(y.^{.x.) * ((1-y.)^{:x.) % ({.x.)beta({:x.)' OK for large x.
NB.
NB. Comments
NB. loggamma takes 2x as long as loggamma0 for scalar y.
NB. but is 5x as fast as loggamma"0 for y.= 90 + i. 2 3 4
NB.
gamma=: ! & <:
beta=: *&gamma % gamma@+
gt100=: | > 100"_
lgbig=: HL2PI"_ - ] - ([: % 12&*) - ([: % 360"_ * ^&3) - -&0.5 * ^.
loggamma0=: ((^.@gamma)`lgbig) @. (| > 100"_)
loggamma=: (] (] merge ^.@:gamma@:(#~ -.) , lgbig@:#~) gt100) underravel
logbeta=: +&:loggamma - loggamma@:+
logcomb=: [: - ^.@:>:@:] + [ logbeta&:>: -~
logfact=: loggamma@:>:
exog=: ([: ^ -@:] + (*^.) - loggamma@:[)`0: @. (] <: 0:)
igp=: dyad def '(x. exog"0 y.) * (((1 1 H. (1,x.)) y.) - 1) % y.'
igqn=: [: >. 4: + (60"_ + 16"_ * [: | [ - 4:) % ]
igq1=: [: (%`+)/ 1: , ] , [: , ([: }.@:i. igqn) ((-{.) ,. 1: ,. [ ,. {:@:]) ,
igq2=: exog * igq1
igq=: (-.@:igp)`igq2 @. ((> (+ 3: * %:))~)"0
erf=: 0.5"_ igp *:
erfc=: 0.5"_ igq *:
x1mxob=: [: ^ +/@:(* ^.@(,-.)) - logbeta/@:[
ibcfn=: [: >. 2: * [: %: 20"_ + >./
ibcf=: dyad define
'a b'=. x. [ m=. i. ibcfn x.
dodd=. -(a+m)*(a+b+m)*y. % (* >:) ap2m=. a (+ +:) m
deven=. m*(b-m)*y. % (* <:) ap2m
((x. x1mxob y.) % a) * (%`+)/ 1 , }. , deven,.1,.dodd,.1
)
ibp=: dyad define
if. y.<:0 do. 0
elseif. y.>:1 do. 1
elseif. 1 do.
if. y. < (>:@{. % >:@+/) x. do. x. ibcf y.
else. -. (|. x.) ibcf -. y.
end.
end.
)
ibq=: |.@[ ibp -.@]
NB. ====================================================================
NB. Pseudorandom Number Generation
NB.
NB. getrl get random link
NB. initrl initialize random link based on time+day
NB. (suffers from year 187213506 bug!)
NB. setrl set random link [y. integer; 0 < y. < _2+2^31]
NB.
getrl=: 9!:0
setrl=: 9!:1
initrl=: setrl @ <. @ ([: (0 60 60 24 31 12"_ #. |.) 6!:0)
NB. ====================================================================
NB. Scaling
NB.
NB. lrscale scale y. where x. = (left lim, right lim)
NB. lsfmv location,scale from mean,variance
NB. lsscale scale y. where x. = (location, scale)
NB. lssd (adverb) location/scale tran for c.d.f.
NB. lssf (adverb) location/scale tran for p.d.f.
NB. lssq (adverb) location/scale tran for quantiles (inverse c.d.f.)
NB.
lrscale=: {.@[ - -/@[ * ]
lsscale=: {.@[ + {:@[ * ]
lssd=: @ ((] - {.@[) % {:@[)
lssf=: lssd % {:@[
lssq=: adverb define
:
({. x.) + ({: x.) * u. y.
)
lsfmv=: ^&1 0.5
NB. ====================================================================
NB. Approximations to N(0,1) c.d.f. and inverse c.d.f.
NB.
NB. n01u upper tail probability of N(0,1) [PFTV]
NB. n01q1 approximate N(0,1) quantiles [A&S]
NB.
n01u=: -: @: erfc @: % %:@:2:
n01q1=: monad define
s=. * y. - 0.5 [ t=. %: _2 * ^. (<. -. ) y.
s * t - (2.515517 0.802853 0.010328 p. t) % 1 1.432788 0.189269 0.001308 p. t
)
NB. ====================================================================
NB. Univariate Discrete Distributions
NB.
NB. bin* binomial, n={.x., p={:x. [PFTV 6.3]
NB. pois* Poisson, mean=x. [PFTV 6.2]
NB.
NB. Subsidiary functions:
NB. binf[01] binf for [small,large] n={.x. or y.
NB. poisf[01] poisf for [small,large] x. or y.
NB.
binc=: (>:@] , {.@[ - ]) ibq {:@[
binf0=: (! {.)~ * (^~ {:)~ * -.@{:@[ ^ (-~ {.)~
binf1=: [: ^ (logcomb {.)~ + (* ^.@{:)~ + (-~ {.)~ * ^.@-.@{:@[
bins=: (>:@] , {.@[ - ]) ibp {:@[
poisc=: igq~
poisf0=: ^ * ^@-@[ % !@]
poisf1=: [: ^ (* ^.)~ - (+ logfact)
poisf=: poisf0`poisf1 @. ([: +./ , > 140"_)
poiss=: igp~
NB. ====================================================================
NB. Univariate Continuous Distributions
NB.
NB. beta* beta (default 1 1)
NB. cauchy* cauchy (i.e. 1&t)
NB. chisq* chi-squared (x. = df; default 1) [PFTV 6.2]
NB. exp* exponential, mean %x.
NB. f* F (x. = dfs; default 1 1) [PFTV 6.3]
NB. gamma* gamma (default 1 1)
NB. norm* Normal (x. = mean,var; default 0 1)
NB. t* t (x. = df; default 1) [PFTV 6.3]
NB. unif* uniform (rectangular) - U(0,1) by default
NB.
NB. normq uses 2 Newton-Raphson iterations
NB.
cauchyc=: (0.5"_ + _3&o. % 1p1"_) : (cauchyc lssd)
cauchys=: -.@cauchyc
chisqc=: -:@[ igp -:@]
chisqf=: verb def ('1 chisqf y.' ; ':' ; '(-: x.,1) gammaf y.')
chisqs=: -:@[ igq -:@]
expr=: -@^.@unifr : (%~ -@^.@unifr)
fc=: -:@|.@[ ibq {:@[ % [: +/ (* ,&1)
ff=: verb define
1 1 ff y.
:
'm2 n2'=. -: x.
y=. 0 >. y.
c=. (m2^m2) * (n2^n2) % m2 beta n2
(y.>:0) * c * (y ^ m2-1) % (n2+m2*y)^m2+n2
)
fs=: -:@|.@[ ibp {:@[ % [: +/ (* ,&1)
gammaf=: verb define
1 1 gammaf y.
:
'a b'=. x.
y=. 0 >. y.
c=. (b^a) % gamma a
(y.>:0) * c * (y^a-1) * ^ -b*y.
)
normc=: (0.5"_ + * * 0.5"_ - n01u@|) : (lsfmv@[ normc lssd ])
normf=: (RSQRT2PI&* @ ^ @ - @ -: @ *:) : (lsfmv@[ normf lssf ])
normq=: ((] + (- normd) % normf@])^:2 n01q1) : (lsfmv@[ normq lssq ])
ttail=: [: -: ([: -: [ , 1:) ibp [ % (+ *:)
tc=: ttail`(1:-ttail) @. (] > 0:)
tf=: verb define
1 tf y.
:
c=. (gamma 0.5*x.+1) % (%: o. x.) * gamma 0.5*x.
c * (1 + (*: y.)%x.) ^ -0.5*x.+1
)
ts=: (1:-ttail)`ttail @. (] > 0:)
unifc=: (0: >. <.&1) : (0: >. 1: <. (- {.)~ % -@:-/@[)
uniff=: (>&0 *. <&1) : (((> {.) *. (< {:))~ % -@-/@[)
unifq=: ] : (lrscale unifq)
unifr=: (([: >:@:? $&2147483647) % 2147483648"_) : (lrscale unifr)
unifs=: 1: - unifc
NB. ====================================================================
NB. TO DO
NB. speed up logcomb etc.
NB. Make cauchy[cs] more accurate for large args
--------------------------------------------------------------------------------------------
Ewart Shaw (Dr. J.E.H.Shaw) Department of Statistics, University of
Warwick.
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
