Nice. Thank you,
-- Raul On Wed, Jun 10, 2015 at 9:36 AM, Shaw, Ewart <[email protected]> wrote: > 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 ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
