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_distBeta](@ref pm_distBeta).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : !%%%%%%%%%%%%%%%%%
28 : #if getBetaPDF_ENABLED
29 : !%%%%%%%%%%%%%%%%%
30 :
31 63520 : if (x == 0._RKC) then
32 3 : if (alpha < 1._RKC) then
33 0 : pdf = +huge(pdf)
34 : else
35 3 : pdf = 0._RKC
36 : end if
37 63517 : elseif (x == 1._RKC) then
38 3 : if (beta < 1._RKC) then
39 0 : pdf = +huge(pdf)
40 : else
41 3 : pdf = 0._RKC
42 : end if
43 : else
44 63514 : call setBetaLogPDF(pdf, x, alpha, beta)
45 63514 : pdf = exp(pdf)
46 : end if
47 :
48 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49 : #elif getBetaLogPDF_ENABLED || setBetaLogPDF_ENABLED
50 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51 :
52 72689 : CHECK_ASSERTION(__LINE__, 0._RKC < x .and. x < 1._RKC, SK_"@setBetaLogPDF(): The condition `0. < x .and. x < 1.` must hold. x = "//getStr(x))
53 72689 : CHECK_ASSERTION(__LINE__, 0._RKC < alpha, SK_"@setBetaLogPDF(): The condition `0. < alpha` must hold. alpha = "//getStr(alpha))
54 72689 : CHECK_ASSERTION(__LINE__, 0._RKC < beta, SK_"@setBetaLogPDF(): The condition `0. < beta` must hold. alpha = "//getStr(beta))
55 : #if LOGBETA_ENABLED
56 6 : CHECK_ASSERTION(__LINE__, abs(logBeta - getLogBeta(alpha, beta)) < sqrt(epsilon(0._RKC)), \
57 : SK_"@setBetaLogPDF(): The condition `abs(logBeta - getLogBeta(alpha, beta)) < sqrt(epsilon(0._RKC)` must hold. logBeta, getLogBeta(alpha, beta) = "\
58 : //getStr([logBeta, getLogBeta(alpha, beta)]))
59 : #endif
60 : logPDF = (alpha - 1._RKC) * log(x) + (beta - 1._RKC) * log(1._RKC - x) - & ! LCOV_EXCL_LINE
61 : #if DEFAULT_ENABLED
62 72687 : getLogBeta(alpha, beta)
63 : #elif LOGBETA_ENABLED
64 2 : logBeta
65 : #else
66 : #error "Unrecognized interface."
67 : #endif
68 :
69 : !%%%%%%%%%%%%%%%%%
70 : #elif getBetaCDF_ENABLED
71 : !%%%%%%%%%%%%%%%%%
72 :
73 4009 : cdf = getBetaInc(x, alpha, beta, signed)
74 :
75 : !%%%%%%%%%%%%%%%%%
76 : #elif setBetaCDF_ENABLED
77 : !%%%%%%%%%%%%%%%%%
78 :
79 4011 : call setBetaInc(cdf, x, alpha, beta, logFuncBeta, signed, info)
80 :
81 : !%%%%%%%%%%%%%%%%%%
82 : #elif setBetaRand_ENABLED
83 : !%%%%%%%%%%%%%%%%%%
84 :
85 : ! Set the URNG.
86 : #if RNGD_ENABLED
87 : #define RNG
88 : #elif RNGF_ENABLED || RNGX_ENABLED
89 : #define RNG rng,
90 : #else
91 : #error "Unrecognized interface."
92 : #endif
93 229 : CHECK_ASSERTION(__LINE__, 0._RKC < alpha, SK_"@setBetaRand(): The condition `0. < alpha` must hold. alpha = "//getStr(alpha))
94 229 : CHECK_ASSERTION(__LINE__, 0._RKC < beta, SK_"@setBetaRand(): The condition `0. < beta` must hold. alpha = "//getStr(beta))
95 229 : if (1._RKC < alpha .or. 1._RKC < beta) then ! Use the algorithm of Johnk.
96 222 : block
97 : #if D0_ENABLED
98 : real(RKC) :: temp
99 : #elif D1_ENABLED
100 3 : real(RKC) :: temp(size(rand, 1, IK))
101 : #else
102 : #error "Unrecognized interface."
103 : #endif
104 222 : call setGammaRand(RNG rand, alpha, 1._RKC)
105 222 : call setGammaRand(RNG temp, beta, 1._RKC)
106 15222 : rand = rand / (rand + temp)
107 : end block
108 : else
109 : block
110 : #if D1_ENABLED
111 : integer(IK) :: irand
112 : #endif
113 : real(RKC) :: temp, invShape(2), dumm
114 7 : invShape(1) = 1._RKC / alpha
115 4 : invShape(2) = 1._RKC / beta
116 : #if D1_ENABLED
117 : #define GET_RAND(i) rand(i)
118 5007 : do irand = 1, size(rand, 1, IK)
119 : #elif D0_ENABLED
120 : #define GET_RAND(i) rand
121 : #endif
122 3 : do
123 6344 : call setUnifRand(RNG GET_RAND(irand))
124 6344 : call setUnifRand(RNG temp)
125 6344 : GET_RAND(irand) = (1._RKC - GET_RAND(irand))**invShape(1)
126 6344 : temp = (1._RKC - temp)**invShape(2)
127 6344 : dumm = GET_RAND(irand) + temp
128 : ! Rejection and cycling happens only if any(unifrnd == 0._RKC),
129 : ! which is approximately 1 in 10^106 odds of occurrence.
130 6344 : if (1._RKC < dumm) cycle
131 5008 : if (0._RKC < dumm) then
132 5008 : GET_RAND(irand) = GET_RAND(irand) / dumm
133 : else
134 0 : GET_RAND(irand) = log(GET_RAND(irand))
135 0 : temp = log(temp)
136 0 : dumm = max(GET_RAND(irand), temp)
137 0 : GET_RAND(irand) = GET_RAND(irand) - dumm
138 0 : temp = temp - dumm
139 0 : GET_RAND(irand) = exp(GET_RAND(irand) - log(exp(GET_RAND(irand)) + exp(temp)))
140 : end if
141 1336 : exit
142 : end do
143 : #if D1_ENABLED
144 : end do
145 : #endif
146 : end block
147 : end if
148 :
149 : #else
150 : !%%%%%%%%%%%%%%%%%%%%%%%%
151 : #error "Unrecognized interface."
152 : !%%%%%%%%%%%%%%%%%%%%%%%%
153 : #endif
154 :
155 : #undef GET_RAND
156 : #undef RNG
|
