Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This include file contains the implementation of procedures in [pm_mathBeta](@ref pm_mathBeta).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : !%%%%%%%%%%%%%%%%%
28 : #if getLogBeta_ENABLED
29 : !%%%%%%%%%%%%%%%%%
30 :
31 163787 : CHECK_ASSERTION(__LINE__, 0._RKC < beta, SK_"@getLogBeta(): The condition `0._RKC < beta` must hold. beta = "//getStr(beta)) ! fpp
32 163787 : CHECK_ASSERTION(__LINE__, 0._RKC < alpha, SK_"@getLogBeta(): The condition `0._RKC < alpha` must hold. alpha = "//getStr(alpha)) ! fpp
33 163787 : logFuncBeta = log_gamma(alpha) + log_gamma(beta) - log_gamma(alpha + beta)
34 :
35 : !%%%%%%%%%%%%%%%%%
36 : #elif getBetaInc_ENABLED
37 : !%%%%%%%%%%%%%%%%%
38 :
39 : real(RKC) , parameter :: reltol = epsilon(0._RKC)**.7, THRESH = 3000._RKC ! Based on the suggestion of Press et al (1992). This may need fine-tuning for modern hardware.
40 : logical(LK) :: dengis
41 : integer(IK) :: info
42 9056 : if (present(signed)) then
43 9 : dengis = signed
44 : else
45 9047 : dengis = .false._LK
46 : end if
47 9056 : if (alpha < THRESH .or. beta < THRESH) then
48 9056 : call setBetaInc(betaInc, x, alpha, beta, getLogBeta(alpha, beta), dengis, info)
49 : else
50 0 : info = 0_IK
51 : end if
52 9056 : if (info /= 0_IK) call setBetaInc(betaInc, x, alpha, beta, getLogBeta(alpha, beta), reltol, dengis, info)
53 : !if (betaInc < 0._RKC) betaInc = betaInc + 1._RKC
54 : if (info /= 0_IK) error stop MODULE_NAME//SK_"@getBetaInc(): Failed to converge in computing the continued fraction representation of the Beta Function." ! LCOV_EXCL_LINE
55 :
56 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57 : #elif setBetaInc_ENABLED && GK21_ENABLED
58 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
59 :
60 : real(RKC), parameter :: EPS = epsilon(0._RKC), lb = 0._RKC
61 : integer(IK), parameter :: nintmax = 100
62 : logical(LK) :: mirrored
63 : integer(IK) :: sindex(nintmax), neval, nint
64 : real(RKC) :: abserr, sinfo(4, nintmax), betaMinus1, alphaMinus1
65 2 : CHECK_ASSERTION(__LINE__, 0._RKC < beta, SK_"@setBetaInc(): The condition `0. < beta` must hold. beta = "//getStr(beta)) ! fpp
66 2 : CHECK_ASSERTION(__LINE__, 0._RKC < alpha, SK_"@setBetaInc(): The condition `0. < alpha` must hold. alpha = "//getStr(alpha)) ! fpp
67 2 : mirrored = signed .and. x < 0._RKC
68 0 : if (mirrored) x = -x
69 2 : if(0._RKC < x .and. x < 1._RKC) then
70 2 : if (.not. mirrored) then
71 2 : mirrored = (alpha + 1._RKC) / (alpha + beta + 2._RKC) < x
72 2 : if (mirrored) x = 1._RKC - x
73 : endif
74 2 : if (mirrored) then ! swap
75 : abserr = beta
76 0 : beta = alpha
77 0 : alpha = abserr
78 : endif
79 : ! Prepare the integrand parameters.
80 2 : betaMinus1 = beta - 1._RKC
81 2 : alphaMinus1 = alpha - 1._RKC
82 : ! integrate.
83 : info =getQuadErr( getFunc & ! LCOV_EXCL_LINE
84 : , lb = lb & ! LCOV_EXCL_LINE
85 : , ub = x & ! LCOV_EXCL_LINE
86 : , abstol = EPS & ! LCOV_EXCL_LINE
87 : , reltol = reltol & ! LCOV_EXCL_LINE
88 : , qrule = qrule & ! LCOV_EXCL_LINE
89 : , integral = betaInc & ! LCOV_EXCL_LINE
90 : , abserr = abserr & ! LCOV_EXCL_LINE
91 : , sinfo = sinfo & ! LCOV_EXCL_LINE
92 : , sindex = sindex & ! LCOV_EXCL_LINE
93 : , neval = neval & ! LCOV_EXCL_LINE
94 : , nint = nint & ! LCOV_EXCL_LINE
95 : , help = weps & ! LCOV_EXCL_LINE
96 2 : )
97 2 : if (mirrored) then
98 0 : if (signed) then
99 0 : betaInc = -betaInc
100 : else
101 0 : betaInc = 1._RKC - betaInc
102 : end if
103 : end if
104 0 : elseif (0._RKC == x .or. x == 1._RKC) then
105 0 : info = 0_IK
106 0 : betaInc = x
107 : else
108 0 : info = 1_IK
109 0 : CHECK_ASSERTION(__LINE__, info == 0_IK, SK_"@setBetaInc(): The condition `0. <= x .and. x <= 1.` must hold. x = "//getStr(x)) ! fpp
110 : endif
111 : contains
112 1722 : pure function getFunc(xx) result(func)
113 : real(RKC), intent(in) :: xx
114 : real(RKC) :: func
115 1722 : func = exp(log(xx) * alphaMinus1 + log(1._RKC - xx) * betaMinus1 - logFuncBeta)
116 1722 : end function
117 :
118 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119 : #elif setBetaInc_ENABLED && Def_ENABLED && 1
120 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121 :
122 : ! This implementation is partially inspired by the proposed approach of Numerical Recipes (1992) by the Great Bill Press et al.
123 : logical(LK) :: mirrored
124 : integer(IK) :: iter, iterTwice
125 : integer(IK) , parameter :: MAX_ITER = 10000_IK
126 : real(RKC) , parameter :: EPS = epsilon(0._RKC), FPMIN = tiny(0._RKC) / EPS
127 : REAL(RKC) :: aa, c, d, delta, sumAlphaBeta, alphamMinusOne, alphamPlusOne
128 728667 : CHECK_ASSERTION(__LINE__, 0._RKC < beta, SK_"@setBetaInc(): The condition `0._RKC < beta` must hold. beta = "//getStr(beta)) ! fpp
129 728667 : CHECK_ASSERTION(__LINE__, 0._RKC < alpha, SK_"@setBetaInc(): The condition `0._RKC < alpha` must hold. alpha = "//getStr(alpha)) ! fpp
130 728667 : mirrored = signed .and. x < 0._RKC
131 0 : if (mirrored) x = -x
132 728667 : if(0._RKC < x .and. x < 1._RKC) then
133 728599 : sumAlphaBeta = alpha + beta
134 728599 : if (.not. mirrored) then
135 728599 : mirrored = (alpha + 1._RKC) / (sumAlphaBeta + 2._RKC) < x
136 728599 : if (mirrored) x = 1._RKC - x
137 : endif
138 728599 : if (mirrored) then
139 : d = beta
140 65507 : beta = alpha
141 65507 : alpha = d
142 : endif
143 728599 : alphamPlusOne = alpha + 1._RKC
144 728599 : alphamMinusOne = alpha - 1._RKC
145 728599 : d = 1._RKC - sumAlphaBeta * x / alphamPlusOne
146 728599 : if (abs(d) < FPMIN) d = FPMIN
147 534422 : c = 1._RKC
148 728599 : d = 1._RKC / d
149 728599 : betaInc = d
150 2851931 : do iter = 1, MAX_ITER
151 2851931 : iterTwice = 2_IK * iter
152 2851931 : aa = iter * (beta - iter) * x / ((alphamMinusOne + iterTwice) * (alpha + iterTwice))
153 2851931 : d = 1._RKC + aa * d
154 2851931 : if (abs(d) < FPMIN) d = FPMIN
155 2851931 : c = 1._RKC + aa / c
156 2851931 : if (abs(c) < FPMIN) c = FPMIN
157 2851931 : d = 1._RKC / d
158 2851931 : betaInc = betaInc * d * c
159 2851931 : aa = -(alpha + iter) * (sumAlphaBeta + iter) * x / ((alpha + iterTwice) * (alphamPlusOne + iterTwice))
160 2851931 : d = 1._RKC + aa * d
161 2851931 : if (abs(d) < FPMIN) d = FPMIN
162 2851931 : c = 1._RKC + aa / c
163 2851931 : if (abs(c) < FPMIN) c = FPMIN
164 2851931 : d = 1._RKC / d
165 2851931 : delta = d * c
166 2851931 : betaInc = betaInc * delta
167 2851931 : if (EPS < abs(delta - 1._RKC)) cycle
168 : !betaInc = betaInc * exp(alpha * log(x) + beta * getLog1p(-x) - logFuncBeta) / alpha
169 728599 : betaInc = betaInc * exp(alpha * log(x) + beta * log(1._RKC - x) - logFuncBeta) / alpha
170 : !if (mirrored) betaInc = 1._RKC - betaInc ! 1._RKC is the regularization term.
171 : !if (mirrored) betaInc = -betaInc ! 1._RKC is the regularization term.
172 728599 : if (mirrored) then
173 65507 : if (signed) then
174 23880 : betaInc = -betaInc
175 : else
176 41627 : betaInc = 1._RKC - betaInc
177 : end if
178 : end if
179 728599 : info = 0_IK
180 2851931 : return
181 : end do
182 0 : info = 1_IK
183 68 : elseif (0._RKC == x .or. x == 1._RKC) then
184 68 : info = 0_IK
185 68 : betaInc = x
186 : else
187 0 : info = 1_IK
188 0 : CHECK_ASSERTION(__LINE__, info == 0_IK, SK_"@setBetaInc(): The condition `0. <= x .and. x <= 1.` must hold. x = "//getStr(x)) ! fpp
189 : endif
190 :
191 : !%%%%%%%%%%%%%%%%%%%%%%
192 : #elif setBetaInc_ENABLED && 1
193 : !%%%%%%%%%%%%%%%%%%%%%%
194 :
195 : ! This implementation is based on netlib BETAIN.
196 : ! As far as evidenced by tests, it requires many cycles for convergence for some precisions other than 64-bit.
197 : ! Its performance and accuracy against the alternative above is yet to be tested.
198 : ! KL Majumder, GP Bhattacharjee. Modifications by John Burkardt.
199 : ! Reference:
200 : ! KL Majumder, GP Bhattacharjee,
201 : ! Algorithm AS 63:
202 : ! The incomplete Beta Integral,
203 : ! Applied Statistics,
204 : ! Volume 22, Number 3, 1973, pages 409-411.
205 : integer(IK) :: ns, iter
206 : logical(LK) :: mirrored
207 : real(RKC) :: ai, sumAlphaBeta, rx, temp, term
208 : real(RKC) , parameter :: EPS = epsilon(0._RKC)**(7._RKC/8._RKC)
209 : integer(IK) , parameter :: MAX_ITER = 10000_IK
210 : CHECK_ASSERTION(__LINE__, 0._RKC < beta, SK_"@setBetaInc(): The condition `0._RKC < beta` must hold. beta = "//getStr(beta)) ! fpp
211 : CHECK_ASSERTION(__LINE__, 0._RKC < alpha, SK_"@setBetaInc(): The condition `0._RKC < alpha` must hold. alpha = "//getStr(alpha)) ! fpp
212 : mirrored = signed .and. x < 0._RKC
213 : if (mirrored) x = -x
214 : if (0._RKC < x .and. x < 1._RKC) then
215 :
216 : ! Change tail if necessary.
217 :
218 : sumAlphaBeta = alpha + beta
219 : if (.not. mirrored) then
220 : !mirrored = (alpha + 1._RKC) / (sumAlphaBeta + 2._RKC) < x ! this is the original netlib condition.
221 : mirrored = alpha < sumAlphaBeta * x ! This is the implementation by press et al.
222 : if (mirrored) x = 1._RKC - x
223 : endif
224 : if (mirrored) then
225 : temp = beta
226 : beta = alpha
227 : alpha = temp
228 : endif
229 :
230 : ai = 1._RKC
231 : term = 1._RKC
232 : betaInc = 1._RKC
233 : ns = int(beta + (1._RKC - x) * sumAlphaBeta, IK)
234 : ! Use the Soper reduction formula.
235 : rx = x / (1._RKC - x)
236 : temp = beta - ai
237 : if (ns == 0_IK) rx = x
238 : do iter = 1, MAX_ITER
239 : term = term * temp * rx / (alpha + ai)
240 : betaInc = betaInc + term
241 : temp = abs(term)
242 : if (temp <= EPS .and. temp <= EPS * betaInc) then
243 : !betaInc = betaInc * exp(alpha * log(x) + (beta - 1._RKC) * getLog1p(-x) - logFuncBeta) / alpha
244 : betaInc = betaInc * exp(alpha * log(x) + (beta - 1._RKC) * log(1._RKC - x) - logFuncBeta) / alpha
245 : if (mirrored) then
246 : if (signed) then
247 : betaInc = -betaInc
248 : else
249 : betaInc = 1._RKC - betaInc
250 : end if
251 : end if
252 : info = 0_IK
253 : return
254 : end if
255 : ai = ai + 1._RKC
256 : ns = ns - 1
257 : if (ns < 0_IK) then
258 : temp = sumAlphaBeta
259 : sumAlphaBeta = sumAlphaBeta + 1._RKC
260 : else
261 : temp = beta - ai
262 : if (ns == 0_IK) rx = x
263 : end if
264 : end do
265 : info = 1_IK
266 : elseif (0._RKC == x .or. x == 1._RKC) then
267 : info = 0_IK
268 : betaInc = x
269 : else
270 : info = 1_IK
271 : CHECK_ASSERTION(__LINE__, info == 0_IK, SK_"@setBetaInc(): The condition `0. <= x .and. x <= 1.` must hold. x = "//getStr(x)) ! fpp
272 : endif
273 :
274 : !%%%%%%%%%%%%%%%%%
275 : #elif getBetaInv_ENABLED
276 : !%%%%%%%%%%%%%%%%%
277 :
278 : integer(IK) :: info
279 : logical(LK) :: dengis
280 50009 : if (present(signed)) then
281 30000 : dengis = signed
282 : else
283 20009 : dengis = .false._LK
284 : end if
285 50009 : call setBetaInv(betaInv, betaInc, alpha, beta, getLogBeta(alpha, beta), dengis, info)
286 : if (info /= 0_IK) error stop MODULE_NAME//SK_"@getBetaInv(): Failed to converge in computing the Regularized Inverse Incomplete Beta Function." ! LCOV_EXCL_LINE
287 :
288 : !%%%%%%%%%%%%%%%%%%%%%%
289 : #elif setBetaInv_ENABLED && 1
290 : !%%%%%%%%%%%%%%%%%%%%%%
291 :
292 : ! This implementation follows the approach of:
293 : ! GW Cran, KJ Martin, GE Thomas,
294 : ! Remark AS R19 and Algorithm AS 109:
295 : ! A Remark on Algorithms AS 63: The Incomplete Beta Integral
296 : ! and AS 64: Inverse of the Incomplete Beta Integeral,
297 : ! Applied Statistics,
298 : ! Volume 26, Number 1, 1977, pages 111-114.
299 : logical(LK) :: mirrored
300 : real(RKC) , parameter :: underflow = tiny(0._RKC)
301 : real(RKC) , parameter :: experr = log10(underflow)
302 : real(RKC) , parameter :: ONE_THIRD = 1._RKC / 3._RKC
303 : real(RKC) , parameter :: ONE_SIXTH = 1._RKC / 6._RKC
304 : real(RKC) , parameter :: FIVE_SIXTH = 5._RKC / 6._RKC
305 : real(RKC) , parameter :: ONE_MEPS = 1._RKC - sqrt(epsilon(0._RKC))
306 : real(RKC) :: tolerance, adj, g, h, prev, r, s, sq, t, tx, w, betaIncOld, betaIncNew
307 70018 : CHECK_ASSERTION(__LINE__, 0._RKC < alpha, SK_"@setBetaInv(): The condition `0. < alpha` must hold. alpha = "//getStr(alpha)) ! fpp
308 70018 : CHECK_ASSERTION(__LINE__, 0._RKC < beta, SK_"@setBetaInv(): The condition `0. < beta` must hold. alpha = "//getStr(beta)) ! fpp
309 70018 : mirrored = signed .and. betaInc < 0._RKC
310 0 : if (mirrored) betaInc = -betaInc
311 70018 : if (0._RKC < betaInc .and. betaInc < 1._RKC) then
312 :
313 : ! Change tail if necessary.
314 :
315 70010 : if (mirrored) then
316 : r = alpha
317 0 : alpha = beta
318 0 : beta = r
319 : else
320 70010 : mirrored = 0.5_RKC < betaInc
321 70010 : if (mirrored) then
322 35002 : betaInc = 1._RKC - betaInc
323 : r = alpha
324 35002 : alpha = beta
325 35002 : beta = r
326 : end if
327 : end if
328 :
329 : ! Compute the initial approximation.
330 :
331 70010 : r = sqrt(-log(betaInc**2))
332 70010 : betaIncNew = r - (2.30753_RKC + 0.27061_RKC * r) / (1._RKC + (0.99229_RKC + 0.04481_RKC * r) * r)
333 70010 : if (1._RKC < alpha .and. 1._RKC < beta) then
334 14007 : r = (betaIncNew * betaIncNew - 3._RKC) * ONE_SIXTH
335 14007 : s = 1._RKC / (alpha + alpha - 1._RKC)
336 14007 : t = 1._RKC / (beta + beta - 1._RKC)
337 14007 : h = 2._RKC / (s + t)
338 14007 : w = betaIncNew * sqrt(h + r) / h - (t - s) * (r + FIVE_SIXTH - 2._RKC / (3._RKC * h))
339 14007 : betaInv = alpha / (alpha + beta * exp(w + w))
340 : else
341 56003 : r = beta + beta
342 56003 : t = 1._RKC / (9._RKC * beta)
343 56003 : t = r * (1._RKC - t + betaIncNew * sqrt(t))**3
344 56003 : if (t <= 0._RKC) then
345 2352 : betaInv = 1._RKC - exp((log((1._RKC - betaInc) * beta) + logFuncBeta) / beta)
346 : else
347 53651 : t = (4._RKC * alpha + r - 2._RKC) / t
348 53651 : if (t <= 1._RKC) then
349 30521 : betaInv = exp((log(betaInc * alpha) + logFuncBeta) / alpha)
350 : else
351 23130 : betaInv = 1._RKC - 2._RKC / (t + 1._RKC)
352 : end if
353 : end if
354 : end if
355 :
356 : ! Solve for X via modified Newton-Raphson method, using the function `setBetaInc()`.
357 :
358 70010 : t = 1._RKC - beta
359 70010 : r = 1._RKC - alpha
360 50000 : betaIncOld = 0._RKC
361 50000 : prev = 1._RKC
362 50000 : sq = 1._RKC
363 70010 : if (betaInv < 0.0001_RKC) betaInv = 0.0001_RKC
364 70010 : if (0.9999_RKC < betaInv) betaInv = 0.9999_RKC
365 70010 : tolerance = 10._RKC**int(max(-5._RKC / alpha**2 - 1._RKC / betaInc**0.2 - 13._RKC, experr), IK)
366 :
367 : ! Begin iteration.
368 :
369 640579 : loopFindRoot: do ! 10
370 710589 : call setBetaInc(betaIncNew, betaInv, alpha, beta, logFuncBeta, signed, info)
371 : if (info /= 0_IK) return ! LCOV_EXCL_LINE
372 710589 : if (betaIncNew < 0._RKC) then
373 : !if (abs(betaIncNew) < epsilon(0._RKC)) error stop "Underflow detected."
374 23880 : betaIncNew = betaIncNew - betaInc
375 : !if (abs(betaIncNew) < epsilon(0._RKC)) error stop "Underflow detected."
376 23880 : betaIncNew = betaIncNew + 1._RKC
377 : else
378 686709 : betaIncNew = betaIncNew - betaInc
379 : end if
380 : !betaIncNew = betaIncNew * exp(logFuncBeta + r * log(betaInv) + t * getLog1p(-betaInv))
381 710589 : betaIncNew = betaIncNew * exp(logFuncBeta + r * log(betaInv) + t * log(1._RKC - betaInv))
382 710589 : if (betaIncNew * betaIncOld <= 0._RKC) prev = max(sq, underflow)
383 520420 : g = 1._RKC
384 : loopRefine: do ! 20
385 1480693 : adj = g * betaIncNew
386 1480693 : sq = adj * adj
387 1480693 : if (sq < prev) then
388 1118986 : tx = betaInv - adj
389 1118986 : if (tx < 0._RKC .or. 1._RKC < tx) then
390 408397 : g = g * ONE_THIRD
391 408397 : cycle loopRefine
392 : end if
393 : else ! 30
394 361707 : g = g * ONE_THIRD
395 361707 : cycle loopRefine
396 : end if
397 : ! Check whether current estimate is acceptable. The change "VALUE = TX" was suggested by Ivan Ukhov.
398 710589 : if (prev <= tolerance .or. betaIncNew * betaIncNew <= tolerance) then ! 40
399 50070 : betaInv = tx
400 50070 : exit loopFindRoot
401 : end if
402 660519 : if (tx /= 0._RKC .and. tx /= 1._RKC) exit loopRefine
403 0 : g = g * ONE_THIRD
404 : end do loopRefine
405 660519 : if (betaInv == tx) exit loopFindRoot
406 170159 : betaInv = tx
407 476640 : betaIncOld = betaIncNew
408 : end do loopFindRoot
409 70010 : if (mirrored) then
410 35002 : if (signed) then
411 15000 : betaInv = -betaInv
412 : else
413 20002 : betaInv = 1._RKC - betaInv
414 : end if
415 : end if
416 :
417 8 : elseif (0._RKC == betaInc .or. betaInc == 1._RKC) then
418 :
419 8 : info = 0_IK
420 8 : betaInv = betaInc
421 :
422 : else
423 :
424 0 : info = 1_IK
425 0 : CHECK_ASSERTION(__LINE__, info == 0_IK, SK_"@setBetaInc(): The condition `0. <= betaInc .and. betaInc <= 1.` must hold. betaInc = "//getStr(betaInc)) ! fpp
426 :
427 : end if
428 :
429 : !%%%%%%%%%%%%%%%%%%%%%%
430 : #elif setBetaInv_ENABLED && 0
431 : !%%%%%%%%%%%%%%%%%%%%%%
432 :
433 : ! This implementation follows the proposed approach of Numerical Recipes by Press et al. (2007) using Halley approximation.
434 : ! No sign corrections are implemented in this algorithm. As such the output is prone to numerical roundup to `1`.
435 : integer(IK) :: iter
436 : real(RKC) , parameter :: EPS = sqrt(epsilon(0._RKC))
437 : real(RKC) :: pp, t, u, betaIncNew, al, h, w, alphaMinus1, betaMinus1
438 : if (signed .and. betaInc < 0._RKC) betaInc = 1 + betaInc
439 : if (0._RKC < betaInc .and. betaInc < 1._RKC) then
440 : betaMinus1 = beta - 1._RKC
441 : alphaMinus1 = alpha - 1._RKC
442 : if (alpha < 1._RKC .or. beta < 1._RKC) then
443 : u = exp(beta * log(beta / (alpha + beta))) / beta
444 : t = exp(alpha * log(alpha / (alpha + beta))) / alpha
445 : w = t + u
446 : if (betaInc < t / w) then
447 : betaInv = (alpha * w * betaInc)**(1._RKC / alpha)
448 : else
449 : betaInv = 1._RKC - (beta * w * (1._RKC - betaInc))**(1._RKC / beta)
450 : end if
451 : else
452 : if (betaInc < 0.5_RKC) then
453 : pp = betaInc
454 : else
455 : pp = 1._RKC - betaInc
456 : end if
457 : t = sqrt(-2._RKC * log(pp))
458 : betaInv = (2.30753_RKC + t * 0.27061_RKC) / (1._RKC + t * (0.99229_RKC + t * 0.04481_RKC)) - t
459 : if (betaInc < .5_RKC) betaInv = -betaInv
460 : al = (betaInv**2 - 3._RKC) / 6._RKC
461 : h = 2._RKC / (1._RKC / (2._RKC * alpha - 1._RKC) + 1._RKC / (2._RKC * beta - 1._RKC))
462 : w = (betaInv * sqrt(al + h) / h) - (1._RKC / (2._RKC * beta - 1._RKC) - 1._RKC / (2._RKC * alpha - 1._RKC)) * (al + 5._RKC / 6._RKC - 2._RKC / (3._RKC * h))
463 : betaInv = alpha / (alpha + beta * exp(2._RKC * w))
464 : end if
465 : loopFindRoot: do iter = 1, 10
466 : if (betaInv == 0._RKC .or. betaInv == 1._RKC) exit loopFindRoot
467 : call setBetaInc(betaIncNew, betaInv, alpha, beta, logFuncBeta, signed, info)
468 : if (info /= 0_IK) return
469 : betaIncNew = betaIncNew - betaInc
470 : t = exp(alphaMinus1 * log(betaInv) + betaMinus1 * log(1._RKC - betaInv) - logFuncBeta)
471 : u = betaIncNew / t
472 : t = u / (1._RKC - .5_RKC * min(1._RKC, u * (alphaMinus1 / betaInv - betaMinus1 / (1._RKC - betaInv))))
473 : betaInv = betaInv - t
474 : if (betaInv <= 0._RKC) betaInv = .5_RKC * (betaInv + t)
475 : if (betaInv >= 1._RKC) betaInv = .5_RKC * (betaInv + t + 1._RKC)
476 : if (abs(t) < EPS * betaInv .and. 1_IK < iter) exit loopFindRoot
477 : end do loopFindRoot
478 : info = 0_IK
479 : elseif (0._RKC == betaInc .or. betaInc == 1._RKC) then
480 : info = 0_IK
481 : betaInv = betaInc
482 : else
483 : info = 1_IK
484 : CHECK_ASSERTION(__LINE__, .false., SK_"@setBetaInv(): The condition `0. <= betaInc .and. betaInc <= 1.` must hold. betaInc = "//getStr(betaInc)) ! fpp
485 : end if
486 :
487 : !%%%%%%%%%%%%%%%%%%%%%%
488 : #elif setBetaInv_ENABLED && 0
489 : !%%%%%%%%%%%%%%%%%%%%%%
490 :
491 : ! Naive root-finding implementation.
492 : ! This implementation is exploratory and should never be used.
493 : integer(IK) :: neval
494 : logical(LK) :: mirrored
495 : real(RKC) :: lb, ub, lf, uf, abstol, h, w, r, s, t, betaIncNew
496 : real(RKC) , parameter :: ONE_SIXTH = 1._RKC / 6._RKC
497 : !real(RKC) , parameter :: abstol = epsilon(0._RKC)
498 : CHECK_ASSERTION(__LINE__, 0._RKC < alpha, SK_"@setBetaInv(): The condition `0. < alpha` must hold. alpha = "//getStr(alpha)) ! fpp
499 : CHECK_ASSERTION(__LINE__, 0._RKC < beta, SK_"@setBetaInv(): The condition `0. < beta` must hold. alpha = "//getStr(beta)) ! fpp
500 : info = 0_IK
501 : mirrored = signed .and. betaInc < 0._RKC
502 : if (mirrored) betaInc = -betaInc
503 : if (0._RKC < betaInc .and. betaInc < 1._RKC) then
504 : if (.not. mirrored) then
505 : mirrored = 0.5_RKC < betaInc
506 : if (mirrored) betaInc = 1._RKC - betaInc
507 : end if
508 : if (mirrored) then
509 : r = alpha
510 : alpha = beta
511 : beta = r
512 : end if
513 : ! Compute the initial approximation.
514 : r = sqrt(-log(betaInc**2))
515 : betaIncNew = r - (2.30753_RKC + 0.27061_RKC * r) / (1._RKC + (0.99229_RKC + 0.04481_RKC * r) * r)
516 : if (1._RKC < alpha .and. 1._RKC < beta) then
517 : r = (betaIncNew * betaIncNew - 3._RKC) * ONE_SIXTH
518 : s = 1._RKC / (alpha + alpha - 1._RKC)
519 : t = 1._RKC / (beta + beta - 1._RKC)
520 : h = 2._RKC / (s + t)
521 : w = betaIncNew * sqrt(h + r) / h - (t - s) * (r + 5._RKC / 6._RKC - 2._RKC / (3._RKC * h))
522 : betaInv = alpha / (alpha + beta * exp(w + w))
523 : else
524 : r = beta + beta
525 : t = 1._RKC / (9._RKC * beta)
526 : t = r * (1._RKC - t + betaIncNew * sqrt(t))**3
527 : if (t <= 0._RKC) then
528 : betaInv = 1._RKC - exp((log((1._RKC - betaInc) * beta) + logFuncBeta) / beta)
529 : else
530 : t = (4._RKC * alpha + r - 2._RKC) / t
531 : if (t <= 1._RKC) then
532 : betaInv = exp((log(betaInc * alpha) + logFuncBeta) / alpha)
533 : else
534 : betaInv = 1._RKC - 2._RKC / (t + 1._RKC)
535 : end if
536 : end if
537 : end if
538 : abstol = epsilon(0._RKC)
539 : if (betaInv < 0._RKC) betaInv = sqrt(abstol)
540 : if (1._RKC < betaInv) betaInv = 1._RKC - sqrt(abstol)
541 : lb = 0._RKC
542 : ub = 1._RKC
543 : lf = betaInc
544 : uf = betaInc - 1._RKC
545 : if (abstol < ub - lb) then
546 : call setRoot(newton, getFunc, betaInv, getFuncDiff, lb, ub, lf, uf, abstol, neval)
547 : if (neval < 0_IK) call setRoot(bisection, getFunc, betaInv, lb, ub, lf, uf, abstol, neval)
548 : if (neval < 0_IK) info = 1_IK
549 : else
550 : if (abs(lf) < abs(uf)) then
551 : betaInv = lb
552 : else
553 : betaInv = ub
554 : end if
555 : end if
556 : if (mirrored) then
557 : if (signed) then
558 : betaInv = -betaInv
559 : else
560 : betaInv = 1._RKC - betaInv
561 : end if
562 : end if
563 :
564 : elseif (0._RKC == betaInc .or. betaInc == 1._RKC) then
565 : info = 0_IK
566 : betaInv = betaInc
567 : else
568 : info = 1_IK
569 : CHECK_ASSERTION(__LINE__, info == 0_IK, SK_"@setBetaInc(): The condition `0. <= betaInc .and. betaInc <= 1.` must hold. betaInc = "//getStr(betaInc)) ! fpp
570 : end if
571 : contains
572 : PURE function getFunc(x) result(func)
573 : real(RKC) , intent(in) :: x ! betaInv guess.
574 : real(RKC) :: func
575 : integer(IK) :: info
576 : call setBetaInc(func, x, alpha, beta, logFuncBeta, signed, info)
577 : if (info /= 0_IK) error stop ! LCOV_EXCL_LINE
578 : func = betaInc - func
579 : end function
580 : PURE function getFuncDiff(x, order) result(betaPDF)
581 : integer(IK) , intent(in) :: order
582 : real(RKC) , intent(in) :: x
583 : real(RKC) , parameter :: SQRT_HUGE = sqrt(huge(x))
584 : real(RKC) :: betaPDF
585 : if(order == 0_IK) then
586 : betaPDF = getFunc(x)
587 : else
588 : if (x /= 0._RKC .and. x /= 1._RKC) then
589 : call setBetaLogPDF(betaPDF, x, alpha, beta)
590 : betaPDF = exp(betaPDF)
591 : elseif ((x == 0._RKC .and. alpha < 1._RKC) .or. (x == 1._RKC .and. beta < 1._RKC)) then
592 : betaPDF = SQRT_HUGE
593 : else
594 : betaPDF = 0._RKC
595 : end if
596 : end if
597 : end function
598 :
599 : #else
600 : !%%%%%%%%%%%%%%%%%%%%%%%%
601 : #error "Unrecognized interface."
602 : !%%%%%%%%%%%%%%%%%%%%%%%%
603 : #endif
|
