64 subroutine cdfgam(x,alpha,eps,iflag,cdfx)
68 real(kind=kind_phys),
intent(out) :: cdfx
69 real(kind=kind_phys),
intent(in) :: x, alpha, eps
71 integer iflag,i,j,k, imax
73 real(kind=kind_phys) dx, dgln, p,u,epsx,pdfl, eta, bl, uflo
74 data imax, uflo / 5000, 1.0e-37 /
88 if (alpha.le.uflo.or.eps.le.uflo)
then
99 pdfl = (alpha-1.0)*log(dx)-dx-dgln
100 if (pdfl.lt.log(uflo))
then
101 if (x.ge.alpha) cdfx = 1.0
108 if (p.le.real(k)) k = k-1
111 bl = (eta-1)*log(dx)-dx-dgln
117 if (u.le.epsx*(p-x))
return
126 if (u.le.epsx*(x-p))
continue
141 real(kind=kind_phys),
intent(in) :: x
142 real(kind=kind_phys),
intent(out) :: dgamlnout
145 real(kind=kind_phys) absacc, b1, b2, b3, b4, b5, b6, b7, b8
146 real(kind=kind_phys) c, dx, q, r, xmin, xn
147 data xmin, absacc / 6.894d0, 1.0e-15 /
148 data c / 0.918938533204672741780329736d0 /
149 data b1 / 0.833333333333333333333333333d-1 /
150 data b2 / - 0.277777777777777777777777778d-2 /
151 data b3 / 0.793650793650793650793650794d-3 /
152 data b4 / - 0.595238095238095238095238095d-3 /
153 data b5 / 0.841750841750841750841750842d-3 /
154 data b6 / - 0.191752691752691752691752692d-2 /
155 data b7 / 0.641025641025641025641025641d-2 /
156 data b8 / - 0.295506535947712418300653595d-1 /
158 if (x.le.0.0) stop
'*** x<=0.0 in function dgamln ***'
160 n = max(0,int(xmin - dx + 1.0d0) )
164 dgamlnout = r*( b1+q*( b2+q*( b3+q*( b4+q*( b5+q*( b6+q*( b7+q*b8 ) ) ) ) ) ) ) +c + (xn-0.5d0)*log(xn)-xn
171 dgamlnout = dgamlnout-log(q)
174 if (dgamlnout + absacc.eq.dgamlnout)
then
175 print *,
' ********* WARNING FROM FUNCTION DGAMLN *********'
176 print *,
' REQUIRED ABSOLUTE ACCURACY NOT ATTAINED FOR X = ',x
187 real(kind=kind_phys),
intent(in) :: pr, p, q
188 real(kind=kind_phys),
intent(out) :: x
190 integer iflag, iter, itmax
191 real(kind=kind_phys) tol, a, b, fa, fb, fc, cdf, tol1
192 real(kind=kind_phys) c, d, e, xm, s, u, v, r, eps
193 data itmax, eps / 50, 1.0e-12 /
211 if (pr.lt.0.0.or.pr.gt.1.)
then
215 if(min(p,q).le.0.)
then
219 if (tol.lt.1.0e-8)
then
227 if (fb*fa.gt.0.0)
then
234 if (fb*fc.gt.0.)
then
240 if (abs(fc).lt.abs(fb))
then
249 tol1 = 2.*eps*abs(b)+0.5*tol
251 if (abs(xm).le.tol1.or.fb.eq.0.0)
then
255 if (abs(e).ge.tol1.and.abs(fa).gt.abs(fb))
then
263 u = s*(2.0*xm*v*(v-r)-(b-a)*(r-1.0))
264 v = (v-1.0)*(r-1.0)*(s-1.0)
268 if (2.0*u.lt.min(3.0*xm*v-abs(tol1*v),abs(e*v)))
then
284 if (abs(d).gt.tol1)
then
289 call cdfbet(b,p,q,eps,iflag,cdf)
290 if (iflag.ne.0)
return
313 real(kind=kind_phys),
intent(in) :: x, p, q, eps
314 real(kind=kind_phys),
intent(out) :: cdfx
316 integer iflag, jmax, j
318 real(kind=kind_phys) dp, dq, gamln, yxeps, w, uflo
319 real(kind=kind_phys) xy, yx, pq, qp, pdfl, u, r, v
320 real(kind=kind_phys) tmp
321 data jmax, w, uflo / 5000, 20.0, 1.0e-30 /
325 if (p.le.uflo.or.q.le.uflo.or.eps.le.uflo)
then
334 ll = (p+w).ge.(p+q+2.0*w)*x
348 dp = (pq-1.)*log(xy)-tmp
350 dq = (qp-1.)*log(yx)-tmp
354 if (pdfl.ge.log(uflo))
then
358 if (u.le.eps*(1.-(pq+qp)*xy/(pq+1.)))
then
359 if (.not.ll) cdfx = 1.-cdfx
371 if (.not.ll) cdfx = 1.-cdfx
376 v = (pq+qp-1.)*xy*v/pq
380 if (.not.ll) cdfx = 1.-cdfx
394 real(kind=kind_phys),
intent(in) :: x
395 real(kind=kind_phys),
intent(out) :: y
398 real(kind=kind_phys) absacc, b1, b2, b3, b4, b5, b6, b7, b8
399 real(kind=kind_phys) c, dx, q, r, xmin, xn
401 data xmin, absacc / 1.357d0, 1.0e-3 /
402 data c / 0.918938533204672741780329736d0 /
403 data b1 / 0.833333333333333333333333333d-1 /
404 data b2 / - 0.277777777777777777777777778d-2 /
405 data b3 / 0.793650793650793650793650794d-3 /
406 data b4 / - 0.595238095238095238095238095d-3 /
407 data b5 / 0.841750841750841750841750842d-3 /
408 data b6 / - 0.191752691752691752691752692d-2 /
409 data b7 / 0.641025641025641025641025641d-2 /
410 data b8 / - 0.295506535947712418300653595d-1 /
412 if (x.le.0.0) stop
'*** x<=0.0 in function gamln ***'
414 n = max(0,int(xmin - dx + 1.0d0) )
418 y = r*( b1+q*( b2+q*( b3+q*( b4+q*( b5+q*( b6+q*( b7+q*b8 ) &
419 & )) ) ) ) ) +c + (xn-0.5d0)*log(xn)-xn
429 if (y + absacc.eq.y)
then
430 print *,
' ********* WARNING FROM FUNCTION GAMLN *********'
431 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,...