65 subroutine cdfgam(x,alpha,eps,iflag,cdfx)
69 real(kind=kind_phys),
intent(out) :: cdfx
70 real(kind=kind_phys),
intent(in) :: x, alpha, eps
72 integer iflag,i,j,k, imax
74 real(kind=kind_phys) dx, dgln, p,u,epsx,pdfl, eta, bl, uflo
75 data imax, uflo / 5000, 1.0e-37 /
89 if (alpha.le.uflo.or.eps.le.uflo)
then
100 pdfl = (alpha-1.0)*log(dx)-dx-dgln
101 if (pdfl.lt.log(uflo))
then
102 if (x.ge.alpha) cdfx = 1.0
109 if (p.le.real(k)) k = k-1
112 bl = (eta-1)*log(dx)-dx-dgln
118 if (u.le.epsx*(p-x))
return
127 if (u.le.epsx*(x-p))
continue
142 real(kind=kind_phys),
intent(in) :: x
143 real(kind=kind_phys),
intent(out) :: dgamlnout
146 real(kind=kind_phys) absacc, b1, b2, b3, b4, b5, b6, b7, b8
147 real(kind=kind_phys) c, dx, q, r, xmin, xn
148 data xmin, absacc / 6.894d0, 1.0e-15 /
149 data c / 0.918938533204672741780329736d0 /
150 data b1 / 0.833333333333333333333333333d-1 /
151 data b2 / - 0.277777777777777777777777778d-2 /
152 data b3 / 0.793650793650793650793650794d-3 /
153 data b4 / - 0.595238095238095238095238095d-3 /
154 data b5 / 0.841750841750841750841750842d-3 /
155 data b6 / - 0.191752691752691752691752692d-2 /
156 data b7 / 0.641025641025641025641025641d-2 /
157 data b8 / - 0.295506535947712418300653595d-1 /
159 if (x.le.0.0) stop
'*** x<=0.0 in function dgamln ***'
161 n = max(0,int(xmin - dx + 1.0d0) )
165 dgamlnout = r*( b1+q*( b2+q*( b3+q*( b4+q*( b5+q*( b6+q*( b7+q*b8 ) ) ) ) ) ) ) +c + (xn-0.5d0)*log(xn)-xn
172 dgamlnout = dgamlnout-log(q)
175 if (dgamlnout + absacc.eq.dgamlnout)
then
176 print *,
' ********* WARNING FROM FUNCTION DGAMLN *********'
177 print *,
' REQUIRED ABSOLUTE ACCURACY NOT ATTAINED FOR X = ',x
188 real(kind=kind_phys),
intent(in) :: pr, p, q
189 real(kind=kind_phys),
intent(out) :: x
191 integer iflag, iter, itmax
192 real(kind=kind_phys) tol, a, b, fa, fb, fc, cdf, tol1
193 real(kind=kind_phys) c, d, e, xm, s, u, v, r, eps
194 data itmax, eps / 50, 1.0e-12 /
212 if (pr.lt.0.0.or.pr.gt.1.)
then
216 if(min(p,q).le.0.)
then
220 if (tol.lt.1.0e-8)
then
228 if (fb*fa.gt.0.0)
then
235 if (fb*fc.gt.0.)
then
241 if (abs(fc).lt.abs(fb))
then
250 tol1 = 2.*eps*abs(b)+0.5*tol
252 if (abs(xm).le.tol1.or.fb.eq.0.0)
then
256 if (abs(e).ge.tol1.and.abs(fa).gt.abs(fb))
then
264 u = s*(2.0*xm*v*(v-r)-(b-a)*(r-1.0))
265 v = (v-1.0)*(r-1.0)*(s-1.0)
269 if (2.0*u.lt.min(3.0*xm*v-abs(tol1*v),abs(e*v)))
then
285 if (abs(d).gt.tol1)
then
290 call cdfbet(b,p,q,eps,iflag,cdf)
291 if (iflag.ne.0)
return
314 real(kind=kind_phys),
intent(in) :: x, p, q, eps
315 real(kind=kind_phys),
intent(out) :: cdfx
317 integer iflag, jmax, j
319 real(kind=kind_phys) dp, dq, gamln, yxeps, w, uflo
320 real(kind=kind_phys) xy, yx, pq, qp, pdfl, u, r, v
321 real(kind=kind_phys) tmp
322 data jmax, w, uflo / 5000, 20.0, 1.0e-30 /
326 if (p.le.uflo.or.q.le.uflo.or.eps.le.uflo)
then
335 ll = (p+w).ge.(p+q+2.0*w)*x
349 dp = (pq-1.)*log(xy)-tmp
351 dq = (qp-1.)*log(yx)-tmp
355 if (pdfl.ge.log(uflo))
then
359 if (u.le.eps*(1.-(pq+qp)*xy/(pq+1.)))
then
360 if (.not.ll) cdfx = 1.-cdfx
372 if (.not.ll) cdfx = 1.-cdfx
377 v = (pq+qp-1.)*xy*v/pq
381 if (.not.ll) cdfx = 1.-cdfx
395 real(kind=kind_phys),
intent(in) :: x
396 real(kind=kind_phys),
intent(out) :: y
399 real(kind=kind_phys) absacc, b1, b2, b3, b4, b5, b6, b7, b8
400 real(kind=kind_phys) c, dx, q, r, xmin, xn
402 data xmin, absacc / 1.357d0, 1.0e-3 /
403 data c / 0.918938533204672741780329736d0 /
404 data b1 / 0.833333333333333333333333333d-1 /
405 data b2 / - 0.277777777777777777777777778d-2 /
406 data b3 / 0.793650793650793650793650794d-3 /
407 data b4 / - 0.595238095238095238095238095d-3 /
408 data b5 / 0.841750841750841750841750842d-3 /
409 data b6 / - 0.191752691752691752691752692d-2 /
410 data b7 / 0.641025641025641025641025641d-2 /
411 data b8 / - 0.295506535947712418300653595d-1 /
413 if (x.le.0.0) stop
'*** x<=0.0 in function gamln ***'
415 n = max(0,int(xmin - dx + 1.0d0) )
419 y = r*( b1+q*( b2+q*( b3+q*( b4+q*( b5+q*( b6+q*( b7+q*b8 ) &
420 & )) ) ) ) ) +c + (xn-0.5d0)*log(xn)-xn
430 if (y + absacc.eq.y)
then
431 print *,
' ********* WARNING FROM FUNCTION GAMLN *********'
432 print *,
' REQUIRED ABSOLUTE ACCURACY NOT ATTAINED FOR X = ',x
subroutine, public ppfbet(pr, p, q, iflag, x)
This subroutine computes the beta distribution value that matches the percentile from the random patt...
subroutine cdfbet(x, p, q, eps, iflag, cdfx)
This subroutine computes the value of the cumulative beta distribution at a single point x,...