Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!!
4 : !!!! MIT License
5 : !!!!
6 : !!!! ParaMonte: plain powerful parallel Monte Carlo library.
7 : !!!!
8 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab
9 : !!!!
10 : !!!! This file is part of the ParaMonte library.
11 : !!!!
12 : !!!! Permission is hereby granted, free of charge, to any person obtaining a
13 : !!!! copy of this software and associated documentation files (the "Software"),
14 : !!!! to deal in the Software without restriction, including without limitation
15 : !!!! the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : !!!! and/or sell copies of the Software, and to permit persons to whom the
17 : !!!! Software is furnished to do so, subject to the following conditions:
18 : !!!!
19 : !!!! The above copyright notice and this permission notice shall be
20 : !!!! included in all copies or substantial portions of the Software.
21 : !!!!
22 : !!!! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
23 : !!!! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
24 : !!!! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
25 : !!!! IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
26 : !!!! DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
27 : !!!! OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
28 : !!!! OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
29 : !!!!
30 : !!!! ACKNOWLEDGMENT
31 : !!!!
32 : !!!! ParaMonte is an honor-ware and its currency is acknowledgment and citations.
33 : !!!! As per the ParaMonte library license agreement terms, if you use any parts of
34 : !!!! this library for any purposes, kindly acknowledge the use of ParaMonte in your
35 : !!!! work (education/research/industry/development/...) by citing the ParaMonte
36 : !!!! library as described on this page:
37 : !!!!
38 : !!!! https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
39 : !!!!
40 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
41 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
42 :
43 : !> \brief This module contains the classes and procedures for various statistical computations.
44 : !> \author Amir Shahmoradi
45 :
46 : module Statistics_mod
47 :
48 : use Constants_mod, only: RK, IK
49 :
50 : implicit none
51 :
52 : !logical, save :: paradramPrintEnabled = .false.
53 : !logical, save :: paradisePrintEnabled = .false.
54 :
55 : character(len=*), parameter :: MODULE_NAME = "@Statistics_mod"
56 :
57 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58 :
59 : interface getSNormCDF
60 : module procedure :: getSNormCDF_SPR, getSNormCDF_DPR
61 : end interface getSNormCDF
62 :
63 : interface getBetaCDF
64 : module procedure :: getBetaCDF_SPR, getBetaCDF_DPR
65 : end interface getBetaCDF
66 :
67 : interface getBetaContinuedFraction
68 : module procedure :: getBetaContinuedFraction_SPR, getBetaContinuedFraction_DPR
69 : end interface getBetaContinuedFraction
70 :
71 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
72 :
73 : interface getMean
74 : module procedure :: getMean_2D
75 : end interface getMean
76 :
77 : interface flatten
78 : module procedure :: flatten_2D
79 : end interface flatten
80 :
81 : interface getNormData
82 : module procedure :: getNormData_2D, normalizeWeightedData_2d
83 : end interface getNormData
84 :
85 : interface getVariance
86 : module procedure :: getVariance_1D, getVariance_2D
87 : end interface getVariance
88 :
89 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
90 :
91 : interface getLogProbLogn
92 : module procedure :: getLogProbLognSP, getLogProbLognMP
93 : end interface getLogProbLogn
94 :
95 : interface getLogProbLognorm
96 : module procedure :: getLogProbLognSP, getLogProbLognMP
97 : end interface getLogProbLognorm
98 :
99 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
100 :
101 : interface getLogProbNormSP
102 : module procedure :: getLogProbNormSP_RK, getLogProbNormSP_CK
103 : end interface getLogProbNormSP
104 :
105 : interface getLogProbNormMP
106 : module procedure :: getLogProbNormMP_RK, getLogProbNormMP_CK
107 : end interface getLogProbNormMP
108 :
109 : interface getLogProbMVNSP
110 : module procedure :: getLogProbMVNSP_RK, getLogProbMVNSP_CK
111 : end interface getLogProbMVNSP
112 :
113 : interface getLogProbMVNMP
114 : module procedure :: getLogProbMVNMP_RK, getLogProbMVNMP_CK
115 : end interface getLogProbMVNMP
116 :
117 : interface getLogProbNorm
118 : module procedure :: getLogProbNormSP_RK, getLogProbNormMP_RK
119 : module procedure :: getLogProbNormSP_CK, getLogProbNormMP_CK
120 : end interface getLogProbNorm
121 :
122 : interface getLogProbMVN
123 : module procedure :: getLogProbMVNSP_RK, getLogProbMVNMP_RK
124 : module procedure :: getLogProbMVNSP_CK, getLogProbMVNMP_CK
125 : end interface getLogProbMVN
126 :
127 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
128 :
129 : interface getLogProbMixNorm
130 : module procedure :: getLogProbMixNormSP_RK, getLogProbMixNormSP_CK
131 : module procedure :: getLogProbMixNormMP_RK, getLogProbMixNormMP_CK
132 : end interface getLogProbMixNorm
133 :
134 : interface getLogProbMixMVN
135 : module procedure :: getLogProbMixMVNSP_RK, getLogProbMixMVNSP_CK
136 : module procedure :: getLogProbMixMVNMP_RK, getLogProbMixMVNMP_CK
137 : end interface getLogProbMixMVN
138 :
139 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
140 :
141 : interface getMahalSqSP
142 : module procedure :: getMahalSqSP_RK, getMahalSqSP_CK
143 : end interface getMahalSqSP
144 :
145 : interface getMahalSqMP
146 : module procedure :: getMahalSqMP_RK, getMahalSqMP_CK
147 : end interface getMahalSqMP
148 :
149 : interface getMahalSq
150 : module procedure :: getMahalSqSP_RK, getMahalSqMP_RK
151 : module procedure :: getMahalSqSP_CK, getMahalSqMP_CK
152 : end interface getMahalSq
153 :
154 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
155 :
156 : contains
157 :
158 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
159 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
160 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
161 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
162 :
163 :
164 : !> \brief
165 : !> Return the square of Mahalanobis distance for a single point. The output is a scalar variable.
166 : !>
167 : !> \param[in] nd : The number of dimensions of the input `Point`.
168 : !> \param[in] MeanVec : The mean vector of the sample.
169 : !> \param[in] InvCovMat : The inverse covariance matrix of the sample.
170 : !> \param[in] Point : The `Point` whose distance from the sample is to computed.
171 : !>
172 : !> \return
173 : !> `mahalSq` : The Mahalanobis distance squared of the point from
174 : !> the sample characterized by the input `MeanVec` and `InvCovMat`.
175 12 : pure function getMahalSqSP_RK(nd,MeanVec,InvCovMat,Point) result(mahalSq)
176 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
177 : !DEC$ ATTRIBUTES DLLEXPORT :: getMahalSqSP_RK
178 : #endif
179 : use Constants_mod, only: IK, RK
180 : implicit none
181 : integer(IK), intent(in) :: nd
182 : real(RK) , intent(in) :: MeanVec(nd)
183 : real(RK) , intent(in) :: InvCovMat(nd,nd) ! Inverse of the covariance matrix
184 : real(RK) , intent(in) :: Point(nd) ! input data points
185 : real(RK) :: mahalSq
186 48 : real(RK) :: NormedPoint(nd)
187 48 : NormedPoint = Point - MeanVec
188 48 : mahalSq = dot_product( NormedPoint , matmul(InvCovMat,NormedPoint) )
189 12 : end function getMahalSqSP_RK
190 :
191 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
192 :
193 : !> \brief
194 : !> Return the square of Mahalanobis distances for an row-wise array of points.
195 : !>
196 : !> \param[in] nd : The number of dimensions of the input `Point` array.
197 : !> \param[in] np : The number of points in the input input `Point` array.
198 : !> \param[in] MeanVec : The mean vector of length `nd` of the sample.
199 : !> \param[in] InvCovMat : The inverse covariance matrix `(nd,nd)` of the sample.
200 : !> \param[in] Point : The `(nd,np)` array of points whose distances squared from the sample are to computed.
201 : !>
202 : !> \return
203 : !> `MahalSq` : A vector of length `np` containing the squares of the Mahalanobis distances
204 : !> of the input points from the sample characterized by the input `MeanVec` and `InvCovMat`.
205 : !>
206 : !> \warning
207 : !> For the sake of preserving the purity and computational efficiency of the function,
208 : !> if the computation fails at any stage, the first element of output will be returned negative.
209 60 : pure function getMahalSqMP_RK(nd,np,MeanVec,InvCovMat,Point) result(MahalSq)
210 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
211 : !DEC$ ATTRIBUTES DLLEXPORT :: getMahalSqMP_RK
212 : #endif
213 12 : use Constants_mod, only: IK, RK
214 : implicit none
215 : integer(IK), intent(in) :: nd,np
216 : real(RK), intent(in) :: MeanVec(nd)
217 : real(RK), intent(in) :: InvCovMat(nd,nd) ! Inverse of the covariance matrix
218 : real(RK), intent(in) :: Point(nd,np) ! input data points
219 : real(RK) :: MahalSq(np) ! function return
220 48 : real(RK) :: NormedPoint(nd)
221 : integer(IK) :: ip
222 36 : do ip = 1,np
223 96 : NormedPoint = Point(1:nd,ip) - MeanVec
224 96 : MahalSq(ip) = dot_product( NormedPoint , matmul(InvCovMat,NormedPoint) )
225 36 : if (MahalSq(ip)<0._RK) then
226 : ! LCOV_EXCL_START
227 : MahalSq(1) = -1._RK
228 : return
229 : end if
230 : ! LCOV_EXCL_STOP
231 : end do
232 24 : end function getMahalSqMP_RK
233 :
234 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
235 :
236 : !> \brief
237 : !> Return the square of Mahalanobis distance for a single complex point. The output is a scalar variable.
238 : !>
239 : !> \param[in] nd : The number of dimensions of the input `Point`.
240 : !> \param[in] MeanVec : The mean vector of the sample.
241 : !> \param[in] InvCovMat : The inverse covariance matrix of the sample.
242 : !> \param[in] Point : The `Point` whose distance from the sample is to computed.
243 : !>
244 : !> \return
245 : !> `mahalSq` : The Mahalanobis distance squared of the point from
246 : !> the sample characterized by the input `MeanVec` and `InvCovMat`.
247 12 : pure function getMahalSqSP_CK(nd,MeanVec,InvCovMat,Point) result(mahalSq)
248 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
249 : !DEC$ ATTRIBUTES DLLEXPORT :: getMahalSqSP_CK
250 : #endif
251 12 : use Constants_mod, only: IK, RK, CK
252 : implicit none
253 : integer(IK), intent(in) :: nd
254 : complex(CK), intent(in) :: MeanVec(nd)
255 : complex(CK), intent(in) :: InvCovMat(nd,nd) ! Inverse of the covariance matrix
256 : complex(CK), intent(in) :: Point(nd) ! input data points
257 : complex(CK) :: mahalSq ! function return
258 84 : mahalSq = sum( (Point-MeanVec) * matmul(InvCovMat,Point-MeanVec) )
259 12 : end function getMahalSqSP_CK
260 :
261 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
262 :
263 : !> \brief
264 : !> Return the square of Mahalanobis distances for an row-wise array of complex-valued points.
265 : !>
266 : !> \param[in] nd : The number of dimensions of the input `Point` array.
267 : !> \param[in] np : The number of points in the input input `Point` array.
268 : !> \param[in] MeanVec : The mean vector of length `nd` of the sample.
269 : !> \param[in] InvCovMat : The inverse covariance matrix `(nd,nd)` of the sample.
270 : !> \param[in] Point : The `(nd,np)` array of points whose distances squared from the sample are to computed.
271 : !>
272 : !> \return
273 : !> `MahalSq` : A vector of length `np` containing the squares of the Mahalanobis distances
274 : !> of the input points from the sample characterized by the input `MeanVec` and `InvCovMat`.
275 : !>
276 : !> \warning
277 : !> For the sake of preserving the purity and computational efficiency of the function,
278 : !> if the computation fails at any stage, the first element of output will be returned negative.
279 60 : pure function getMahalSqMP_CK(nd,np,MeanVec,InvCovMat,Point) result(MahalSq)
280 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
281 : !DEC$ ATTRIBUTES DLLEXPORT :: getMahalSqMP_CK
282 : #endif
283 12 : use Constants_mod, only: IK, RK, CK
284 : implicit none
285 : integer(IK), intent(in) :: nd,np
286 : complex(CK), intent(in) :: MeanVec(nd)
287 : complex(CK), intent(in) :: InvCovMat(nd,nd) ! Inverse of the covariance matrix
288 : complex(CK), intent(in) :: Point(nd,np) ! input data points
289 : complex(CK) :: MahalSq(np) ! function return
290 : integer(IK) :: ip
291 36 : do ip = 1,np
292 48 : MahalSq(ip) = sum( (Point(1:nd,ip)-MeanVec) * &
293 192 : matmul(InvCovMat,Point(1:nd,ip)-MeanVec) )
294 36 : if (real(MahalSq(ip))<0._RK) then
295 : ! LCOV_EXCL_START
296 : MahalSq(1) = (-1._RK, -1._RK)
297 : return
298 : end if
299 : ! LCOV_EXCL_STOP
300 : end do
301 24 : end function getMahalSqMP_CK
302 :
303 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
304 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
305 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
306 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
307 :
308 9 : pure function getLogProbNormSP_RK(mean,inverseVariance,logSqrtInverseVariance,point) result(logProbNorm)
309 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
310 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbNormSP_RK
311 : #endif
312 12 : use Constants_mod, only: RK, LOGINVSQRT2PI
313 : implicit none
314 : real(RK), intent(in) :: mean,inverseVariance,logSqrtInverseVariance,point
315 : real(RK) :: logProbNorm
316 9 : logProbNorm = LOGINVSQRT2PI + logSqrtInverseVariance - 0.5_RK * inverseVariance * (point-mean)**2
317 18 : end function getLogProbNormSP_RK
318 :
319 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
320 :
321 45 : pure function getLogProbNormMP_RK(np,mean,inverseVariance,logSqrtInverseVariance,Point) result(logProbNorm)
322 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
323 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbNormMP_RK
324 : #endif
325 9 : use Constants_mod, only: LOGINVSQRT2PI
326 : implicit none
327 : integer(IK), intent(in) :: np
328 : real(RK) , intent(in) :: mean,inverseVariance,logSqrtInverseVariance,Point(np)
329 : real(RK) :: logProbNorm(np)
330 27 : logProbNorm = LOGINVSQRT2PI + logSqrtInverseVariance - 0.5_RK * inverseVariance * (Point-mean)**2
331 9 : end function getLogProbNormMP_RK
332 :
333 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
334 :
335 : ! NOTE: if MahalSq computation fails, output probability will be returned as NullVal%RK from module Constant_mod.
336 9 : pure function getLogProbMVNSP_RK(nd,MeanVec,InvCovMat,logSqrtDetInvCovMat,Point) result(logProbNorm)
337 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
338 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMVNSP_RK
339 : #endif
340 9 : use Constants_mod, only: LOGINVSQRT2PI, NullVal
341 : implicit none
342 : integer(IK), intent(in) :: nd
343 : real(RK) , intent(in) :: MeanVec(nd)
344 : real(RK) , intent(in) :: InvCovMat(nd,nd)
345 : real(RK) , intent(in) :: logSqrtDetInvCovMat
346 : real(RK) , intent(in) :: Point(nd)
347 9 : real(RK) :: logProbNorm, dummy
348 9 : dummy = getMahalSqSP_RK(nd,MeanVec,InvCovMat,Point)
349 9 : if (dummy<0._RK) then
350 0 : logProbNorm = NullVal%RK
351 : else
352 9 : logProbNorm = nd*LOGINVSQRT2PI + logSqrtDetInvCovMat - 0.5_RK * dummy
353 : end if
354 9 : end function getLogProbMVNSP_RK
355 :
356 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
357 :
358 : ! NOTE: if MahalSq computation fails, output probability will be returned as NullVal%RK from module Constant_mod.
359 45 : pure function getLogProbMVNMP_RK(nd,np,MeanVec,InvCovMat,logSqrtDetInvCovMat,Point) result(logProbNorm)
360 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
361 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMVNMP_RK
362 : #endif
363 9 : use Constants_mod, only: LOGINVSQRT2PI, NullVal
364 : implicit none
365 : integer(IK), intent(in) :: nd,np
366 : real(RK) , intent(in) :: MeanVec(nd)
367 : real(RK) , intent(in) :: InvCovMat(nd,nd)
368 : real(RK) , intent(in) :: logSqrtDetInvCovMat
369 : real(RK) , intent(in) :: Point(nd,np)
370 36 : real(RK) :: logProbNorm(np), Dummy(np)
371 9 : Dummy = getMahalSqMP_RK(nd,np,MeanVec,InvCovMat,Point)
372 9 : if (Dummy(1)<0._RK) then
373 0 : logProbNorm = NullVal%RK
374 : else
375 27 : logProbNorm = nd*LOGINVSQRT2PI + logSqrtDetInvCovMat - 0.5_RK * Dummy
376 : end if
377 9 : end function getLogProbMVNMP_RK
378 :
379 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
380 :
381 9 : function getLogProbNormSP_CK(mean,inverseVariance,logSqrtInverseVariance,point) result(logProbNorm)
382 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
383 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbNormSP_CK
384 : #endif
385 9 : use Constants_mod, only: RK, CK, LOGINVSQRT2PI
386 : implicit none
387 : complex(CK), intent(in) :: mean,inverseVariance,logSqrtInverseVariance,point
388 : complex(CK) :: logProbNorm
389 9 : logProbNorm = LOGINVSQRT2PI + logSqrtInverseVariance - 0.5_RK * inverseVariance * (point-mean)**2
390 18 : end function getLogProbNormSP_CK
391 :
392 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
393 :
394 45 : function getLogProbNormMP_CK(np,mean,inverseVariance,logSqrtInverseVariance,Point) result(logProbNorm)
395 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
396 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbNormMP_CK
397 : #endif
398 9 : use Constants_mod, only: IK, RK, CK, LOGINVSQRT2PI
399 : implicit none
400 : integer(IK), intent(in) :: np
401 : complex(CK) , intent(in) :: mean,inverseVariance,logSqrtInverseVariance,Point(np)
402 : complex(CK) :: logProbNorm(np)
403 27 : logProbNorm = LOGINVSQRT2PI + logSqrtInverseVariance - 0.5_RK * inverseVariance * (Point-mean)**2
404 9 : end function getLogProbNormMP_CK
405 :
406 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
407 :
408 9 : function getLogProbMVNSP_CK(nd,MeanVec,InvCovMat,logSqrtDetInvCovMat,Point) result(logProbMVN)
409 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
410 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMVNSP_CK
411 : #endif
412 9 : use Constants_mod, only: IK, RK, CK, LOGINVSQRT2PI, NullVal
413 : implicit none
414 : integer(IK), intent(in) :: nd
415 : complex(CK), intent(in) :: MeanVec(nd)
416 : complex(CK), intent(in) :: InvCovMat(nd,nd)
417 : complex(CK), intent(in) :: logSqrtDetInvCovMat
418 : complex(CK), intent(in) :: Point(nd)
419 : complex(CK) :: logProbMVN, dummy
420 9 : dummy = getMahalSqSP(nd,MeanVec,InvCovMat,Point)
421 9 : if (real(dummy)<0._RK) then
422 0 : logProbMVN = NullVal%RK
423 : else
424 9 : logProbMVN = nd*LOGINVSQRT2PI + logSqrtDetInvCovMat - 0.5_RK * dummy
425 : end if
426 9 : end function getLogProbMVNSP_CK
427 :
428 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
429 :
430 45 : function getLogProbMVNMP_CK(nd,np,MeanVec,InvCovMat,logSqrtDetInvCovMat,Point) result(logProbMVN)
431 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
432 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMVNMP_CK
433 : #endif
434 9 : use Constants_mod, only: IK, RK, CK, LOGINVSQRT2PI, NullVal
435 : implicit none
436 : integer(IK), intent(in) :: nd,np
437 : complex(CK), intent(in) :: MeanVec(nd)
438 : complex(CK), intent(in) :: InvCovMat(nd,nd)
439 : complex(CK), intent(in) :: logSqrtDetInvCovMat
440 : complex(CK), intent(in) :: Point(nd,np)
441 36 : complex(CK) :: logProbMVN(np), Dummy(np)
442 9 : Dummy = getMahalSqMP(nd,np,MeanVec,InvCovMat,Point)
443 9 : if (real(Dummy(1))<0._RK) then
444 0 : logProbMVN = NullVal%RK
445 : else
446 27 : logProbMVN = nd*LOGINVSQRT2PI + logSqrtDetInvCovMat - 0.5_RK * Dummy
447 : end if
448 9 : end function getLogProbMVNMP_CK
449 :
450 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
451 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
452 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
453 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
454 :
455 : ! SDSP stands for Single-Dimensional Gaussian mixture with Single Point input
456 : ! For a proper probability normalization, the sum of the amplitudes must equal one.
457 3 : function getLogProbMixNormSP_RK(nmode,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,point) result(logProbMixNorm)
458 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
459 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixNormSP_RK
460 : #endif
461 9 : use Constants_mod, only: IK, RK, LOGTINY_RK
462 : implicit none
463 : integer(IK), intent(in) :: nmode
464 : real(RK) , intent(in) :: LogAmplitude(nmode),MeanVec(nmode)
465 : real(RK) , intent(in) :: InvCovMat(nmode),LogSqrtDetInvCovMat(nmode)
466 : real(RK) , intent(in) :: point
467 : real(RK) :: logProbMixNorm
468 9 : real(RK) :: normFac,LogProb(nmode)
469 : integer(IK) :: imode
470 9 : do imode = 1, nmode
471 9 : LogProb(imode) = LogAmplitude(imode) + getLogProbNormSP_RK(MeanVec(imode),InvCovMat(imode),LogSqrtDetInvCovMat(imode),point)
472 : end do
473 12 : normFac = maxval(LogProb)
474 9 : LogProb = LogProb - normFac
475 33 : where(LogProb<LOGTINY_RK)
476 3 : LogProb = 0._RK
477 : elsewhere
478 3 : LogProb = exp(LogProb)
479 : end where
480 9 : logProbMixNorm = normFac + log(sum(LogProb))
481 3 : end function getLogProbMixNormSP_RK
482 :
483 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
484 :
485 15 : function getLogProbMixNormMP_RK(nmode,np,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,Point) result(logProbMixNorm)
486 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
487 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixNormMP_RK
488 : #endif
489 3 : use Constants_mod, only: IK, RK, LOGTINY_RK
490 : implicit none
491 : integer(IK), intent(in) :: nmode,np
492 : real(RK) , intent(in) :: LogAmplitude(nmode),MeanVec(nmode)
493 : real(RK) , intent(in) :: InvCovMat(nmode),LogSqrtDetInvCovMat(nmode)
494 : real(RK) , intent(in) :: Point(np)
495 : real(RK) :: logProbMixNorm(np)
496 30 : real(RK) :: NormFac(np),LogProb(nmode,np)
497 : integer(IK) :: imode, ip
498 9 : do imode = 1, nmode
499 21 : LogProb(imode,1:np) = LogAmplitude(imode) + getLogProbNormMP_RK(np,MeanVec(imode),InvCovMat(imode),LogSqrtDetInvCovMat(imode),Point)
500 : end do
501 3 : NormFac = maxval(LogProb,dim=1)
502 9 : do ip = 1,np
503 18 : LogProb(1:nmode,ip) = LogProb(1:nmode,ip) - NormFac(ip)
504 18 : do imode = 1,nmode
505 18 : if ( LogProb(imode,ip) < LOGTINY_RK ) then
506 0 : LogProb(imode,ip) = 0._RK
507 : else
508 12 : LogProb(imode,ip) = exp( LogProb(imode,ip) )
509 : end if
510 : end do
511 21 : logProbMixNorm(ip) = NormFac(ip) + log(sum(LogProb(1:nmode,ip)))
512 : end do
513 3 : end function getLogProbMixNormMP_RK
514 :
515 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
516 :
517 3 : function getLogProbMixMVNSP_RK(nmode,nd,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,Point) result(logProbMixMVN)
518 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
519 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixMVNSP_RK
520 : #endif
521 3 : use Constants_mod, only: IK, RK, LOGTINY_RK
522 : implicit none
523 : integer(IK), intent(in) :: nmode,nd
524 : real(RK) , intent(in) :: LogAmplitude(nmode), MeanVec(nd,nmode)
525 : real(RK) , intent(in) :: InvCovMat(nd,nd,nmode), LogSqrtDetInvCovMat(nmode)
526 : real(RK) , intent(in) :: Point(nd)
527 : real(RK) :: logProbMixMVN
528 9 : real(RK) :: normFac,LogProb(nmode)
529 : integer(IK) :: imode
530 9 : do imode = 1, nmode
531 9 : LogProb(imode) = LogAmplitude(imode) + getLogProbMVNSP_RK(nd,MeanVec(1:nd,imode),InvCovMat(1:nd,1:nd,imode),LogSqrtDetInvCovMat(imode),Point)
532 : end do
533 12 : normFac = maxval(LogProb)
534 9 : LogProb = LogProb - normFac
535 33 : where(LogProb<LOGTINY_RK)
536 3 : LogProb = 0._RK
537 : elsewhere
538 3 : LogProb = exp(LogProb)
539 : end where
540 9 : logProbMixMVN = normFac + log(sum(LogProb))
541 3 : end function getLogProbMixMVNSP_RK
542 :
543 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
544 :
545 15 : function getLogProbMixMVNMP_RK(nmode,nd,np,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,Point) result(logProbMixMVN)
546 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
547 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixMVNMP_RK
548 : #endif
549 3 : use Constants_mod, only: IK, RK, LOGTINY_RK
550 : implicit none
551 : integer(IK), intent(in) :: nmode,nd,np
552 : real(RK) , intent(in) :: LogAmplitude(nmode),MeanVec(nd,nmode)
553 : real(RK) , intent(in) :: InvCovMat(nd,nd,nmode), LogSqrtDetInvCovMat(nmode)
554 : real(RK) , intent(in) :: Point(nd,np)
555 : real(RK) :: logProbMixMVN(np)
556 30 : real(RK) :: NormFac(np),LogProb(nmode,np)
557 : integer(IK) :: imode, ip
558 9 : do imode = 1, nmode
559 12 : LogProb(imode,1:np) = LogAmplitude(imode) + &
560 33 : getLogProbMVNMP_RK(nd,np,MeanVec(1:nd,imode),InvCovMat(1:nd,1:nd,imode),LogSqrtDetInvCovMat(imode),Point)
561 : end do
562 3 : NormFac = maxval(LogProb,dim=1)
563 9 : do ip = 1,np
564 18 : LogProb(1:nmode,ip) = LogProb(1:nmode,ip) - NormFac(ip)
565 18 : do imode = 1,nmode
566 18 : if ( LogProb(imode,ip)<LOGTINY_RK ) then
567 6 : LogProb(imode,ip) = 0._RK
568 : else
569 6 : LogProb(imode,ip) = exp( LogProb(imode,ip) )
570 : end if
571 : end do
572 21 : logProbMixMVN(ip) = NormFac(ip) + log(sum(LogProb(1:nmode,ip)))
573 : end do
574 3 : end function getLogProbMixMVNMP_RK
575 :
576 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
577 :
578 : ! SDSP stands for 1-dimensional Gaussian mixture with scalar input point
579 3 : function getLogProbMixNormSP_CK(nmode,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,point) result(logProbMixNorm)
580 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
581 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixNormSP_CK
582 : #endif
583 3 : use Constants_mod, only: IK, RK, CK, LOGTINY_RK
584 : implicit none
585 : integer(IK), intent(in) :: nmode
586 : complex(CK), intent(in) :: LogAmplitude(nmode),MeanVec(nmode)
587 : complex(CK), intent(in) :: InvCovMat(nmode),LogSqrtDetInvCovMat(nmode)
588 : complex(CK), intent(in) :: point
589 : complex(CK) :: logProbMixNorm
590 9 : complex(CK) :: normFac,LogProb(nmode)
591 : integer(IK) :: imode
592 9 : do imode = 1, nmode
593 9 : LogProb(imode) = LogAmplitude(imode) + getLogProbNorm(MeanVec(imode),InvCovMat(imode),LogSqrtDetInvCovMat(imode),point)
594 : end do
595 12 : normFac = maxval(real(LogProb))
596 9 : LogProb = LogProb - normFac
597 33 : where(real(LogProb)<LOGTINY_RK)
598 3 : LogProb = 0._RK
599 : elsewhere
600 3 : LogProb = exp(LogProb)
601 : end where
602 9 : logProbMixNorm = normFac + log(sum(LogProb))
603 3 : end function getLogProbMixNormSP_CK
604 :
605 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
606 :
607 15 : function getLogProbMixNormMP_CK(nmode,np,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,Point) result(logProbMixNorm)
608 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
609 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixNormMP_CK
610 : #endif
611 3 : use Constants_mod, only: IK, RK, CK, LOGTINY_RK
612 : implicit none
613 : integer(IK), intent(in) :: nmode,np
614 : complex(CK), intent(in) :: LogAmplitude(nmode),MeanVec(nmode)
615 : complex(CK), intent(in) :: InvCovMat(nmode),LogSqrtDetInvCovMat(nmode)
616 : complex(CK), intent(in) :: Point(np)
617 : complex(CK) :: logProbMixNorm(np)
618 30 : complex(CK) :: normFac(np),LogProb(nmode,np)
619 : integer(IK) :: imode, ip
620 9 : do imode = 1, nmode
621 21 : LogProb(imode,1:np) = LogAmplitude(imode) + getLogProbNorm(np,MeanVec(imode),InvCovMat(imode),LogSqrtDetInvCovMat(imode),Point)
622 : end do
623 27 : normFac = maxval(real(LogProb),dim=1)
624 9 : do ip = 1,np
625 18 : LogProb(1:nmode,ip) = LogProb(1:nmode,ip) - normFac(ip)
626 18 : do imode = 1,nmode
627 18 : if ( real(LogProb(imode,ip)) < LOGTINY_RK ) then
628 0 : LogProb(imode,ip) = 0._RK
629 : else
630 12 : LogProb(imode,ip) = exp( LogProb(imode,ip) )
631 : end if
632 : end do
633 21 : logProbMixNorm(ip) = normFac(ip) + log(sum(LogProb(1:nmode,ip)))
634 : end do
635 3 : end function getLogProbMixNormMP_CK
636 :
637 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
638 :
639 3 : function getLogProbMixMVNSP_CK(nmode,nd,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,Point) result(logProbMixMVN)
640 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
641 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixMVNSP_CK
642 : #endif
643 3 : use Constants_mod, only: IK, RK, CK, LOGTINY_RK
644 : implicit none
645 : integer(IK), intent(in) :: nmode,nd
646 : complex(CK), intent(in) :: LogAmplitude(nmode), MeanVec(nd,nmode)
647 : complex(CK), intent(in) :: InvCovMat(nd,nd,nmode), LogSqrtDetInvCovMat(nmode)
648 : complex(CK), intent(in) :: Point(nd)
649 : complex(CK) :: logProbMixMVN
650 9 : complex(CK) :: normFac,LogProb(nmode)
651 : integer(IK) :: imode
652 9 : do imode = 1, nmode
653 9 : LogProb(imode) = LogAmplitude(imode) + getLogProbMVN(nd,MeanVec(1:nd,imode),InvCovMat(1:nd,1:nd,imode),LogSqrtDetInvCovMat(imode),Point)
654 : end do
655 12 : normFac = maxval(real(LogProb))
656 9 : LogProb = LogProb - normFac
657 33 : where(real(LogProb)<LOGTINY_RK)
658 3 : LogProb = 0._RK
659 : elsewhere
660 3 : LogProb = exp(LogProb)
661 : end where
662 9 : logProbMixMVN = normFac + log(sum(LogProb))
663 3 : end function getLogProbMixMVNSP_CK
664 :
665 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
666 :
667 15 : function getLogProbMixMVNMP_CK(nmode,nd,np,LogAmplitude,MeanVec,InvCovMat,LogSqrtDetInvCovMat,Point) result(logProbMixMVN)
668 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
669 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMixMVNMP_CK
670 : #endif
671 3 : use Constants_mod, only: IK, RK, CK, LOGTINY_RK
672 : implicit none
673 : integer(IK), intent(in) :: nmode,nd,np
674 : complex(CK), intent(in) :: LogAmplitude(nmode),MeanVec(nd,nmode)
675 : complex(CK), intent(in) :: InvCovMat(nd,nd,nmode), LogSqrtDetInvCovMat(nmode)
676 : complex(CK), intent(in) :: Point(nd,np)
677 : complex(CK) :: logProbMixMVN(np)
678 30 : complex(CK) :: normFac(np),LogProb(nmode,np)
679 : integer(IK) :: imode, ip
680 9 : do imode = 1, nmode
681 21 : LogProb(imode,1:np) = LogAmplitude(imode) + getLogProbMVN(nd,np,MeanVec(1:nd,imode),InvCovMat(1:nd,1:nd,imode),LogSqrtDetInvCovMat(imode),Point)
682 : end do
683 27 : normFac = maxval(real(LogProb),dim=1)
684 9 : do ip = 1,np
685 18 : LogProb(1:nmode,ip) = LogProb(1:nmode,ip) - normFac(ip)
686 18 : do imode = 1,nmode
687 18 : if ( real(LogProb(imode,ip))<LOGTINY_RK ) then
688 6 : LogProb(imode,ip) = 0._RK
689 : else
690 6 : LogProb(imode,ip) = exp( LogProb(imode,ip) )
691 : end if
692 : end do
693 21 : logProbMixMVN(ip) = normFac(ip) + log(sum(LogProb(1:nmode,ip)))
694 : end do
695 3 : end function getLogProbMixMVNMP_CK
696 :
697 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
698 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
699 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
700 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
701 :
702 : !include "Statistics_mod@MahalSq_RK.inc.f90"
703 : !include "Statistics_mod@MahalSq_CK.inc.f90"
704 : !include "Statistics_mod@LogProbGaus_RK.inc.f90"
705 : !include "Statistics_mod@LogProbGaus_CK.inc.f90"
706 : !include "Statistics_mod@LogProbGausMix_RK.inc.f90"
707 : !include "Statistics_mod@LogProbGausMix_CK.inc.f90"
708 :
709 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
710 :
711 : !> \brief
712 : !> Flatten the input `Point` array such that each element of the output
713 : !> `FlattenedPoint` array has the same unity weight as elements in the array.
714 : !>
715 : !> \param[in] nd : The number of dimensions of the input sample.
716 : !> \param[in] np : The number of points in the sample.
717 : !> \param[in] Point : The array of shape `(nd,np)` containing the sample.
718 : !> \param[in] Weight : The vector of length `np` containing the weights of points in the sample.
719 : !> The values of elements of Weight are allowed to be negative, in which case,
720 : !> the corresponding elements will be excluded from the output `FlattenedPoint`.
721 : !>
722 : !> \warning
723 : !> Note the shape of the input argument `Point(nd,np)`.
724 : !>
725 : !> \return
726 : !> `FlattenedPoint` : The flattened array whose elements all have the same weight.
727 11 : pure function flatten_2D(nd,np,Point,Weight) result(FlattenedPoint)
728 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
729 : !DEC$ ATTRIBUTES DLLEXPORT :: flatten_2D
730 : #endif
731 : ! the implementation for one-dimension is very concise and nice: mean = sum(Weight*Point) / sum(Weight)
732 : implicit none
733 : integer(IK), intent(in) :: np,nd ! np: number of observations, nd: number of variables for each observation
734 : real(RK) , intent(in) :: Point(nd,np) ! Point is the data matrix
735 : integer(IK), intent(in) :: Weight(np) ! sample weight
736 : integer(IK) :: ip, iweight, sumWeight, counter
737 : real(RK), allocatable :: FlattenedPoint(:,:)
738 11 : sumWeight = 0_IK
739 103 : do ip = 1, np
740 103 : if (Weight(ip)>0_IK) sumWeight = sumWeight + Weight(ip)
741 : end do
742 11 : allocate(FlattenedPoint(nd,sumWeight))
743 11 : counter = 0_IK
744 103 : do ip = 1, np
745 192 : do iweight = 1, Weight(ip)
746 89 : counter = counter + 1_IK
747 439 : FlattenedPoint(1:nd,counter) = Point(1:nd,ip)
748 : end do
749 : end do
750 3 : end function flatten_2D
751 :
752 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
753 :
754 : !> \brief
755 : !> Return the mean of a sample of multidimensional points.
756 : !>
757 : !> \param[in] nd : The number of dimensions of the input sample.
758 : !> \param[in] np : The number of points in the sample.
759 : !> \param[in] Point : The array of shape `(nd,np)` containing the sample.
760 : !> \param[in] Weight : The vector of length `np` containing the weights of points in the sample (**optional**, default = vector of ones).
761 : !>
762 : !> \warning
763 : !> Note the shape of the input argument `Point(nd,np)`.
764 : !>
765 : !> \return
766 : !> `Mean` : The output mean vector of length `nd`.
767 80 : pure function getMean_2D(nd,np,Point,Weight) result(Mean)
768 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
769 : !DEC$ ATTRIBUTES DLLEXPORT :: getMean_2D
770 : #endif
771 : ! the implementation for one-dimension is very concise and nice: mean = sum(Weight*Point) / sum(Weight)
772 : implicit none
773 : integer(IK), intent(in) :: np,nd ! np: number of observations, nd: number of variables for each observation
774 : real(RK) , intent(in) :: Point(nd,np) ! Point is the data matrix
775 : integer(IK), intent(in), optional :: Weight(np) ! sample weight
776 : real(RK) :: Mean(nd) ! output mean vector
777 : integer(IK) :: ip, sumWeight
778 50 : Mean = 0._RK
779 15 : if (present(Weight)) then
780 12 : sumWeight = 0_IK
781 51 : do ip = 1,np
782 39 : sumWeight = sumWeight + Weight(ip)
783 143 : Mean = Mean + Weight(ip) * Point(1:nd,ip)
784 : end do
785 38 : Mean = Mean / sumWeight
786 : else
787 18 : do ip = 1,np
788 63 : Mean = Mean + Point(1:nd,ip)
789 : end do
790 12 : Mean = Mean / real(np,kind=RK)
791 : end if
792 26 : end function getMean_2D
793 :
794 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
795 :
796 : !> \brief
797 : !> Return the normalized data with respect to the input mean vector of length `nd`.
798 : !>
799 : !> \param[in] nd : The number of dimensions of the input sample.
800 : !> \param[in] np : The number of points in the sample.
801 : !> \param[in] Mean : The mean vector of length `nd`.
802 : !> \param[in] Point : The array of shape `(nd,np)` containing the sample.
803 : !>
804 : !> \return
805 : !> `NormData` : The output normalized points array of shape `(np,nd)`.
806 : !>
807 : !> \remark
808 : !> Note the difference in the shape of the input `Point` vs. the output `NormData`.
809 63 : pure function getNormData_2D(nd,np,Mean,Point) result(NormData)
810 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
811 : !DEC$ ATTRIBUTES DLLEXPORT :: getNormData_2D
812 : #endif
813 : implicit none
814 : integer(IK), intent(in) :: np,nd ! np is the number of observations, nd is the number of parameters for each observation
815 : real(RK) , intent(in) :: Mean(nd) ! Mean vector
816 : real(RK) , intent(in) :: Point(nd,np) ! Point is the matrix of the data, CovMat contains the elements of the sample covariance matrix
817 : real(RK) :: NormData(np,nd)
818 : integer(IK) :: i
819 18 : do i = 1,np
820 63 : NormData(i,1:nd) = Point(1:nd,i) - Mean
821 : end do
822 18 : end function getNormData_2D
823 :
824 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
825 :
826 : !> \brief
827 : !> Return the normalized 2D data of size `(nd,np)` with respect to the mean of the data along the second dimension of length `np`.
828 : !>
829 : !> \param[in] nd : The length of the input `Data` matrix along the first dimension.
830 : !> \param[in] np : The length of the input `Data` matrix along the second dimension.
831 : !> \param[in] Data : The input data series data vector.
832 : !> \param[in] Weight : The vector of weights of the input data points (**optional**, default = array of ones).
833 : !> \param[in] tenabled : A logical value that, if `.true.` will cause the output `NormData` to have transposed shape
834 : !> of the input `Point(nd,np)` matrix, that is `(np,nd)` (**optional**, default = `.false.`).
835 : !>
836 : !> \return
837 : !> `NormData` : The integrated autocorrelation (IAC) via the BatchMeans method.
838 : !>
839 : !> \remark
840 : !> Note that np must be large enough to get a meaningful answer.
841 57 : pure function normalizeWeightedData_2D(nd, np, Data, Weight, tenabled) result(NormData)
842 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
843 : !DEC$ ATTRIBUTES DLLEXPORT :: normalizeWeightedData_2D
844 : #endif
845 3 : use Constants_mod, only: IK, RK
846 : implicit none
847 : integer(IK) , intent(in) :: nd, np
848 : real(RK) , intent(in) :: Data(nd,np)
849 : integer(IK) , intent(in), optional :: Weight(np)
850 : logical , intent(in), optional :: tenabled
851 : real(RK) , allocatable :: NormData(:,:)
852 : logical :: tenabledDefault
853 153 : real(RK) :: Mean(nd)
854 51 : real(RK) :: sumWeight
855 : integer(IK) :: id, ip
856 :
857 51 : tenabledDefault = .false.
858 51 : if (present(tenabled)) tenabledDefault = tenabled
859 :
860 102 : Mean = 0._RK
861 51 : if (present(Weight)) then
862 6 : sumWeight = 0._RK
863 30261 : do ip = 1, np
864 30255 : sumWeight = sumWeight + Weight(ip)
865 60516 : Mean = Mean + Weight(ip) * Data(1:nd,ip)
866 : end do
867 : else
868 45 : sumWeight = np
869 449370 : do ip = 1, np
870 898695 : Mean = Mean + Data(1:nd,ip)
871 : end do
872 : end if
873 102 : Mean = Mean / sumWeight
874 :
875 51 : if (tenabledDefault) then
876 6 : allocate(NormData(np,nd))
877 6 : do concurrent(id = 1:nd)
878 30267 : NormData(1:np,id) = Data(id,1:np) - Mean(id)
879 : end do
880 : else
881 45 : allocate(NormData(nd,np))
882 45 : do concurrent(id = 1:nd)
883 449415 : NormData(id,1:np) = Data(id,1:np) - Mean(id)
884 : end do
885 : end if
886 :
887 102 : end function normalizeWeightedData_2D
888 :
889 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
890 :
891 : !> \brief
892 : !> Return the variance of the input vector of values of length `np`.
893 : !>
894 : !> \param[in] np : The number of points in the sample.
895 : !> \param[in] mean : The mean of the vector.
896 : !> \param[in] Point : The array of shape `np` containing the sample.
897 : !> \param[in] Weight : The vector of weights of the points in the sample (**optional**).
898 : !> \param[in] sumWeight : The sum of the vector of weights (**optional**, if `Weight` is also missing).
899 : !>
900 : !> \return
901 : !> `variance` : The variance of the input sample.
902 : !>
903 : !> \warning
904 : !> If `Weight` is present, then `sumWeight` must be also present.
905 : !> Why? if mean is already given, it implies that sumWeight is also computed a priori.
906 : !>
907 : !> \remark
908 : !> One also use the concise Fortran syntax to achieve the same goal as this function:
909 : !> ```
910 : !> mean = sum(Weight*Point) / sum(Weight)
911 : !> variance = sum( (Weight*(Point-mean))**2 ) / (sum(Weight)-1)
912 : !> ```
913 : !> But the above concise version will be slightly slower as it involves three loops instead of two.
914 9 : pure function getVariance_1D(np,mean,Point,Weight,sumWeight) result(variance)
915 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
916 : !DEC$ ATTRIBUTES DLLEXPORT :: getVariance_1D
917 : #endif
918 : implicit none
919 : integer(IK), intent(in) :: np ! np is the number of observations (points) whose variance is to be computed
920 : real(RK) , intent(in) :: mean ! the mean value of the vector Point
921 : real(RK) , intent(in) :: Point(np) ! Point is the vector of data
922 : integer(IK), intent(in), optional :: Weight(np), sumWeight ! sample weight
923 : real(RK) :: variance ! output mean vector
924 : integer(IK) :: ip
925 6 : variance = 0._RK
926 6 : if (present(Weight)) then
927 48 : do ip = 1,np
928 48 : variance = variance + Weight(ip) * ( Point(ip) - mean )**2
929 : end do
930 3 : variance = variance / real(sumWeight-1_IK,kind=RK)
931 : else
932 48 : do ip = 1,np
933 48 : variance = variance + ( Point(ip) - mean )**2
934 : end do
935 3 : variance = variance / real(np-1_IK,kind=RK)
936 : end if
937 57 : end function getVariance_1D
938 :
939 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
940 :
941 : !> \brief
942 : !> Return the vector of variance of a set of `np` points of `nd` dimensions.
943 : !>
944 : !> \param[in] nd : The number of dimensions of the input sample.
945 : !> \param[in] np : The number of points in the sample.
946 : !> \param[in] Mean : The mean vector of the sample.
947 : !> \param[in] Point : The array of shape `(nd,np)` containing the sample.
948 : !> \param[in] Weight : The vector of weights of the points in the sample of shape `(nd,np)` (**optional**).
949 : !>
950 : !> \return
951 : !> `Variance` : The vector of length `nd` of the variances of the input sample.
952 36 : pure function getVariance_2D(nd,np,Mean,Point,Weight) result(Variance)
953 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
954 : !DEC$ ATTRIBUTES DLLEXPORT :: getVariance_2D
955 : #endif
956 : ! returns the variance of each row in Point weighted by the corresponding Weight if provided
957 : ! pass the Mean vector by calling getMean() or getMean_2D()
958 : implicit none
959 : integer(IK), intent(in) :: nd, np ! np is the number of observations (points) whose variance is to be computed
960 : real(RK) , intent(in) :: Mean(nd) ! the Mean value of the vector Point
961 : real(RK) , intent(in) :: Point(nd,np) ! Point is the vector of data
962 : integer(IK), intent(in), optional :: Weight(np) ! sample weight
963 : real(RK) :: Variance(nd) ! output Mean vector
964 : integer(IK) :: ip, sumWeight
965 24 : Variance = 0._RK
966 6 : if (present(Weight)) then
967 3 : sumWeight = 0_IK
968 18 : do ip = 1,np
969 15 : sumWeight = sumWeight + Weight(ip)
970 63 : Variance = Variance + Weight(ip) * ( Point(1:nd,ip) - Mean )**2
971 : end do
972 12 : Variance = Variance / real(sumWeight-1_IK,kind=RK)
973 : else
974 18 : do ip = 1,np
975 63 : Variance = Variance + ( Point(1:nd,ip) - Mean )**2
976 : end do
977 12 : Variance = Variance / real(np-1_IK,kind=RK)
978 : end if
979 12 : end function getVariance_2D
980 :
981 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
982 :
983 : !> \brief
984 : !> Return the lower triangle Cholesky Factor of the covariance matrix of a set of points in the lower part of `CholeskyLower`.
985 : ! The upper part of `CholeskyLower`, including the diagonal elements of it, will contain the covariance matrix of the sample.
986 : ! The output argument `Diagonal`, contains the diagonal terms of Cholesky Factor.
987 : !>
988 : !> \param[in] nd : The number of dimensions of the input sample.
989 : !> \param[in] np : The number of points in the sample.
990 : !> \param[in] Mean : The mean vector of the sample.
991 : !> \param[in] Point : The array of shape `(nd,np)` containing the sample.
992 : !> \param[out] CholeskyLower : The output matrix of shape `(nd,nd)` whose lower triangle contains elements of the Cholesky factor.
993 : !> The upper triangle of the matrix contains the covariance matrix of the sample.
994 : !> \param[out] Diagonal : The diagonal elements of the Cholesky factor.
995 : !>
996 : !> \todo
997 : !> The efficiency of this code can be further improved.
998 3 : subroutine getSamCholFac(nd,np,Mean,Point,CholeskyLower,Diagonal)
999 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1000 : !DEC$ ATTRIBUTES DLLEXPORT :: getSamCholFac
1001 : #endif
1002 6 : use Matrix_mod, only: getCholeskyFactor
1003 : implicit none
1004 : integer(IK), intent(in) :: nd,np ! np is the number of observations, nd is the number of parameters for each observation
1005 : real(RK) , intent(in) :: Mean(nd) ! Mean vector
1006 : real(RK) , intent(in) :: Point(nd,np) ! Point is the matrix of the data, CovMat contains the elements of the sample covariance matrix
1007 : real(RK) , intent(out) :: CholeskyLower(nd,nd) ! Lower Cholesky Factor of the covariance matrix
1008 : real(RK) , intent(out) :: Diagonal(nd) ! Diagonal elements of the Cholesky factorization
1009 57 : real(RK) :: NormData(np,nd), npMinusOneInverse
1010 : integer(IK) :: i,j
1011 :
1012 18 : do i = 1,np
1013 63 : NormData(i,1:nd) = Point(1:nd,i) - Mean
1014 : end do
1015 :
1016 : ! Only upper half of CovMat is needed
1017 3 : npMinusOneInverse = 1._RK / real(np-1,kind=RK)
1018 12 : do j = 1,nd
1019 30 : do i = 1,j
1020 : ! Get the covariance matrix elements: only the upper half of CovMat is needed
1021 117 : CholeskyLower(i,j) = dot_product( NormData(1:np,i) , NormData(1:np,j) ) * npMinusOneInverse
1022 : end do
1023 : end do
1024 :
1025 : ! XXX: The efficiency can be improved by merging it with the above loops
1026 3 : call getCholeskyFactor(nd,CholeskyLower,Diagonal)
1027 :
1028 3 : end subroutine getSamCholFac
1029 :
1030 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1031 :
1032 : !> \brief
1033 : !> Return the sample mean, covariance matrix, the Mahalanobis distances squared of the points with respect to the sample,
1034 : !> and the square root of the determinant of the inverse covariance matrix of the sample.
1035 : !>
1036 : !> \param[in] nd : The number of dimensions of the input sample.
1037 : !> \param[in] np : The number of points in the sample.
1038 : !> \param[in] Point : The array of shape `(np,nd)` containing the sample.
1039 : !> \param[out] CovMat : The output matrix of shape `(nd,nd)` representing the covariance matrix of the input data.
1040 : !> \param[out] Mean : The output mean vector of the sample.
1041 : !> \param[out] MahalSq : The output Mahalanobis distances squared of the points with respect to the sample (**optional**).
1042 : !> \param[out] InvCovMat : The output inverse covariance matrix of the input data (**optional**).
1043 : !> \param[out] sqrtDetInvCovMat : The output square root of the determinant of the inverse covariance matrix of the sample (**optional**).
1044 : !>
1045 : !> \warning
1046 : !> If `sqrtDetInvCovMat` is present, then `MahalSq` and `InvCovMat` must be also present.
1047 : !>
1048 : !> \remark
1049 : !> Note the shape of the input `Point(np,nd)`.
1050 : !>
1051 : !> \remark
1052 : !> See also, [getSamCovMeanTrans](@ref getsamcovmeantrans).
1053 : !>
1054 : !> \remark
1055 : !> For more information, see Geisser & Cornfield (1963) "Posterior distributions for multivariate normal parameters".
1056 : !> Also, Box and Tiao (1973), "Bayesian Inference in Statistical Analysis" Page 421.
1057 : !>
1058 : !> \author
1059 : !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
1060 9 : subroutine getSamCovMean(np,nd,Point,CovMat,Mean,MahalSq,InvCovMat,sqrtDetInvCovMat)
1061 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1062 : !DEC$ ATTRIBUTES DLLEXPORT :: getSamCovMean
1063 : #endif
1064 3 : use Matrix_mod, only: getInvPosDefMatSqrtDet
1065 : implicit none
1066 : integer(IK), intent(in) :: np,nd ! np is the number of observations, nd is the number of parameters for each observation
1067 : real(RK) , intent(in) :: Point(np,nd) ! Point is the matrix of the data, CovMat contains the elements of the sample covariance matrix
1068 : real(RK) , intent(out) :: CovMat(nd,nd) ! Covariance matrix of the input data
1069 : real(RK) , intent(out) :: Mean(nd) ! Mean vector
1070 : real(RK) , intent(out), optional :: MahalSq(np) ! Vector of Mahalanobis Distances Squared, with respect to the mean position of the sample
1071 : real(RK) , intent(out), optional :: InvCovMat(nd,nd) ! Inverse Covariance matrix of the input data
1072 : real(RK) , intent(out), optional :: sqrtDetInvCovMat ! sqrt determinant of the inverse covariance matrix
1073 57 : real(RK) :: NormData(np,nd)
1074 15 : real(RK) :: DummyVec(nd)
1075 : integer(IK) :: i,j
1076 :
1077 12 : do j = 1,nd
1078 54 : Mean(j) = sum(Point(1:np,j)) / real(np,kind=RK)
1079 57 : NormData(1:np,j) = Point(1:np,j) - Mean(j)
1080 : end do
1081 12 : do i = 1,nd
1082 39 : do j = 1,nd
1083 171 : CovMat(i,j) = dot_product(NormData(1:np,i),NormData(1:np,j))/real(np-1,kind=RK)
1084 : end do
1085 : end do
1086 :
1087 3 : if (present(sqrtDetInvCovMat)) then ! Calculate InCovMat, MahalSq, and sqrt Determinant of InCovMat
1088 12 : do j = 1,nd
1089 30 : do i = 1,j
1090 27 : InvCovMat(i,j) = CovMat(i,j) ! Only the upper half of CovMat is needed
1091 : end do
1092 : end do
1093 3 : call getInvPosDefMatSqrtDet(nd,InvCovMat,sqrtDetInvCovMat)
1094 18 : do i = 1,np
1095 60 : do j = 1,nd
1096 195 : DummyVec(j) = dot_product(InvCovMat(1:nd,j),NormData(i,1:nd))
1097 : end do
1098 63 : MahalSq(i) = dot_product(NormData(i,1:nd),DummyVec)
1099 : end do
1100 : end if
1101 :
1102 3 : end subroutine getSamCovMean
1103 :
1104 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1105 :
1106 : !> \brief
1107 : !> Return the sample mean, covariance matrix, the Mahalanobis distances squared of the points with respect to the sample,
1108 : !> and the square root of the determinant of the inverse covariance matrix of the sample.
1109 : !>
1110 : !> \param[in] nd : The number of dimensions of the input sample.
1111 : !> \param[in] np : The number of points in the sample.
1112 : !> \param[in] Point : The array of shape `(nd,np)` containing the sample.
1113 : !> \param[out] CovMat : The output matrix of shape `(nd,nd)` representing the covariance matrix of the input data.
1114 : !> \param[out] Mean : The output mean vector of the sample.
1115 : !> \param[out] MahalSq : The output Mahalanobis distances squared of the points with respect to the sample (**optional**).
1116 : !> \param[out] InvCovMat : The output inverse covariance matrix of the input data (**optional**).
1117 : !> \param[out] sqrtDetInvCovMat : The output square root of the determinant of the inverse covariance matrix of the sample (**optional**).
1118 : !>
1119 : !> \warning
1120 : !> If `sqrtDetInvCovMat` is present, then `MahalSq` and `InvCovMat` must be also present.
1121 : !>
1122 : !> \remark
1123 : !> Note the shape of the input `Point(nd,np)`.
1124 : !>
1125 : !> \remark
1126 : !> This subroutine has the same functionality as [getSamCovMean](@ref getsamcovmean), with the only difference that input data is transposed here on input.
1127 : !> Based on the preliminary benchmarks with Intel 2017 ifort, `getSamCovMean()` is slightly faster than `getSamCovMeanTrans()`.
1128 : !>
1129 : !> \remark
1130 : !> For more information, see Geisser & Cornfield (1963) "Posterior distributions for multivariate normal parameters".
1131 : !> Also, Box and Tiao (1973), "Bayesian Inference in Statistical Analysis" Page 421.
1132 : !>
1133 : !> \author
1134 : !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
1135 18 : subroutine getSamCovMeanTrans(np,nd,Point,CovMat,Mean,MahalSq,InvCovMat,sqrtDetInvCovMat)
1136 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1137 : !DEC$ ATTRIBUTES DLLEXPORT :: getSamCovMeanTrans
1138 : #endif
1139 :
1140 3 : use Matrix_mod, only: getInvPosDefMatSqrtDet
1141 : implicit none
1142 : integer(IK), intent(in) :: np,nd ! np is the number of observations, nd is the number of parameters for each observation
1143 : real(RK) , intent(in) :: Point(nd,np) ! Point is the matrix of the data, CovMat contains the elements of the sample covariance matrix
1144 : real(RK) , intent(out) :: CovMat(nd,nd) ! Covariance matrix of the input data
1145 : real(RK) , intent(out) :: Mean(nd) ! Mean vector
1146 : real(RK) , intent(out), optional :: MahalSq(np) ! Vector of Mahalanobis Distances Squared, with respect to the mean position of the sample
1147 : real(RK) , intent(out), optional :: InvCovMat(nd,nd) ! Inverse Covariance matrix of the input data
1148 : real(RK) , intent(out), optional :: sqrtDetInvCovMat ! sqrt determinant of the inverse covariance matrix
1149 60 : real(RK) , dimension(nd) :: DummyVec
1150 432 : real(RK) , dimension(nd,np) :: NormData
1151 : integer(IK) :: i,j
1152 :
1153 48 : Mean = 0._RK
1154 117 : do i = 1,np
1155 432 : do j = 1,nd
1156 420 : Mean(j) = Mean(j) + Point(j,i)
1157 : end do
1158 : end do
1159 48 : Mean = Mean / real(np,kind=RK)
1160 :
1161 117 : do i = 1,np
1162 432 : NormData(1:nd,i) = Point(1:nd,i) - Mean
1163 : end do
1164 :
1165 48 : do i = 1,nd
1166 156 : do j = 1,nd
1167 1089 : CovMat(i,j) = dot_product(NormData(i,1:np),NormData(j,1:np)) / real(np-1,kind=RK)
1168 : end do
1169 : end do
1170 :
1171 12 : if (present(sqrtDetInvCovMat)) then ! Calculate InCovMat, MahalSq, and sqrt Determinant of InCovMat
1172 12 : do j = 1,nd
1173 30 : do i = 1,j
1174 27 : InvCovMat(i,j) = CovMat(i,j) ! Only the upper half of CovMat is needed
1175 : end do
1176 : end do
1177 3 : call getInvPosDefMatSqrtDet(nd,InvCovMat,sqrtDetInvCovMat)
1178 18 : do i = 1,np
1179 60 : do j = 1,nd
1180 195 : DummyVec(j) = dot_product(InvCovMat(1:nd,j),NormData(1:nd,i))
1181 : end do
1182 63 : MahalSq(i) = dot_product(NormData(1:nd,i),DummyVec)
1183 : !MahalSq = dot_product(NormData(1:nd,i),DummyVec)
1184 : !if (maxMahal<MahalSq) maxMahal = MahalSq
1185 : end do
1186 : !maxMahal = maxval(MahalSq)
1187 : end if
1188 :
1189 12 : end subroutine getSamCovMeanTrans
1190 :
1191 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1192 :
1193 : !> \brief
1194 : !> Return the sample mean and the upper triangle of the covariance matrix of the input sample.
1195 : !>
1196 : !> \param[in] nd : The number of dimensions of the input sample.
1197 : !> \param[in] np : The number of points in the sample.
1198 : !> \param[in] Point : The array of shape `(nd,np)` containing the sample.
1199 : !> \param[out] CovMatUpper : The output matrix of shape `(nd,nd)` whose upper triangle represents the covariance matrix of the input data.
1200 : !> \param[out] Mean : The output mean vector of the sample.
1201 : !>
1202 : !> \remark
1203 : !> Note the shape of the input `Point(nd,np)`.
1204 : !>
1205 : !> \remark
1206 : !> This subroutine has the same functionality as [getSamCovMeanTrans](@ref getsamcovmeantrans), with the only difference
1207 : !> only the upper triangle of the covariance matrix is returned. Also, optional arguments are not available.
1208 : !> This subroutine is specifically optimized for use in the ParaMonte samplers.
1209 : !>
1210 : !> \remark
1211 : !> For more information, see Geisser & Cornfield (1963) "Posterior distributions for multivariate normal parameters".
1212 : !> Also, Box and Tiao (1973), "Bayesian Inference in Statistical Analysis" Page 421.
1213 : !>
1214 : !> \author
1215 : !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
1216 63 : subroutine getSamCovUpperMeanTrans(np,nd,Point,CovMatUpper,Mean)
1217 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1218 : !DEC$ ATTRIBUTES DLLEXPORT :: getSamCovUpperMeanTrans
1219 : #endif
1220 12 : use Matrix_mod, only: getInvPosDefMatSqrtDet
1221 : implicit none
1222 : integer(IK), intent(in) :: np,nd ! np is the number of observations, nd is the number of parameters for each observation
1223 : real(RK) , intent(in) :: Point(nd,np) ! Point is the matrix of the data, CovMatUpper contains the elements of the sample covariance matrix
1224 : real(RK) , intent(out) :: CovMatUpper(nd,nd) ! Covariance matrix of the input data
1225 : real(RK) , intent(out) :: Mean(nd) ! Mean vector
1226 1165 : real(RK) :: npMinusOneInvReal, NormData(nd,np)
1227 : integer(IK) :: i,j
1228 :
1229 190 : Mean = 0._RK
1230 386 : do i = 1,np
1231 1165 : do j = 1,nd
1232 1102 : Mean(j) = Mean(j) + Point(j,i)
1233 : end do
1234 : end do
1235 190 : Mean = Mean / real(np,kind=RK)
1236 :
1237 386 : do i = 1,np
1238 1165 : NormData(1:nd,i) = Point(1:nd,i) - Mean
1239 : end do
1240 :
1241 63 : npMinusOneInvReal = 1._RK / real(np-1,kind=RK)
1242 190 : do j = 1,nd
1243 402 : do i = 1,j
1244 1769 : CovMatUpper(i,j) = dot_product(NormData(i,1:np),NormData(j,1:np)) * npMinusOneInvReal
1245 : end do
1246 : end do
1247 :
1248 63 : end subroutine getSamCovUpperMeanTrans
1249 :
1250 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1251 :
1252 : !> \brief
1253 : !> Return the mean and the upper triangle of the covariance matrix of the input *weighted* sample.
1254 : !>
1255 : !> \param[in] nd : The number of dimensions of the input sample.
1256 : !> \param[in] sumWeight : The sum of all sample weights.
1257 : !> \param[in] np : The number of points in the sample.
1258 : !> \param[in] Point : The array of shape `(nd,np)` containing the sample.
1259 : !> \param[in] Weight : The integer array of length `np`, representing the weights of individual points in the sample.
1260 : !> \param[out] CovMatUpper : The output matrix of shape `(nd,nd)` whose upper triangle represents the covariance matrix of the input data.
1261 : !> \param[out] Mean : The output mean vector of the sample.
1262 : !>
1263 : !> \warning
1264 : !> Note the shape of the input argument `Point(nd,np)`.
1265 : !>
1266 : !> \remark
1267 : !> This subroutine has the same functionality as [getSamCovUpperMeanTrans](@ref getsamcovuppermeantrans), with the only difference
1268 : !> only the upper triangle of the covariance matrix is returned. Also, optional arguments are not available.
1269 : !> This subroutine is specifically optimized for use in the ParaMonte samplers.
1270 : !>
1271 : !> \remark
1272 : !> For more information, see Geisser & Cornfield (1963) "Posterior distributions for multivariate normal parameters".
1273 : !> Also, Box and Tiao (1973), "Bayesian Inference in Statistical Analysis" Page 421.
1274 : !>
1275 : !> \author
1276 : !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
1277 30904 : subroutine getWeiSamCovUppMeanTrans(np,sumWeight,nd,Point,Weight,CovMatUpper,Mean)
1278 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1279 : !DEC$ ATTRIBUTES DLLEXPORT :: getWeiSamCovUppMeanTrans
1280 : #endif
1281 :
1282 63 : use Matrix_mod, only: getInvPosDefMatSqrtDet
1283 : implicit none
1284 : integer(IK), intent(in) :: np,nd,sumWeight ! np is the number of observations, nd is the number of parameters for each observation
1285 : integer(IK), intent(in) :: Weight(np) ! weights of the points
1286 : real(RK) , intent(in) :: Point(nd,np) ! Point is the matrix of the data, CovMatUpper contains the elements of the sample covariance matrix
1287 : real(RK) , intent(out) :: CovMatUpper(nd,nd) ! Covariance matrix of the input data
1288 : real(RK) , intent(out) :: Mean(nd) ! Mean vector
1289 30904 : real(RK) :: sumWeightMinusOneInvReal
1290 694940 : real(RK) :: NormData(nd,np)
1291 : integer(IK) :: i,j,ip
1292 :
1293 80805 : Mean = 0._RK
1294 277546 : do i = 1,np
1295 694940 : do j = 1,nd
1296 664036 : Mean(j) = Mean(j) + Weight(i)*Point(j,i)
1297 : end do
1298 : end do
1299 80805 : Mean = Mean / real(sumWeight,kind=RK)
1300 :
1301 277546 : do i = 1,np
1302 694940 : NormData(1:nd,i) = Point(1:nd,i) - Mean
1303 : end do
1304 :
1305 30904 : sumWeightMinusOneInvReal = 1._RK / real(sumWeight-1,kind=RK)
1306 80805 : do j = 1,nd
1307 149706 : do i = 1,j
1308 68901 : CovMatUpper(i,j) = 0
1309 657062 : do ip = 1,np
1310 657062 : CovMatUpper(i,j) = CovMatUpper(i,j) + Weight(ip)*NormData(i,ip)*NormData(j,ip)
1311 : end do
1312 118802 : CovMatUpper(i,j) = CovMatUpper(i,j) * sumWeightMinusOneInvReal
1313 : end do
1314 : end do
1315 :
1316 30904 : end subroutine getWeiSamCovUppMeanTrans
1317 :
1318 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1319 :
1320 : !> \brief
1321 : !> Given two input sample means and covariance matrices, return the combination of them as a single mean and covariance matrix.
1322 : !>
1323 : !> \param[in] nd : The number of dimensions of the input sample.
1324 : !> \param[in] npA : The number of points in sample `A`.
1325 : !> \param[in] MeanVecA : The mean vector of sample `A`.
1326 : !> \param[in] CovMatA : The covariance matrix of sample `A`.
1327 : !> \param[in] npB : The number of points in sample `B`.
1328 : !> \param[in] MeanVecB : The mean vector of sample `B`.
1329 : !> \param[in] CovMatB : The covariance matrix of sample `B`.
1330 : !> \param[out] MeanVecAB : The output mean vector of the combined sample.
1331 : !> \param[out] CovMatAB : The output covariance matrix of the combined sample.
1332 : !>
1333 : !> \author
1334 : !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
1335 : !>
1336 : !> \remark
1337 : !> An exact implementation of this algorithm which needs only the upper triangles of the input matrices and
1338 : !> yields only the upper triangle of the covariance matrix is given in [mergeMeanCovUpper](@ref mergemeancovupper).
1339 : !> The alternative implementation is much more efficient, by a factor of 6-7 with all compiler optimization flags on.
1340 : !>
1341 : ! This subroutine uses a recursion equation similar to http://stats.stackexchange.com/questions/97642/how-to-combine-sample-means-and-sample-variances
1342 : ! See also: O’Neill, (2014), "Some Useful Moment Results in Sampling Problems".
1343 : ! See also: https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Covariance
1344 : ! See also: https://stats.stackexchange.com/questions/43159/how-to-calculate-pooled-variance-of-two-groups-given-known-group-variances-mean
1345 3 : subroutine mergeMeanCov(nd,npA,MeanVecA,CovMatA,npB,MeanVecB,CovMatB,MeanVecAB,CovMatAB)
1346 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1347 : !DEC$ ATTRIBUTES DLLEXPORT :: mergeMeanCov
1348 : #endif
1349 : implicit none
1350 : integer(IK), intent(in) :: nd
1351 : integer(IK), intent(in) :: npA,npB
1352 : real(RK) , intent(in) :: MeanVecA(nd),CovMatA(nd,nd)
1353 : real(RK) , intent(in) :: MeanVecB(nd),CovMatB(nd,nd)
1354 : real(RK) , intent(out) :: MeanVecAB(nd),CovMatAB(nd,nd)
1355 42 : real(RK) , dimension(nd,1) :: MeanMatA,MeanMatB,MeanMat
1356 42 : real(RK) :: DistanceSq(nd,nd)
1357 3 : real(RK) :: npABinverse
1358 : integer(IK) :: npAB
1359 :
1360 3 : npAB = npA + npB
1361 3 : npABinverse = 1._RK / real(npAB, kind=RK)
1362 12 : MeanMatA(1:nd,1) = MeanVecA
1363 12 : MeanMatB(1:nd,1) = MeanVecB
1364 :
1365 : ! Compute the new Mean
1366 :
1367 12 : MeanVecAB = ( npA * MeanVecA + npB * MeanVecB ) * npABinverse
1368 12 : MeanMat(1:nd,1) = MeanVecAB
1369 :
1370 : ! Compute the new Covariance matrix
1371 :
1372 69 : DistanceSq = matmul( (MeanMatA-MeanMatB), transpose((MeanMatA-MeanMatB)) ) * npA * npB * npABinverse
1373 39 : CovMatAB = ( (npA-1) * CovMatA + (npB-1) * CovMatB + DistanceSq ) / real(npAB-1, kind=RK)
1374 :
1375 30904 : end subroutine mergeMeanCov
1376 :
1377 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1378 :
1379 : ! !> \brief
1380 : ! !> Given two input sample means and covariance matrices, return the combination of them as a single mean and covariance matrix.
1381 : ! !>
1382 : ! !> \param[in] nd : The number of dimensions of the input sample.
1383 : ! !> \param[in] npA : The number of points in sample `A`.
1384 : ! !> \param[in] MeanVecA : The mean vector of sample `A`.
1385 : ! !> \param[in] CovMatUpperA : The covariance matrix of sample `A`.
1386 : ! !> \param[in] npB : The number of points in sample `B`.
1387 : ! !> \param[in] MeanVecB : The mean vector of sample `B`.
1388 : ! !> \param[in] CovMatUpperB : The covariance matrix of sample `B`.
1389 : ! !> \param[out] MeanVecAB : The output mean vector of the combined sample.
1390 : ! !> \param[out] CovMatUpperAB : The output covariance matrix of the combined sample.
1391 : ! !>
1392 : ! !> \todo
1393 : ! !> The efficiency of this algorithm might still be improved by converting the upper triangle covariance matrix to a packed vector.
1394 : ! !>
1395 : ! !> \remark
1396 : ! !> This subroutine is the same as [mergeMeanCov](@ref mergeMeanCov), with the **important difference** that only the
1397 : ! !> upper triangles and diagonals of the input covariance matrices need to be given by the user: `CovMatUpperA`, `CovMatUpperB`
1398 : ! !> This alternative implementation is 6-7 times faster, with all compiler optimization flags on.
1399 : ! !>
1400 : ! !> \author
1401 : ! !> Amir Shahmoradi, Nov 24, 2020, 4:19 AM, Dallas, TX
1402 : ! subroutine mergeMeanCovUpperSlow(nd,npA,MeanVecA,CovMatUpperA,npB,MeanVecB,CovMatUpperB,MeanVecAB,CovMatUpperAB)
1403 : !#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1404 : ! !DEC$ ATTRIBUTES DLLEXPORT :: mergeMeanCovUpperSlow
1405 : !#endif
1406 : ! implicit none
1407 : ! integer(IK), intent(in) :: nd
1408 : ! integer(IK), intent(in) :: npA,npB
1409 : ! real(RK) , intent(in) :: MeanVecA(nd),CovMatUpperA(nd,nd)
1410 : ! real(RK) , intent(in) :: MeanVecB(nd),CovMatUpperB(nd,nd)
1411 : ! real(RK) , intent(out) :: MeanVecAB(nd),CovMatUpperAB(nd,nd)
1412 : ! real(RK) :: npABinverse, npAnpB2npAB
1413 : ! real(RK) :: npA2npAB, npB2npAB
1414 : ! real(RK) :: MeanVecDiffAB(nd)
1415 : ! integer(IK) :: npAB, i, j
1416 : !
1417 : ! npAB = npA + npB
1418 : ! npABinverse = 1._RK / real(npAB, kind=RK)
1419 : ! npAnpB2npAB = npA * npB * npABinverse
1420 : ! npA2npAB = npA * npABinverse
1421 : ! npB2npAB = npB * npABinverse
1422 : !
1423 : ! ! Compute the new Mean and Covariance matrix
1424 : !
1425 : ! do j = 1, nd
1426 : ! MeanVecDiffAB(j) = MeanVecA(j) - MeanVecB(j)
1427 : ! !MeanVecAB(j) = ( npA * MeanVecA(j) + npB * MeanVecB(j) ) * npABinverse
1428 : ! MeanVecAB(j) = npA2npAB * MeanVecA(j) + npB2npAB * MeanVecB(j)
1429 : ! do i = 1, j
1430 : ! CovMatUpperAB(i,j) = ( (npA-1) * CovMatUpperA(i,j) + (npB-1) * CovMatUpperB(i,j) + MeanVecDiffAB(i) * MeanVecDiffAB(j) * npAnpB2npAB ) / real(npAB-1, kind=RK)
1431 : ! end do
1432 : ! end do
1433 : !
1434 : !
1435 : ! end subroutine mergeMeanCovUpperSlow
1436 :
1437 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1438 :
1439 : !> \brief
1440 : !> Given two input sample means and covariance matrices, return the combination of them as a single mean and covariance matrix.
1441 : !>
1442 : !> \param[in] nd : The number of dimensions of the input sample.
1443 : !> \param[in] npA : The number of points in sample `A`.
1444 : !> \param[in] MeanVecA : The mean vector of sample `A`.
1445 : !> \param[in] CovMatUpperA : The covariance matrix of sample `A`.
1446 : !> \param[in] npB : The number of points in sample `B`.
1447 : !> \param[in] MeanVecB : The mean vector of sample `B`.
1448 : !> \param[in] CovMatUpperB : The covariance matrix of sample `B`.
1449 : !> \param[out] MeanVecAB : The output mean vector of the combined sample.
1450 : !> \param[out] CovMatUpperAB : The output covariance matrix of the combined sample.
1451 : !>
1452 : !> \remark
1453 : !> This subroutine is the same as [mergeMeanCov](@ref mergemeancov), with the **important difference** that only the
1454 : !> upper triangles and diagonals of the input covariance matrices need to be given by the user: `CovMatUpperA`, `CovMatUpperB`
1455 : !> This alternative implementation is 6-7 times faster, with all compiler optimization flags on.
1456 : !> In addition, all computational coefficients are predefined in this implementation,
1457 : !> resulting in an extra 10%-15% efficiency gain.
1458 : !>
1459 : !> \todo
1460 : !> The efficiency of this algorithm might still be improved by converting the upper triangle covariance matrix to a packed vector.
1461 : !>
1462 : !> \author
1463 : !> Amir Shahmoradi, Nov 24, 2020, 4:19 AM, Dallas, TX
1464 30173 : subroutine mergeMeanCovUpper(nd,npA,MeanVecA,CovMatUpperA,npB,MeanVecB,CovMatUpperB,MeanVecAB,CovMatUpperAB)
1465 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1466 : !DEC$ ATTRIBUTES DLLEXPORT :: mergeMeanCovUpper
1467 : #endif
1468 : implicit none
1469 : integer(IK), intent(in) :: nd
1470 : integer(IK), intent(in) :: npA,npB
1471 : real(RK) , intent(in) :: MeanVecA(nd),CovMatUpperA(nd,nd)
1472 : real(RK) , intent(in) :: MeanVecB(nd),CovMatUpperB(nd,nd)
1473 : real(RK) , intent(out) :: MeanVecAB(nd),CovMatUpperAB(nd,nd)
1474 78908 : real(RK) :: MeanVecDiffAB(nd)
1475 30173 : real(RK) :: npAnpB2npAB2npABMinusOne
1476 30173 : real(RK) :: npAMinusOne2npABMinusOne
1477 30173 : real(RK) :: npBMinusOne2npABMinusOne
1478 30173 : real(RK) :: npABMinusOneInverse
1479 30173 : real(RK) :: npABinverse
1480 30173 : real(RK) :: npA2npAB
1481 30173 : real(RK) :: npB2npAB
1482 : integer(IK) :: npAB, i, j
1483 :
1484 30173 : npAB = npA + npB
1485 30173 : npABinverse = 1._RK / real(npAB, kind=RK)
1486 30173 : npABMinusOneInverse = 1._RK / real(npAB-1, kind=RK)
1487 30173 : npAnpB2npAB2npABMinusOne = npA * npB * npABinverse * npABMinusOneInverse
1488 30173 : npA2npAB = npA * npABinverse
1489 30173 : npB2npAB = npB * npABinverse
1490 30173 : npAMinusOne2npABMinusOne = (npA - 1_IK) * npABMinusOneInverse
1491 30173 : npBMinusOne2npABMinusOne = (npB - 1_IK) * npABMinusOneInverse
1492 :
1493 : ! Compute the new Mean and Covariance matrix
1494 :
1495 78908 : do j = 1, nd
1496 48735 : MeanVecDiffAB(j) = MeanVecA(j) - MeanVecB(j)
1497 48735 : MeanVecAB(j) = npA2npAB * MeanVecA(j) + npB2npAB * MeanVecB(j)
1498 146208 : do i = 1, j
1499 67300 : CovMatUpperAB(i,j) = npAMinusOne2npABMinusOne * CovMatUpperA(i,j) &
1500 67300 : + npBMinusOne2npABMinusOne * CovMatUpperB(i,j) &
1501 116035 : + npAnpB2npAB2npABMinusOne * MeanVecDiffAB(i) * MeanVecDiffAB(j)
1502 : end do
1503 : end do
1504 :
1505 :
1506 3 : end subroutine mergeMeanCovUpper
1507 :
1508 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1509 :
1510 : !> \brief
1511 : !> Given two input sample means and covariance matrices, return the combination of them as a single mean and covariance matrix.
1512 : !>
1513 : !> \param[in] nd : The number of dimensions of the input sample.
1514 : !> \param[in] npA : The number of points in sample `A`.
1515 : !> \param[in] MeanVecA : The mean vector of sample `A`.
1516 : !> \param[in] CovMatUpperA : The covariance matrix of sample `A`.
1517 : !> \param[in] npB : The number of points in sample `B`.
1518 : !> \param[inout] MeanVecB : The mean vector of sample `B`.
1519 : !> \param[inout] CovMatUpperB : The covariance matrix of sample `B`.
1520 : !>
1521 : !> \remark
1522 : !> This subroutine is the same as [mergeMeanCovUpper](@ref mergemeancovupper), with the **important difference** that
1523 : !> the resulting output mean and covariance matrices are written to the input arguments `MeanVecB`, `CovMatUpperB`.
1524 : !> This alternative implementation results in another extra 15%-20% efficiency gain. This result is based on
1525 : !> the benchmarks with Intel Fortran compiler 19.4 with all compiler optimization flags on.
1526 : !>
1527 : !> \todo
1528 : !> The efficiency of this algorithm might still be improved by converting the upper triangle covariance matrix to a packed vector.
1529 : !>
1530 : !> \author
1531 : !> Amir Shahmoradi, Nov 25, 2020, 1:00 AM, Dallas, TX
1532 6 : subroutine mergeMeanCovUpperDense(nd,npA,MeanVecA,CovMatUpperA,npB,MeanVecB,CovMatUpperB)
1533 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1534 : !DEC$ ATTRIBUTES DLLEXPORT :: mergeMeanCovUpperDense
1535 : #endif
1536 : implicit none
1537 : integer(IK), intent(in) :: nd
1538 : integer(IK), intent(in) :: npA,npB
1539 : real(RK) , intent(in) :: MeanVecA(nd),CovMatUpperA(nd,nd)
1540 : real(RK) , intent(inout) :: MeanVecB(nd),CovMatUpperB(nd,nd)
1541 24 : real(RK) :: MeanVecDiffAB(nd)
1542 6 : real(RK) :: npAnpB2npAB2npABMinusOne
1543 6 : real(RK) :: npAMinusOne2npABMinusOne
1544 6 : real(RK) :: npBMinusOne2npABMinusOne
1545 6 : real(RK) :: npABMinusOneInverse
1546 6 : real(RK) :: npABinverse
1547 6 : real(RK) :: npA2npAB
1548 6 : real(RK) :: npB2npAB
1549 : integer(IK) :: npAB, i, j
1550 :
1551 6 : npAB = npA + npB
1552 6 : npABinverse = 1._RK / real(npAB, kind=RK)
1553 6 : npABMinusOneInverse = 1._RK / real(npAB-1, kind=RK)
1554 6 : npAnpB2npAB2npABMinusOne = npA * npB * npABinverse * npABMinusOneInverse
1555 6 : npA2npAB = npA * npABinverse
1556 6 : npB2npAB = npB * npABinverse
1557 6 : npAMinusOne2npABMinusOne = (npA - 1_IK) * npABMinusOneInverse
1558 6 : npBMinusOne2npABMinusOne = (npB - 1_IK) * npABMinusOneInverse
1559 :
1560 : ! Compute the new Mean and Covariance matrix
1561 :
1562 24 : do j = 1, nd
1563 18 : MeanVecDiffAB(j) = MeanVecA(j) - MeanVecB(j)
1564 18 : MeanVecB(j) = npA2npAB * MeanVecA(j) + npB2npAB * MeanVecB(j)
1565 60 : do i = 1, j
1566 36 : CovMatUpperB(i,j) = npAMinusOne2npABMinusOne * CovMatUpperA(i,j) &
1567 36 : + npBMinusOne2npABMinusOne * CovMatUpperB(i,j) &
1568 54 : + npAnpB2npAB2npABMinusOne * MeanVecDiffAB(i) * MeanVecDiffAB(j)
1569 : end do
1570 : end do
1571 :
1572 30173 : end subroutine mergeMeanCovUpperDense
1573 :
1574 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1575 :
1576 : !> \brief
1577 : !> Return a random standard Gaussian deviate with zero mean and unit variance.
1578 : !>
1579 : !> \return
1580 : !> `randGaus` : The random standard Gaussian deviate with zero mean and unit variance.
1581 : !>
1582 : !> \remark
1583 : !> See also, Numerical Recipes in Fortran, by Press et al. (1990)
1584 : !>
1585 : !> \author
1586 : !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
1587 898151 : function getRandGaus() result(randGaus)
1588 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1589 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandGaus
1590 : #endif
1591 :
1592 : implicit none
1593 : integer(IK), save :: iset=0
1594 : real(RK) , save :: gset
1595 898151 : real(RK) :: fac,rsq,vec(2)
1596 : real(RK) :: randGaus
1597 :
1598 898151 : if (iset == 0_IK) then
1599 122746 : do
1600 : !block
1601 : !integer :: i, n
1602 : !real(RK) :: unifrnd(30)
1603 : !!integer, dimension(:), allocatable :: seed
1604 : !if (paradramPrintEnabled .or. paradisePrintEnabled) then
1605 : ! !do i = 1, 22
1606 : ! call random_number(unifrnd)
1607 : ! write(*,"(*(g0,:,'"//new_line("a")//"'))") unifrnd
1608 : ! !end do
1609 : ! !call random_seed(size = n); allocate(seed(n))
1610 : ! !call random_seed(get = seed)
1611 : ! !write(*,"(*(g0,:,' '))") seed
1612 : ! !write(*,"(*(g0,:,' '))") StateOld
1613 : ! !write(*,"(*(g0,:,' '))") StateNew
1614 : ! !write(*,"(*(g0,:,' '))") CholeskyLower
1615 : ! !write(*,"(*(g0,:,' '))") domainCheckCounter
1616 : ! paradisePrintEnabled = .false.
1617 : ! paradramPrintEnabled = .false.
1618 : !end if
1619 : !end block
1620 571822 : call random_number(vec)
1621 1715470 : vec = 2._RK*vec - 1._RK
1622 571822 : rsq = vec(1)**2 + vec(2)**2
1623 571822 : if ( rsq > 0._RK .and. rsq < 1._RK ) exit
1624 : end do
1625 449076 : fac = sqrt(-2._RK*log(rsq)/rsq)
1626 449076 : gset = vec(1)*fac
1627 449076 : randGaus = vec(2)*fac
1628 449076 : iset = 1_IK
1629 : else
1630 449075 : randGaus = gset
1631 449075 : iset = 0_IK
1632 : endif
1633 :
1634 898157 : end function getRandGaus
1635 :
1636 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1637 :
1638 : !> \brief
1639 : !> Return a random Gaussian deviate with the given mean and standard deviation.
1640 : !>
1641 : !> \param[in] mean : The mean of the Gaussian distribution.
1642 : !> \param[in] std : The standard deviation of the Gaussian distribution. It must be a positive real number.
1643 : !>
1644 : !> \return
1645 : !> `randNorm` : A normally distributed deviate with the given mean and standard deviation.
1646 : !>
1647 : !> \author
1648 : !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
1649 3 : function getRandNorm(mean,std) result(randNorm)
1650 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1651 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandNorm
1652 : #endif
1653 : implicit none
1654 : real(RK), intent(in) :: mean, std
1655 : real(RK) :: randNorm
1656 6 : randNorm = mean + std * getRandGaus()
1657 898154 : end function getRandNorm
1658 :
1659 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1660 :
1661 : !> \brief
1662 : !> Return a log-normally distributed deviate with the given mean and standard deviation.
1663 : !>
1664 : !> \param[in] mean : The mean of the Lognormal distribution.
1665 : !> \param[in] std : The standard deviation of the Lognormal distribution. It must be a positive real number.
1666 : !>
1667 : !> \return
1668 : !> `randLogn` : A Lognormally distributed deviate with the given mean and standard deviation.
1669 : !>
1670 : !> \author
1671 : !> Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan
1672 3 : function getRandLogn(mean,std) result(randLogn)
1673 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1674 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandLogn
1675 : #endif
1676 : implicit none
1677 : real(RK), intent(in) :: mean, std
1678 : real(RK) :: randLogn
1679 6 : randLogn = exp( mean + std*getRandGaus() )
1680 6 : end function getRandLogn
1681 :
1682 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1683 :
1684 : ! This subroutine is legacy and slow. use getRandMVN() in this same module.
1685 : ! Given the mean vector MeanVec and the covariance matrix CovMat, this subroutine generates a random vector x (of length nd>=2)
1686 : ! from an nd-dimensional multivariate normal distribution.
1687 : ! First a vector of nd standard normal random deviates is generated,
1688 : ! then this vector is multiplied by the Cholesky decomposition of the covariance matrix,
1689 : ! then the vector MeanVec is added to the resulting vector, and is stored in the output vector x.
1690 : ! ATTENTION: Only the upper half of the covariance matrix (plus the diagonal terms) need to be given in the input.
1691 : ! in the ouput, the upper half and diagonal part will still be the covariance matrix, while the lower half will be
1692 : ! the Cholesky decomposition elements (excluding its diagonal terms that are provided only in the vector Diagonal).
1693 : ! USES choldc.f90, getRandGaus.f90
1694 : !> Amir Shahmoradi, March 22, 2012, 2:21 PM, IFS, UTEXAS
1695 3 : subroutine getMVNDev(nd,MeanVec,CovMat,X)
1696 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1697 : !DEC$ ATTRIBUTES DLLEXPORT :: getMVNDev
1698 : #endif
1699 :
1700 : use iso_fortran_env, only: output_unit
1701 3 : use Matrix_mod, only: getCholeskyFactor
1702 :
1703 : implicit none
1704 : integer(IK), intent(in) :: nd
1705 : real(RK) , intent(in) :: MeanVec(nd), CovMat(nd,nd)
1706 : real(RK) , intent(out) :: X(nd)
1707 33 : real(RK) :: CholeskyLower(nd,nd), Diagonal(nd), DummyVec(nd)
1708 : integer(IK) :: i
1709 :
1710 21 : CholeskyLower = CovMat
1711 3 : call getCholeskyFactor(nd,CholeskyLower,Diagonal)
1712 3 : if (Diagonal(1)<0._RK) then
1713 : ! LCOV_EXCL_START
1714 : write(output_unit,"(A)") "getCholeskyFactor() failed in getMVNDev()"
1715 : error stop
1716 : end if
1717 : ! LCOV_EXCL_STOP
1718 9 : do i=1,nd
1719 6 : DummyVec(i) = getRandGaus()
1720 9 : x(i) = DummyVec(i) * Diagonal(i)
1721 : end do
1722 6 : do i=2,nd
1723 9 : x(i) = x(i) + dot_product(CholeskyLower(i,1:i-1),DummyVec(1:i-1))
1724 : end do
1725 9 : x = x + MeanVec
1726 :
1727 3 : end subroutine getMVNDev
1728 :
1729 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1730 :
1731 : ! This subroutine is legacy and slow. use getRandMVU() in this same module.
1732 : ! Given the mean vector MeanVec and the covariance matrix CovMat, this subroutine generates a random vector X (of length nd>=2)
1733 : ! from an nd-dimensional multivariate ellipsoidal uniform distribution, such that getMVUDev() is randomly distributed inside the nd-dimensional ellipsoid.
1734 : ! ATTENTION: Only the upper half of the covariance matrix (plus the diagonal terms) need to be given in the input.
1735 : ! in the ouput, the upper half and diagonal part will still be the covariance matrix, while the lower half will be
1736 : ! the Cholesky decomposition elements (excluding its diagonal terms that are provided only in the vector Diagonal).
1737 : ! USES getCholeskyFactor.f90, getRandGaus.f90
1738 : !> Amir Shahmoradi, April 25, 2016, 2:21 PM, IFS, UTEXAS
1739 3 : subroutine getMVUDev(nd,MeanVec,CovMat,X)
1740 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1741 : !DEC$ ATTRIBUTES DLLEXPORT :: getMVUDev
1742 : #endif
1743 :
1744 3 : use Matrix_mod, only: getCholeskyFactor
1745 :
1746 : implicit none
1747 : integer(IK), intent(in) :: nd
1748 : real(RK) , intent(in) :: MeanVec(nd)
1749 : real(RK) , intent(in) :: CovMat(nd,nd)
1750 : real(RK) , intent(out) :: X(nd)
1751 33 : real(RK) :: Diagonal(nd), DummyVec(nd), CholeskyLower(nd,nd), dummy
1752 : integer(IK) :: i
1753 :
1754 21 : CholeskyLower = CovMat
1755 3 : call getCholeskyFactor(nd,CholeskyLower,Diagonal)
1756 3 : if (Diagonal(1)<0._RK) then
1757 : ! LCOV_EXCL_START
1758 : error stop
1759 : !call abortProgram( output_unit , 1 , 1 , 'Statitistics@getMVUDev()@getCholeskyFactor() failed.' )
1760 : end if
1761 : ! LCOV_EXCL_STOP
1762 9 : do i=1,nd
1763 9 : DummyVec(i) = getRandGaus()
1764 : end do
1765 3 : call random_number(dummy)
1766 9 : dummy = (dummy**(1._RK/real(nd,kind=RK)))/norm2(DummyVec) ! Now DummyVec is a uniformly drawn point from inside of nd-D sphere.
1767 9 : DummyVec = dummy*DummyVec
1768 :
1769 : ! Now transform this point to a point inside the ellipsoid.
1770 9 : do i=1,nd
1771 9 : X(i) = DummyVec(i)*Diagonal(i)
1772 : end do
1773 :
1774 6 : do i=2,nd
1775 9 : X(i) = X(i) + dot_product(CholeskyLower(i,1:i-1),DummyVec(1:i-1))
1776 : end do
1777 :
1778 9 : X = X + MeanVec
1779 :
1780 3 : end subroutine getMVUDev
1781 :
1782 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1783 :
1784 : !> \brief
1785 : !> Return a MultiVariate Normal (MVN) random vector with the given mean and
1786 : !> covariance matrix represented by the the input Cholesky factorization.
1787 : !>
1788 : !> \param[in] nd : The number of dimensions of the MVN distribution.
1789 : !> \param[in] MeanVec : The mean vector of the MVN distribution.
1790 : !> \param[in] CholeskyLower : The Cholesky lower triangle of the covariance matrix of the MVN distribution.
1791 : !> \param[in] Diagonal : The Diagonal elements of the Cholesky lower triangle of the covariance matrix of the MVN distribution.
1792 : !>
1793 : !> \return
1794 : !> `RandMVN` : The randomly generated MVN vector.
1795 : !>
1796 : !> \author
1797 : !> Amir Shahmoradi, April 23, 2017, 12:36 AM, ICES, UTEXAS
1798 1592980 : function getRandMVN(nd,MeanVec,CholeskyLower,Diagonal) result(RandMVN)
1799 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1800 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandMVN
1801 : #endif
1802 : implicit none
1803 : integer(IK), intent(in) :: nd
1804 : real(RK) , intent(in) :: MeanVec(nd)
1805 : real(RK) , intent(in) :: CholeskyLower(nd,nd), Diagonal(nd) ! Cholesky lower triangle and its diagonal terms, calculated from the input CovMat.
1806 338415 : real(RK) :: RandMVN(nd), dummy
1807 : integer(IK) :: j,i
1808 916148 : RandMVN = MeanVec
1809 916148 : do j = 1,nd
1810 577733 : dummy = getRandGaus()
1811 577733 : RandMVN(j) = RandMVN(j) + Diagonal(j) * dummy
1812 1155470 : do i = j+1,nd
1813 817051 : RandMVN(i) = RandMVN(i) + CholeskyLower(i,j) * dummy
1814 : end do
1815 : end do
1816 338418 : end function getRandMVN
1817 :
1818 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1819 :
1820 : ! ! Given an input Mean vector of length nd, Covariance Matrix of dimension (nd,nd), as well as a vector of integers representing
1821 : ! ! the variables (corresponding to CovMat columns) that are given
1822 : ! ! This subroutine gives out a conditional Multivariate Normal Random deviate.
1823 : ! ! random p-tivariate normal deviate, given that the first pg variables x1 are given (i.e. fixed).
1824 : ! ! For a review of Multivariate Normal distribution: Applied Multivariate Statistical Analysis, Johnson, Wichern, 1998, 4th ed.
1825 : ! !> Amir Shahmoradi, Oct 20, 2009, 9:12 PM, MTU
1826 : ! function getCondRandMVN(nd,MeanVec,CovMat,nIndIndx,IndIndx) result(CondRandMVN)
1827 : ! use Matrix_mod, only: getRegresCoef
1828 : ! implicit none
1829 : ! integer(IK), intent(in) :: nd, nIndIndx, IndIndx(nIndIndx)
1830 : ! real(RK) , intent(in) :: MeanVec(nd), CovMat(nd,nd)
1831 : ! real(RK) :: CondRandMVN(nd),
1832 : ! integer(IK) :: j, i
1833 : ! end function getCondRandMVN
1834 :
1835 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1836 :
1837 : !> \brief
1838 : !> Given the Cholesky Lower triangle and diagonals of a given covariance matrix, this function return one point uniformly
1839 : !> randomly drawn from inside of an `nd`-ellipsoid, whose `nd` elements `RandMVU(i), i=1:nd` are guaranteed
1840 : !> to be in the range:
1841 : !> ```
1842 : !> MeanVec(i) - sqrt(CovMat(i,i)) < RandMVU(i) < MeanVec(i) + sqrt(CovMat(i,i))
1843 : !> ```
1844 : !>
1845 : !> \param[in] nd : The number of dimensions of the MVU distribution.
1846 : !> \param[in] MeanVec : The mean vector of the MVU distribution.
1847 : !> \param[in] CholeskyLower : The Cholesky lower triangle of the covariance matrix of the MVU distribution.
1848 : !> \param[in] Diagonal : The Diagonal elements of the Cholesky lower triangle of the covariance matrix of the MVU distribution.
1849 : !>
1850 : !> \return
1851 : !> `RandMVU` : The randomly generated MVU vector.
1852 : !>
1853 : !> \author
1854 : !> Amir Shahmoradi, April 23, 2017, 1:36 AM, ICES, UTEXAS
1855 727694 : function getRandMVU(nd,MeanVec,CholeskyLower,Diagonal) result(RandMVU)
1856 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1857 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandMVU
1858 : #endif
1859 : implicit none
1860 : integer(IK), intent(in) :: nd
1861 : real(RK) , intent(in) :: MeanVec(nd)
1862 : real(RK) , intent(in) :: CholeskyLower(nd,nd) ! Cholesky lower triangle, calculated from the input CovMat.
1863 : real(RK) , intent(in) :: Diagonal(nd) ! Cholesky diagonal terms, calculated from the input CovMat.
1864 581419 : real(RK) :: RandMVU(nd), dummy, DummyVec(nd), sumSqDummyVec
1865 : integer(IK) :: i,j
1866 146275 : sumSqDummyVec = 0._RK
1867 435144 : do j=1,nd
1868 288869 : DummyVec(j) = getRandGaus()
1869 435144 : sumSqDummyVec = sumSqDummyVec + DummyVec(j)**2
1870 : end do
1871 146275 : call random_number(dummy)
1872 146275 : dummy = dummy**(1._RK/nd) / sqrt(sumSqDummyVec)
1873 435144 : DummyVec = DummyVec * dummy ! a uniform random point from inside of nd-sphere
1874 435144 : RandMVU = MeanVec
1875 435144 : do j = 1,nd
1876 288869 : RandMVU(j) = RandMVU(j) + Diagonal(j) * DummyVec(j)
1877 577738 : do i = j+1,nd
1878 431463 : RandMVU(i) = RandMVU(i) + CholeskyLower(i,j) * DummyVec(j)
1879 : end do
1880 : end do
1881 484690 : end function getRandMVU
1882 :
1883 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1884 :
1885 : !> \brief
1886 : !> Return `.true.` if the input `NormedPoint` (normalized with respect to the center of the target ellipsoid) is
1887 : !> within or on the boundary of the ellipsoid whose boundary is described by the representative matrix
1888 : !> \f$ \Sigma \f$ (`RepMat`), such that,
1889 : !> \f{equation}{
1890 : !> X^T ~ \Sigma^{-1} ~ X = 1 ~,
1891 : !> \f}
1892 : !> for all \f$X\f$ on the boundary.
1893 : !>
1894 : !> \param[in] nd : The number of dimensions of the ellipsoid (the size of `NormedPoint`).
1895 : !> \param[in] NormedPoint : The input point, normalized with respect to the center of the ellipsoid,
1896 : !> whose location with respect to the ellipsoid boundary is to be determined.
1897 : !> \param[in] InvRepMat : The inverse of the representative matrix of the target ellipsoid.
1898 : !>
1899 : !> \return
1900 : !> `isInsideEllipsoid` : The logical value indicating whether the input point is inside or on the boundary of the target ellipsoid.
1901 : !>
1902 : !> \remark
1903 : !> Note that the input matrix is the inverse of `RepMat`: `InvRepMat`.
1904 : !>
1905 : !> \author
1906 : !> Amir Shahmoradi, April 23, 2017, 1:36 AM, ICES, UTEXAS
1907 15 : pure function isInsideEllipsoid(nd,NormedPoint,InvRepMat)
1908 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1909 : !DEC$ ATTRIBUTES DLLEXPORT :: isInsideEllipsoid
1910 : #endif
1911 146275 : use Math_mod, only: getLogVolEllipsoid
1912 : implicit none
1913 : integer(IK), intent(in) :: nd
1914 : real(RK) , intent(in) :: NormedPoint(nd)
1915 : real(RK) , intent(in) :: InvRepMat(nd,nd)
1916 : logical :: isInsideEllipsoid
1917 45 : isInsideEllipsoid = dot_product(NormedPoint,matmul(InvRepMat,NormedPoint)) <= 1._RK
1918 15 : end function isInsideEllipsoid
1919 :
1920 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1921 :
1922 : !> \brief
1923 : !> Return the natural logarithm of the probability density function value of a point uniformly distributed within an ellipsoid,
1924 : !> whose logarithm of the square root of the determinant of its representative covariance matrix is given by `logSqrtDetCovMat`.
1925 : !>
1926 : !> \param[in] nd : The number of dimensions of the MVU distribution.
1927 : !> \param[in] logSqrtDetCovMat : The logarithm of the square root of the determinant of
1928 : !> the inverse of the representative covariance matrix of the ellipsoid.
1929 : !>
1930 : !> \return
1931 : !> `logProbMVU` : The natural logarithm of the probability density function value
1932 : !> of a point uniformly distributed within the target ellipsoid.
1933 : !>
1934 : !> \author
1935 : !> Amir Shahmoradi, April 23, 2017, 1:36 AM, ICES, UTEXAS
1936 3 : pure function getLogProbMVU(nd,logSqrtDetCovMat) result(logProbMVU)
1937 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1938 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbMVU
1939 : #endif
1940 15 : use Math_mod, only: getLogVolEllipsoid
1941 : implicit none
1942 : integer(IK), intent(in) :: nd
1943 : real(RK) , intent(in) :: logSqrtDetCovMat
1944 : real(RK) :: logProbMVU
1945 3 : logProbMVU = -getLogVolEllipsoid(nd,logSqrtDetCovMat)
1946 6 : end function getLogProbMVU
1947 :
1948 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1949 :
1950 : !> \brief
1951 : !> Return a random point on the target ellipsoid by projecting a random point uniformly distributed within the ellipsoid on its boundary.
1952 : !>
1953 : !> \param[in] nd : The number of dimensions of the ellipsoid.
1954 : !> \param[in] MeanVec : The mean vector (center) of the ellipsoid.
1955 : !> \param[in] CholeskyLower : The Cholesky lower triangle of the representative covariance matrix of the ellipsoid.
1956 : !> \param[in] Diagonal : The Diagonal elements of the Cholesky lower triangle of the representative covariance matrix of the ellipsoid.
1957 : !>
1958 : !> \return
1959 : !> `RandPointOnEllipsoid` : A random point on the target ellipsoid by projecting a random
1960 : !> point uniformly distributed within the ellipsoid on its boundary.
1961 : !>
1962 : !> \remark
1963 : !> This is algorithm is similar to [getRandMVU](@ref getrandmvu), with the only difference that
1964 : !> points are drawn randomly from the surface of the ellipsoid instead of inside of its interior.
1965 : !>
1966 : !> \remark
1967 : !> Note that the distribution of points on the surface of the ellipsoid is **NOT** uniform.
1968 : !> Regions of high curvature will have more points randomly sampled from them.
1969 : !> Generating uniform random points on arbitrary-dimension ellipsoids is not a task with trivial solution!
1970 : !>
1971 : !> \todo
1972 : !> The performance of this algorithm can be further improved.
1973 : !>
1974 : !> \author
1975 : !> Amir Shahmoradi, April 23, 2017, 1:36 AM, ICES, UTEXAS
1976 15 : function getRandPointOnEllipsoid(nd,MeanVec,CholeskyLower,Diagonal) result(RandPointOnEllipsoid)
1977 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1978 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandPointOnEllipsoid
1979 : #endif
1980 : implicit none
1981 : integer(IK), intent(in) :: nd
1982 : real(RK) , intent(in) :: MeanVec(nd)
1983 : real(RK) , intent(in) :: CholeskyLower(nd,nd) ! Cholesky lower triangle, calculated from the MVN CovMat.
1984 : real(RK) , intent(in) :: Diagonal(nd) ! Cholesky diagonal terms, calculated from the MVN CovMat.
1985 12 : real(RK) :: RandPointOnEllipsoid(nd), DummyVec(nd), sumSqDummyVec
1986 : integer(IK) :: i,j
1987 3 : sumSqDummyVec = 0._RK
1988 9 : do j=1,nd
1989 6 : DummyVec(j) = getRandGaus()
1990 9 : sumSqDummyVec = sumSqDummyVec + DummyVec(j)**2
1991 : end do
1992 9 : DummyVec = DummyVec / sqrt(sumSqDummyVec) ! DummyVec is a random point on the surface of nd-sphere.
1993 9 : RandPointOnEllipsoid = 0._RK
1994 9 : do j = 1,nd
1995 6 : RandPointOnEllipsoid(j) = RandPointOnEllipsoid(j) + Diagonal(j) * DummyVec(j)
1996 12 : do i = j+1,nd
1997 9 : RandPointOnEllipsoid(i) = RandPointOnEllipsoid(i) + CholeskyLower(i,j) * DummyVec(j)
1998 : end do
1999 : end do
2000 9 : RandPointOnEllipsoid = RandPointOnEllipsoid + MeanVec
2001 6 : end function getRandPointOnEllipsoid
2002 :
2003 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2004 :
2005 : !> \brief
2006 : !> Return the natural logarithm of the Lognormal probability density function.
2007 : !>
2008 : !> \param[in] logMean : The mean of the Lognormal distribution.
2009 : !> \param[in] inverseVariance : The inverse variance of the Lognormal distribution.
2010 : !> \param[in] logSqrtInverseVariance : The natural logarithm of the square root of the inverse variance of the Lognormal distribution.
2011 : !> \param[in] logPoint : The natural logarithm of the point at which the Lognormal PDF must be computed.
2012 : !>
2013 : !> \return
2014 : !> `logProbLogn` : The natural logarithm of the Lognormal probability density function.
2015 : !>
2016 : !> \author
2017 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2018 3894 : pure function getLogProbLognSP(logMean,inverseVariance,logSqrtInverseVariance,logPoint) result(logProbLogn)
2019 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2020 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbLognSP
2021 : #endif
2022 3 : use Constants_mod, only: LOGINVSQRT2PI
2023 : implicit none
2024 : real(RK), intent(in) :: logMean,inverseVariance,logSqrtInverseVariance,logPoint
2025 : real(RK) :: logProbLogn
2026 3894 : logProbLogn = LOGINVSQRT2PI + logSqrtInverseVariance - logPoint - 0.5_RK * inverseVariance * (logPoint-logMean)**2
2027 7788 : end function getLogProbLognSP
2028 :
2029 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2030 :
2031 : !> \brief
2032 : !> Return the natural logarithm of the Lognormal probability density function.
2033 : !>
2034 : !> \param[in] np : The size of the input vector of points represented by `LogPoint`.
2035 : !> \param[in] logMean : The mean of the Lognormal distribution.
2036 : !> \param[in] inverseVariance : The inverse variance of the Lognormal distribution.
2037 : !> \param[in] logSqrtInverseVariance : The natural logarithm of the square root of the inverse variance of the Lognormal distribution.
2038 : !> \param[in] LogPoint : The natural logarithm of the vector of points at which the Lognormal PDF must be computed.
2039 : !>
2040 : !> \return
2041 : !> `logProbLogn` : The natural logarithm of the Lognormal probability density function.
2042 : !>
2043 : !> \author
2044 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2045 36 : pure function getLogProbLognMP(np,logMean,inverseVariance,logSqrtInverseVariance,LogPoint) result(logProbLogn)
2046 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2047 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbLognMP
2048 : #endif
2049 3894 : use Constants_mod, only: LOGINVSQRT2PI
2050 : implicit none
2051 : integer(IK), intent(in) :: np
2052 : real(RK) , intent(in) :: logMean,inverseVariance,logSqrtInverseVariance,LogPoint(np)
2053 : real(RK) :: logProbLogn(np)
2054 24 : logProbLogn = LOGINVSQRT2PI + logSqrtInverseVariance - LogPoint - 0.5_RK * inverseVariance * (LogPoint-logMean)**2
2055 6 : end function getLogProbLognMP
2056 :
2057 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2058 :
2059 : !> \brief
2060 : !> Return a single-precision uniformly-distributed random real-valued number in the range `[0,1]`.
2061 : !>
2062 : !> \param[inout] idum : The input integer random seed with the `save` attribute.
2063 : !>
2064 : !> \return
2065 : !> `randRealLecuyer` : A single-precision uniformly-distributed random real-valued number in the range `[0,1]`.
2066 : !>
2067 : !> \warning
2068 : !> Do not change the value of `idum` between calls.
2069 : !>
2070 : !> \remark
2071 : !> This routine is guaranteed to random numbers with priodicity larger than `~2*10**18` random numbers.
2072 : !> For more info see Press et al. (1990) Numerical Recipes.
2073 : !>
2074 : !> \remark
2075 : !> This routine is solely kept for backward compatibility reasons.
2076 : !> The Fortran intrinsic subroutine `random_number()` is recommended to be used against this function.
2077 : !>
2078 : !> \author
2079 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2080 600 : function getRandRealLecuyer(idum) result(randRealLecuyer)
2081 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2082 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandRealLecuyer
2083 : #endif
2084 : implicit none
2085 : integer(IK), intent(inout) :: idum
2086 : integer(IK), parameter :: im1=2147483563, im2=2147483399, imm1=im1-1, ia1=40014, ia2=40692
2087 : integer(IK), parameter :: iq1=53668, iq2=52774, ir1=12211, ir2=3791, ntab=32, ndiv=1+imm1/ntab
2088 : real(RK) , parameter :: am=1._RK/real(im1,kind=RK), eps=1.2e-7_RK, rnmx=1._RK-eps
2089 : real(RK) :: randRealLecuyer
2090 : integer(IK) :: idum2,j,k,iv(ntab),iy
2091 : save :: iv, iy, idum2
2092 : data idum2/123456789/, iv/ntab*0/, iy/0/
2093 600 : if (idum <= 0) then
2094 0 : idum = max(-idum,1)
2095 0 : idum2 = idum
2096 0 : do j = ntab+8,1,-1
2097 0 : k = idum/iq1
2098 0 : idum = ia1*(idum-k*iq1)-k*ir1
2099 0 : if (idum < 0) idum = idum+im1
2100 0 : if (j <= ntab) iv(j) = idum
2101 : end do
2102 0 : iy = iv(1)
2103 : endif
2104 600 : k = idum/iq1
2105 600 : idum = ia1*(idum-k*iq1)-k*ir1
2106 600 : if (idum < 0) idum=idum+im1
2107 600 : k = idum2/iq2
2108 600 : idum2 = ia2*(idum2-k*iq2)-k*ir2
2109 600 : if (idum2 < 0) idum2=idum2+im2
2110 600 : j = 1+iy/ndiv
2111 600 : iy = iv(j)-idum2
2112 600 : iv(j) = idum
2113 600 : if(iy < 1)iy = iy+imm1
2114 600 : randRealLecuyer = min(am*iy,rnmx)
2115 606 : end function getRandRealLecuyer
2116 :
2117 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2118 :
2119 : !> \brief
2120 : !> Return an integer uniformly-distributed random integer-valued number in the range `[lowerBound , upperBound]`.
2121 : !>
2122 : !> \param[in] lowerBound : The inclusive integer lower bound of the integer uniform distribution.
2123 : !> \param[in] upperBound : The inclusive integer upper bound of the integer uniform distribution.
2124 : !> \param[inout] idum : The input integer random seed with the `save` attribute.
2125 : !>
2126 : !> \return
2127 : !> `randRealLecuyer` : A uniformly-distributed random integer-valued number in the range `[lowerBound , upperBound]`.
2128 : !>
2129 : !> \warning
2130 : !> Do not change the value of `idum` between calls.
2131 : !>
2132 : !> \remark
2133 : !> This routine is guaranteed to random numbers with priodicity larger than `~2*10**18` random numbers.
2134 : !> For more info see Press et al. (1990) Numerical Recipes.
2135 : !>
2136 : !> \remark
2137 : !> This routine is solely kept for backward compatibility reasons.
2138 : !> The [getRandInt](@ref getrandint) is recommended to be used against this routine.
2139 : !>
2140 : !> \author
2141 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2142 300 : function getRandIntLecuyer(lowerBound,upperBound,idum)
2143 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2144 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandIntLecuyer
2145 : #endif
2146 : implicit none
2147 : integer(IK), intent(in) :: lowerBound,upperBound
2148 : integer(IK), intent(inout) :: idum
2149 : integer(IK) :: getRandIntLecuyer
2150 600 : getRandIntLecuyer = lowerBound + nint( getRandRealLecuyer(idum)*real(upperBound-lowerBound,kind=RK) )
2151 900 : end function getRandIntLecuyer
2152 :
2153 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2154 :
2155 : !> \brief
2156 : !> Return an integer uniformly-distributed random integer-valued number in the range `[lowerBound , upperBound]`.
2157 : !>
2158 : !> \param[in] lowerBound : The lower bound of the integer uniform distribution.
2159 : !> \param[in] upperBound : The upper bound of the integer uniform distribution.
2160 : !>
2161 : !> \return
2162 : !> `randInt` : A uniformly-distributed random integer-valued number in the range `[lowerBound , upperBound]`.
2163 : !>
2164 : !> \author
2165 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2166 60113 : function getRandInt(lowerBound,upperBound) result(randInt)
2167 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2168 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandInt
2169 : #endif
2170 : implicit none
2171 : integer(IK), intent(in) :: lowerBound,upperBound
2172 60113 : real(RK) :: dummy
2173 : integer(IK) :: randInt
2174 60113 : call random_number(dummy)
2175 60113 : randInt = lowerBound + nint( dummy*real(upperBound-lowerBound,kind=RK) )
2176 60413 : end function getRandInt
2177 :
2178 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2179 :
2180 : !> \brief
2181 : !> Return an integer uniformly-distributed random integer-valued number in the range `[lowerBound , upperBound]` using
2182 : !> the built-in random number generator of Fortran.
2183 : !>
2184 : !> \param[in] lowerBound : The lower bound of the integer uniform distribution.
2185 : !> \param[in] upperBound : The upper bound of the integer uniform distribution.
2186 : !>
2187 : !> \return
2188 : !> `randUniform` : A uniformly-distributed random real-valued number in the range `[lowerBound , upperBound]`.
2189 : !>
2190 : !> \author
2191 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2192 300 : function getRandUniform(lowerBound,upperBound) result(randUniform)
2193 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2194 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandUniform
2195 : #endif
2196 : implicit none
2197 : real(RK), intent(in) :: lowerBound, upperBound
2198 : real(RK) :: randUniform
2199 300 : call random_number(randUniform)
2200 300 : randUniform = lowerBound + randUniform * (upperBound - lowerBound)
2201 60413 : end function getRandUniform
2202 :
2203 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2204 :
2205 : !> \brief
2206 : !> Return a Gamma-distributed random number, following the prescription in the GSL library.
2207 : !>
2208 : !> \param[in] alpha : The shape parameter of the Gamma distribution.
2209 : !>
2210 : !> \return
2211 : !> `randGamma` : A Gamma-distributed random real-valued number in the range `[0,+Infinity]`.
2212 : !>
2213 : !> \warning
2214 : !> The condition `alpha > 0` must hold, otherwise the negative value `-1` will be returned to indicate the occurrence of error.
2215 : !>
2216 : !> \author
2217 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2218 1503 : function getRandGamma(alpha) result(randGamma)
2219 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2220 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandGamma
2221 : #endif
2222 : implicit none
2223 : real(RK), intent(in) :: alpha
2224 : real(RK) :: randGamma
2225 1503 : real(RK) :: c,u,v,z
2226 1503 : if (alpha<=0._RK) then ! illegal value of alpha
2227 3 : randGamma = -1._RK
2228 3 : return
2229 : else
2230 1500 : randGamma = alpha
2231 1500 : if (randGamma<1._RK) randGamma = randGamma + 1._RK
2232 1500 : randGamma = randGamma - 0.3333333333333333_RK
2233 1500 : c = 3._RK*sqrt(randGamma)
2234 1500 : c = 1._RK / c
2235 22 : do
2236 0 : do
2237 1522 : z = getRandGaus()
2238 1522 : v = 1._RK + c*z
2239 1522 : if (v<=0._RK) cycle
2240 1522 : exit
2241 : end do
2242 1522 : v = v**3
2243 1522 : call random_number(u)
2244 1522 : if ( log(u) >= 0.5_RK * z**2 + randGamma * ( 1._RK - v + log(v) ) ) cycle
2245 1500 : randGamma = randGamma * v
2246 1522 : exit
2247 : end do
2248 1500 : if (alpha<1._RK) then
2249 0 : call random_number(u)
2250 0 : randGamma = randGamma * u**(1._RK/alpha)
2251 : end if
2252 : end if
2253 1803 : end function getRandGamma
2254 :
2255 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2256 :
2257 : !> \brief
2258 : !> Return a Gamma-distributed random number, whose shape parameter `alpha` is an integer.
2259 : !>
2260 : !> \param[in] alpha : The shape integer parameter of the Gamma distribution.
2261 : !>
2262 : !> \return
2263 : !> `randGamma` : A Gamma-distributed random real-valued number in the range `[0,+Infinity]`.
2264 : !>
2265 : !> \warning
2266 : !> The condition `alpha > 1` must hold, otherwise the negative value `-1` will be returned to indicate the occurrence of error.
2267 : !>
2268 : !> \author
2269 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2270 303 : function getRandGammaIntShape(alpha) result(randGamma)
2271 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2272 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandGammaIntShape
2273 : #endif
2274 : implicit none
2275 : integer(IK), intent(in) :: alpha
2276 : real(RK) :: randGamma
2277 303 : real(RK) :: am,e,h,s,x,y,Vector(2),Array(5)
2278 303 : if (alpha < 1) then ! illegal value of alpha
2279 3 : randGamma = -1._RK
2280 3 : return
2281 300 : elseif (alpha < 6) then
2282 300 : call random_number(Array(1:alpha))
2283 900 : x = -log(product(Array(1:alpha)))
2284 : else ! use rejection sampling
2285 0 : do
2286 0 : call random_number(Vector)
2287 0 : Vector(2) = 2._RK*Vector(2)-1._RK
2288 0 : if (dot_product(Vector,Vector) > 1._RK) cycle
2289 0 : y = Vector(2) / Vector(1)
2290 0 : am = alpha - 1
2291 0 : s = sqrt(2._RK*am + 1._RK)
2292 0 : x = s*y + am
2293 0 : if (x <= 0.0) cycle
2294 0 : e = (1._RK+y**2) * exp(am*log(x/am)-s*y)
2295 0 : call random_number(h)
2296 0 : if (h <= e) exit
2297 : end do
2298 : end if
2299 300 : randGamma = x
2300 1806 : end function getRandGammaIntShape
2301 :
2302 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2303 :
2304 : !> \brief
2305 : !> Return a random Beta-distributed variable.
2306 : !>
2307 : !> \param[in] alpha : The first shape parameter of the Beta distribution.
2308 : !> \param[in] beta : The second shape parameter of the Beta distribution.
2309 : !>
2310 : !> \return
2311 : !> `randBeta` : A Beta-distributed random real-valued number in the range `[0,1]`.
2312 : !>
2313 : !> \warning
2314 : !> The conditions `alpha > 0` and `beta > 0` must hold, otherwise the negative
2315 : !> value `-1` will be returned to indicate the occurrence of error.
2316 : !>
2317 : !> \author
2318 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2319 609 : function getRandBeta(alpha,beta) result(randBeta)
2320 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2321 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandBeta
2322 : #endif
2323 : implicit none
2324 : real(RK), intent(in) :: alpha,beta
2325 : real(RK) :: randBeta
2326 609 : real(RK) :: x
2327 609 : if ( alpha>0._RK .and. beta>0._RK ) then
2328 600 : x = getRandGamma(alpha)
2329 600 : randBeta = x / ( x + getRandGamma(beta) )
2330 : else ! illegal value of alpha or beta
2331 9 : randBeta = -1._RK
2332 : end if
2333 912 : end function getRandBeta
2334 :
2335 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2336 :
2337 : !> \brief
2338 : !> Return a random Exponential-distributed value whose inverse mean is given as input.
2339 : !>
2340 : !> \param[in] invMean : The inverse of the mean of the Exponential distribution.
2341 : !>
2342 : !> \return
2343 : !> `randExp` : An Exponential-distributed random real-valued number in the range `[0,+Infinity]` with mean `1 / invMean`.
2344 : !>
2345 : !> \warning
2346 : !> It is the user's onus to ensure `invMean > 0` on input. This condition will NOT be checked by this routine.
2347 : !>
2348 : !> \author
2349 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2350 300 : function getRandExpWithInvMean(invMean) result(randExp)
2351 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2352 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandExpWithInvMean
2353 : #endif
2354 : implicit none
2355 : real(RK), intent(in) :: invMean
2356 : real(RK) :: randExp
2357 300 : call random_number(randExp)
2358 300 : randExp = -log(randExp) * invMean
2359 909 : end function getRandExpWithInvMean
2360 :
2361 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2362 :
2363 : !> \brief
2364 : !> Return a random Exponential-distributed value whose mean is \f$\lambda = 1\f$.
2365 : !>
2366 : !> \return
2367 : !> `randExp` : A random Exponential-distributed value whose mean \f$\lambda = 1\f$.
2368 : !>
2369 : !> \author
2370 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2371 300 : function getRandExp() result(randExp)
2372 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2373 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandExp
2374 : #endif
2375 : implicit none
2376 : real(RK) :: randExp
2377 300 : call random_number(randExp)
2378 300 : randExp = -log(randExp)
2379 600 : end function getRandExp
2380 :
2381 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2382 :
2383 : !> \brief
2384 : !> Return a random correlation matrix.
2385 : !>
2386 : !> \param[in] nd : The rank of the correlation matrix.
2387 : !> \param[in] eta : The parameter that roughly represents the shape parameters of the Beta distribution.
2388 : !> The larger the value of `eta` is, the more homogeneous the correlation matrix will look.
2389 : !> In other words, set this parameter to some small number to generate strong random correlations
2390 : !> in the output random correlation matrix.
2391 : !>
2392 : !> \return
2393 : !> `RandCorMat` : A random correlation matrix.
2394 : !>
2395 : !> \warning
2396 : !> The conditions `nd > 1` and `eta > 0.0` must hold, otherwise the first element of
2397 : !> output, `getRandCorMat(1,1)`, will be set to `-1` to indicate the occurrence of an error.
2398 : !>
2399 : !> \author
2400 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2401 2781 : function getRandCorMat(nd,eta) result(RandCorMat)
2402 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2403 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandCorMat
2404 : #endif
2405 300 : use Matrix_mod, only: getCholeskyFactor
2406 : implicit none
2407 : integer(IK), intent(in) :: nd
2408 : real(RK) , intent(in) :: eta
2409 309 : real(RK) :: RandCorMat(nd,nd), dummy
2410 1545 : real(RK) :: beta,sumSqDummyVec,DummyVec(nd-1),W(nd-1),Diagonal(nd-1)
2411 : integer(IK) :: m,j,i
2412 :
2413 309 : if (nd<2_IK .or. eta<=0._RK) then ! illegal value for eta.
2414 9 : RandCorMat(1,1) = -1._RK
2415 9 : return
2416 : end if
2417 :
2418 900 : do m = 1,nd
2419 900 : RandCorMat(m,m) = 1._RK
2420 : end do
2421 300 : beta = eta + 0.5_RK*(nd-2._RK)
2422 300 : dummy = getRandBeta(beta,beta)
2423 300 : if (dummy<=0._RK .or. dummy>=1._RK) then
2424 : ! LCOV_EXCL_START
2425 : error stop
2426 : !call abortProgram( output_unit , 1 , 1 , 'Statitistics@getRandCorMat() failed. Random Beta variable out of bound: ' // num2str(dummy) )
2427 : end if
2428 : ! LCOV_EXCL_STOP
2429 300 : RandCorMat(1,2) = 2._RK * dummy - 1._RK ! for the moment, only the upper half of RandCorMat is needed, the lower half will contain cholesky lower triangle.
2430 :
2431 300 : do m = 2,nd-1
2432 0 : beta = beta - 0.5_RK
2433 0 : sumSqDummyVec = 0._RK
2434 0 : do j=1,m
2435 0 : DummyVec(j) = getRandGaus()
2436 0 : sumSqDummyVec = sumSqDummyVec + DummyVec(j)**2
2437 : end do
2438 0 : DummyVec(1:m) = DummyVec(1:m) / sqrt(sumSqDummyVec) ! DummyVec is now a uniform random point from inside of m-sphere
2439 0 : dummy = getRandBeta(0.5e0_RK*m,beta)
2440 0 : W(1:m) = sqrt(dummy) * DummyVec(1:m)
2441 0 : call getCholeskyFactor(m,RandCorMat(1:m,1:m),Diagonal(1:m))
2442 0 : if (Diagonal(1)<0._RK) then
2443 0 : error stop
2444 : !call abortProgram( output_unit , 1 , 1 , 'Statitistics@getRandCorMat()@getCholeskyFactor() failed.' )
2445 : end if
2446 0 : DummyVec(1:m) = 0._RK
2447 0 : do j = 1,m
2448 0 : DummyVec(j) = DummyVec(j) + Diagonal(j) * W(j)
2449 0 : do i = j+1,m
2450 0 : DummyVec(i) = DummyVec(i) + RandCorMat(i,j) * DummyVec(j)
2451 : end do
2452 : end do
2453 300 : RandCorMat(1:m,m+1) = DummyVec(1:m)
2454 : end do
2455 600 : do i=1,nd-1
2456 900 : RandCorMat(i+1:nd,i) = RandCorMat(i,i+1:nd)
2457 : end do
2458 618 : end function getRandCorMat
2459 :
2460 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2461 :
2462 : ! function getRandCorMat(nd,eta) ! based on the idea of LKJ (2007). But there is something wrong with this routine
2463 : ! use Matrix_mod, only: getCholeskyFactor
2464 : ! implicit none
2465 : ! !integer, intent(in) :: nd,eta
2466 : ! integer, intent(in) :: nd
2467 : ! real(RK), intent(in) :: eta
2468 : ! integer :: m,mNew,j,i
2469 : ! real(RK) :: getRandCorMat(nd,nd), dummy, failureCounter
2470 : ! real(RK) :: beta,sumSqDummyVec,DummyVec(nd-1),W(nd-1),Diagonal(nd-1)
2471 : !
2472 : ! if (nd<2 .or. eta<=0._RK) then ! illegal value for eta. set getRandCorMat=0, return
2473 : ! getRandCorMat = -1._RK
2474 : ! return
2475 : ! end if
2476 : !
2477 : ! do m = 1,nd
2478 : ! getRandCorMat(m,m) = 1._RK
2479 : ! end do
2480 : ! beta = eta + 0.5_RK*(nd-2._RK)
2481 : !
2482 : ! do
2483 : ! dummy = getRandBeta(beta,beta)
2484 : ! if (dummy>0._RK .and. dummy<1._RK) exit
2485 : ! write(*,*) "**Warning** random Beta variable out of bound.", dummy
2486 : ! write(*,*) "Something is wrong with getRandBeta()."
2487 : ! cycle
2488 : ! end do
2489 : ! getRandCorMat(1,2) = 2._RK * dummy - 1._RK ! for the moment, only the upper half of getRandCorMat is needed, the lower half will contain cholesky lower triangle.
2490 : !
2491 : ! m = 2
2492 : ! call getCholeskyFactor(m,getRandCorMat(1:m,1:m),Diagonal(1:m))
2493 : !
2494 : ! failureCounter = 0
2495 : ! onionLayer: do
2496 : !
2497 : ! beta = beta - 0.5_RK
2498 : !
2499 : ! sumSqDummyVec = 0._RK
2500 : ! do j=1,m
2501 : ! DummyVec(j) = getRandGaus()
2502 : ! sumSqDummyVec = sumSqDummyVec + DummyVec(j)**2
2503 : ! end do
2504 : ! DummyVec(1:m) = DummyVec(1:m) / sqrt(sumSqDummyVec) ! DummyVec is now a uniform random point from inside of m-sphere
2505 : !
2506 : ! mNew = m + 1
2507 : ! posDefCheck: do
2508 : !
2509 : ! do
2510 : ! dummy = getRandBeta(0.5_RK*m,beta)
2511 : ! if (dummy>0._RK .and. dummy<1._RK) exit
2512 : ! write(*,*) "**Warning** random Beta variable out of bound.", dummy
2513 : ! write(*,*) "Something is wrong with getRandBeta()."
2514 : ! read(*,*)
2515 : ! cycle
2516 : ! end do
2517 : ! W(1:m) = sqrt(dummy) * DummyVec(1:m)
2518 : !
2519 : ! getRandCorMat(1:m,mNew) = 0._RK
2520 : ! do j = 1,m
2521 : ! getRandCorMat(j,mNew) = getRandCorMat(j,mNew) + Diagonal(j) * W(j)
2522 : ! do i = j+1,m
2523 : ! getRandCorMat(i,mNew) = getRandCorMat(i,mNew) + getRandCorMat(i,j) * getRandCorMat(j,mNew)
2524 : ! end do
2525 : ! end do
2526 : !
2527 : !
2528 : ! call getCholeskyFactor(mNew,getRandCorMat(1:mNew,1:mNew),Diagonal(1:mNew)) ! Now check if the new matrix is positive-definite, then proceed with the next layer
2529 : ! if (Diagonal(1)<0._RK) then
2530 : ! failureCounter = failureCounter + 1
2531 : ! cycle posDefCheck
2532 : ! !write(*,*) "Cholesky factorization failed in getRandCorMat()."
2533 : ! !write(*,*) m
2534 : ! !write(*,*) getRandCorMat(1:m,1:m)
2535 : ! !stop
2536 : ! end if
2537 : ! exit posDefCheck
2538 : !
2539 : ! end do posDefCheck
2540 : !
2541 : ! if (mNew==nd) exit onionLayer
2542 : ! m = mNew
2543 : !
2544 : ! end do onionLayer
2545 : !
2546 : ! if (failureCounter>0) write(*,*) 'failureRatio: ', dble(failureCounter)/dble(nd-2)
2547 : ! do i=1,nd-1
2548 : ! getRandCorMat(i+1:nd,i) = getRandCorMat(i,i+1:nd)
2549 : ! end do
2550 : !
2551 : ! end function getRandCorMat
2552 :
2553 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2554 :
2555 : !> \brief
2556 : !> Returns a random correlation matrix using Monte Carlo rejection method.
2557 : !>
2558 : !> \param[in] nd : The rank of the correlation matrix.
2559 : !> \param[in] minRho : The minimum correlation coefficient to be expected in the output random correlation matrix.
2560 : !> \param[in] maxRho : The maximum correlation coefficient to be expected in the output random correlation matrix.
2561 : !>
2562 : !> \return
2563 : !> `RandCorMat` : A random correlation matrix. Only the upper half of
2564 : !> `RandCorMat` is the correlation matrix, lower half is NOT set on output.
2565 : !>
2566 : !> \warning
2567 : !> The conditions `nd >= 1` and `maxRho < minRho` must hold, otherwise, `RandCorMat(1,1) = -1._RK` will be returned.
2568 : !>
2569 : !> \remark
2570 : !> This subroutine is very slow for high matrix dimensions ( `nd >~ 10` ).
2571 : !>
2572 : !> \author
2573 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2574 2754 : function getRandCorMatRejection(nd,minRho,maxRho) result(RandCorMat)
2575 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2576 : !DEC$ ATTRIBUTES DLLEXPORT :: getRandCorMatRejection
2577 : #endif
2578 :
2579 309 : use Matrix_mod, only: isPosDef
2580 : implicit none
2581 : integer(IK), intent(in) :: nd
2582 : real(RK) , intent(in) :: minRho,maxRho
2583 918 : real(RK) :: RandCorMat(nd,nd), RhoVec(nd*(nd-1))
2584 : integer(IK) :: i,j,irho
2585 306 : if (maxRho<minRho .or. nd<1_IK) then
2586 6 : RandCorMat(1,1) = -1._RK
2587 6 : return
2588 : end if
2589 300 : if (nd==1_IK) then
2590 0 : RandCorMat = 1._RK
2591 : else
2592 0 : rejection: do
2593 300 : call random_number(RhoVec)
2594 900 : RhoVec = minRho + RhoVec*(maxRho-minRho)
2595 300 : irho = 0
2596 900 : do j=1,nd
2597 600 : RandCorMat(j,j) = 1._RK
2598 1200 : do i=1,j-1
2599 300 : irho = irho + 1
2600 900 : RandCorMat(i,j) = RhoVec(irho)
2601 : end do
2602 : end do
2603 300 : if (isPosDef(nd,RandCorMat)) exit rejection
2604 300 : cycle rejection
2605 : end do rejection
2606 : end if
2607 600 : do j=1,nd-1
2608 900 : RandCorMat(j+1:nd,j) = RandCorMat(j,j+1:nd)
2609 : end do
2610 612 : end function getRandCorMatRejection
2611 :
2612 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2613 :
2614 : !> \brief
2615 : !> Convert the upper-triangle covariance matrix to the upper-triangle correlation matrix.
2616 : !>
2617 : !> \param[in] nd : The rank of the covariance matrix.
2618 : !> \param[in] CovMatUpper : The upper-triangle covariance matrix. The lower-triangle will not be used.
2619 : !>
2620 : !> \return
2621 : !> `CorMatUpper` : An upper-triangle correlation matrix. The lower-triangle will NOT be set.
2622 : !>
2623 : !> \author
2624 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2625 3826 : pure function getCorMatUpperFromCovMatUpper(nd,CovMatUpper) result(CorMatUpper)
2626 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2627 : !DEC$ ATTRIBUTES DLLEXPORT :: getCorMatUpperFromCovMatUpper
2628 : #endif
2629 : implicit none
2630 : integer(IK) , intent(in) :: nd
2631 : real(RK) , intent(in) :: CovMatUpper(nd,nd)
2632 : real(RK) :: CorMatUpper(nd,nd)
2633 1345 : real(RK) :: InverseStdVec(nd)
2634 : integer(IK) :: i,j
2635 1345 : do j = 1, nd
2636 827 : InverseStdVec(j) = 1._RK / sqrt(CovMatUpper(j,j))
2637 2481 : do i = 1, j
2638 1963 : CorMatUpper(i,j) = CovMatUpper(i,j) * InverseStdVec(j) * InverseStdVec(i)
2639 : end do
2640 : end do
2641 824 : end function getCorMatUpperFromCovMatUpper
2642 :
2643 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2644 :
2645 : !> \brief
2646 : !> Convert the upper-triangle correlation matrix to the upper-triangle covariance matrix.
2647 : !>
2648 : !> \param[in] nd : The rank of the correlation matrix.
2649 : !> \param[in] StdVec : The input standard deviation vector of length `nd`.
2650 : !> \param[in] CorMatUpper : The upper-triangle correlation matrix. The lower-triangle will not be used.
2651 : !>
2652 : !> \return
2653 : !> `CovMatUpper` : An upper-triangle covariance matrix. The lower-triangle will NOT be set.
2654 : !>
2655 : !> \author
2656 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2657 27 : pure function getCovMatUpperFromCorMatUpper(nd,StdVec,CorMatUpper) result(CovMatUpper)
2658 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2659 : !DEC$ ATTRIBUTES DLLEXPORT :: getCovMatUpperFromCorMatUpper
2660 : #endif
2661 : implicit none
2662 : integer(IK) , intent(in) :: nd
2663 : real(RK) , intent(in) :: StdVec(nd), CorMatUpper(nd,nd)
2664 : real(RK) :: CovMatUpper(nd,nd)
2665 : integer(IK) :: i,j
2666 9 : do j=1,nd
2667 18 : do i=1,j
2668 15 : CovMatUpper(i,j) = CorMatUpper(i,j) * StdVec(j) * StdVec(i)
2669 : end do
2670 : end do
2671 521 : end function getCovMatUpperFromCorMatUpper
2672 :
2673 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2674 :
2675 : !> \brief
2676 : !> Convert the lower-triangle correlation matrix to the upper-triangle covariance matrix.
2677 : !>
2678 : !> \param[in] nd : The rank of the correlation matrix.
2679 : !> \param[in] StdVec : The input standard deviation vector of length `nd`.
2680 : !> \param[in] CorMatLower : The lower-triangle correlation matrix. The upper-triangle will not be used.
2681 : !>
2682 : !> \return
2683 : !> `CovMatUpper` : An upper-triangle covariance matrix. The lower-triangle will NOT be set.
2684 : !>
2685 : !> \author
2686 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2687 27 : pure function getCovMatUpperFromCorMatLower(nd,StdVec,CorMatLower) result(CovMatUpper)
2688 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2689 : !DEC$ ATTRIBUTES DLLEXPORT :: getCovMatUpperFromCorMatLower
2690 : #endif
2691 : implicit none
2692 : integer(IK) , intent(in) :: nd
2693 : real(RK) , intent(in) :: StdVec(nd), CorMatLower(nd,nd)
2694 : real(RK) :: CovMatUpper(nd,nd)
2695 : integer(IK) :: i,j
2696 9 : do j=1,nd
2697 6 : CovMatUpper(j,j) = StdVec(j)**2
2698 12 : do i=1,j-1
2699 9 : CovMatUpper(i,j) = CorMatLower(j,i) * StdVec(j) * StdVec(i)
2700 : end do
2701 : end do
2702 6 : end function getCovMatUpperFromCorMatLower
2703 :
2704 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2705 :
2706 : !> \brief
2707 : !> Convert the upper-triangle correlation matrix to the lower-triangle covariance matrix.
2708 : !>
2709 : !> \param[in] nd : The rank of the correlation matrix.
2710 : !> \param[in] StdVec : The input standard deviation vector of length `nd`.
2711 : !> \param[in] CorMatUpper : The upper-triangle correlation matrix. The lower-triangle will not be used.
2712 : !>
2713 : !> \return
2714 : !> `CovMatLower` : An lower-triangle covariance matrix. The upper-triangle will NOT be set.
2715 : !>
2716 : !> \author
2717 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2718 27 : pure function getCovMatLowerFromCorMatUpper(nd,StdVec,CorMatUpper) result(CovMatLower)
2719 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2720 : !DEC$ ATTRIBUTES DLLEXPORT :: getCovMatLowerFromCorMatUpper
2721 : #endif
2722 : implicit none
2723 : integer(IK) , intent(in) :: nd
2724 : real(RK) , intent(in) :: StdVec(nd), CorMatUpper(nd,nd)
2725 : real(RK) :: CovMatLower(nd,nd)
2726 : integer(IK) :: i,j
2727 9 : do j=1,nd
2728 6 : CovMatLower(j,j) = StdVec(j)**2
2729 12 : do i=1,j-1
2730 9 : CovMatLower(j,i) = CorMatUpper(i,j) * StdVec(j) * StdVec(i)
2731 : end do
2732 : end do
2733 6 : end function getCovMatLowerFromCorMatUpper
2734 :
2735 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2736 :
2737 : !> \brief
2738 : !> Convert the lower-triangle correlation matrix to the lower-triangle covariance matrix.
2739 : !>
2740 : !> \param[in] nd : The rank of the correlation matrix.
2741 : !> \param[in] StdVec : The input standard deviation vector of length `nd`.
2742 : !> \param[in] CorMatLower : The lower-triangle correlation matrix. The upper-triangle will not be used.
2743 : !>
2744 : !> \return
2745 : !> `CovMatLower` : An lower-triangle covariance matrix. The upper-triangle will NOT be set.
2746 : !>
2747 : !> \author
2748 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2749 27 : pure function getCovMatLowerFromCorMatLower(nd,StdVec,CorMatLower) result(CovMatLower)
2750 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2751 : !DEC$ ATTRIBUTES DLLEXPORT :: getCovMatLowerFromCorMatLower
2752 : #endif
2753 : implicit none
2754 : integer(IK) , intent(in) :: nd
2755 : real(RK) , intent(in) :: StdVec(nd), CorMatLower(nd,nd)
2756 : real(RK) :: CovMatLower(nd,nd)
2757 : integer(IK) :: i,j
2758 9 : do j=1,nd
2759 6 : CovMatLower(j,j) = StdVec(j)**2
2760 12 : do i=1,j-1
2761 9 : CovMatLower(j,i) = CorMatLower(j,i) * StdVec(j) * StdVec(i)
2762 : end do
2763 : end do
2764 6 : end function getCovMatLowerFromCorMatLower
2765 :
2766 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2767 :
2768 : !> \brief
2769 : !> Convert the input correlation matrix to the output covariance matrix.
2770 : !>
2771 : !> \param[in] nd : The rank of the correlation matrix.
2772 : !> \param[in] StdVec : The input standard deviation vector of length `nd`.
2773 : !> \param[in] CorMatUpper : The upper-triangle correlation matrix. The lower-triangle will not be used.
2774 : !>
2775 : !> \return
2776 : !> `CovMatFull` : The full covariance matrix.
2777 : !>
2778 : !> \author
2779 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2780 7098 : pure function getCovMatFromCorMatUpper(nd,StdVec,CorMatUpper) result(CovMatFull)
2781 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2782 : !DEC$ ATTRIBUTES DLLEXPORT :: getCovMatFromCorMatUpper
2783 : #endif
2784 : implicit none
2785 : integer(IK) , intent(in) :: nd
2786 : real(RK) , intent(in) :: StdVec(nd), CorMatUpper(nd,nd) ! only upper half needed
2787 : real(RK) :: CovMatFull(nd,nd)
2788 : integer(IK) :: i,j
2789 2562 : do j=1,nd
2790 1512 : CovMatFull(j,j) = StdVec(j)**2
2791 3024 : do i=1,j-1
2792 462 : CovMatFull(i,j) = CorMatUpper(i,j) * StdVec(j) * StdVec(i)
2793 1974 : CovMatFull(j,i) = CovMatFull(i,j)
2794 : end do
2795 : end do
2796 1053 : end function getCovMatFromCorMatUpper
2797 :
2798 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2799 :
2800 : ! !> \brief
2801 : ! !> Return the Geometric distribution PDF values for a range of trials, starting at index `1`.
2802 : ! !> If the probability of success on each trial is `successProb`, then the probability that
2803 : ! !> the `k`th trial (out of `k` trials) is the first success is `GeoLogPDF(k)`.
2804 : ! !>
2805 : ! !> \param[in] successProb : The probability of success.
2806 : ! !> \param[in] logPdfPrecision : The precision value below which the PDF is practically considered to be zero (**optional**).
2807 : ! !> \param[in] minSeqLen : The minimum length of the range of `k` values for which the PDF will be computed (**optional**).
2808 : ! !> \param[in] seqLen : The length of the range of `k` values for which the PDF will be computed (**optional**).
2809 : ! !> If provided, it will overwrite the the output sequence length as inferred from
2810 : ! !> the combination of `minSeqLen` and `logPdfPrecision`.
2811 : ! !>
2812 : ! !> \return
2813 : ! !> `GeoLogPDF` : An allocatable representing the geometric PDF over a range of `k` values, whose length is
2814 : ! !> `seqLen`, or if not provided, is determined from the values of `logPdfPrecision` and `minSeqLen`.
2815 : ! !>
2816 : ! !> \author
2817 : ! !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2818 : ! function getGeoLogPDF_old(successProb,logPdfPrecision,minSeqLen,seqLen) result(GeoLogPDF)
2819 : !#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2820 : ! !DEC$ ATTRIBUTES DLLEXPORT :: getGeoLogPDF_old
2821 : !#endif
2822 : ! use Constants_mod, only: IK, RK
2823 : ! implicit none
2824 : ! real(RK) , intent(in) :: successProb
2825 : ! real(RK) , intent(in), optional :: logPdfPrecision
2826 : ! integer(IK) , intent(in), optional :: minSeqLen
2827 : ! integer(IK) , intent(in), optional :: seqLen
2828 : ! real(RK) , allocatable :: GeoLogPDF(:)
2829 : ! real(RK) , parameter :: LOG_PDF_PRECISION = log(0.001_RK)
2830 : ! real(RK) :: logProbFailure
2831 : ! integer(IK) :: lenGeoLogPDF, i
2832 : ! logProbFailure = log(1._RK - successProb)
2833 : ! if (present(seqLen)) then
2834 : ! lenGeoLogPDF = seqLen
2835 : ! else
2836 : ! if (present(logPdfPrecision)) then
2837 : ! lenGeoLogPDF = ceiling( logPdfPrecision / logProbFailure)
2838 : ! else
2839 : ! lenGeoLogPDF = ceiling(LOG_PDF_PRECISION / logProbFailure)
2840 : ! end if
2841 : ! if (present(minSeqLen)) lenGeoLogPDF = max(minSeqLen,lenGeoLogPDF)
2842 : ! end if
2843 : ! allocate(GeoLogPDF(lenGeoLogPDF))
2844 : ! GeoLogPDF(1) = log(successProb)
2845 : ! do i = 2, lenGeoLogPDF
2846 : ! GeoLogPDF(i) = GeoLogPDF(i-1) + logProbFailure
2847 : ! end do
2848 : ! end function getGeoLogPDF_old
2849 :
2850 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2851 :
2852 : !> \brief
2853 : !> Return the Geometric distribution PDF values for the input number of trials,
2854 : !> the trials at which first success happens, and the success probability.
2855 : !>
2856 : !> \param[in] numTrial : The number of trials.
2857 : !> \param[in] SuccessStep : The vector of trials of length `numTrial` at which the first success happens.
2858 : !> \param[in] successProb : The probability of success.
2859 : !>
2860 : !> \return
2861 : !> `LogProbGeo` : An output vector of length `numTrial` representing the geometric PDF values.
2862 : !>
2863 : !> \author
2864 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2865 39 : pure function getLogProbGeo(numTrial, SuccessStep, successProb) result(LogProbGeo)
2866 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2867 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbGeo
2868 : #endif
2869 1050 : use Constants_mod, only: IK, RK
2870 : implicit none
2871 : integer(IK) , intent(in) :: numTrial
2872 : integer(IK) , intent(in) :: SuccessStep(numTrial)
2873 : real(RK) , intent(in) :: successProb
2874 : real(RK) :: LogProbGeo(numTrial)
2875 3 : real(RK) :: logProbSuccess, logProbFailure
2876 3 : logProbSuccess = log(successProb)
2877 3 : logProbFailure = log(1._RK - successProb)
2878 33 : LogProbGeo = logProbSuccess + (SuccessStep-1_IK) * logProbFailure
2879 3 : end function getLogProbGeo
2880 :
2881 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2882 :
2883 : !> \brief
2884 : !> Compute the natural logarithm of the Geometric distribution PDF of a limited range of Bernoulli trials,
2885 : !> starting at index `1` up to `maxNumTrial`. In other words, upon reaching the trial `maxNumTrial`,
2886 : !> the Bernoulli trials count restart from index `1`. This Cyclic Geometric distribution is
2887 : !> particularly useful in the parallelization studies of Monte Carlo simulation.
2888 : !>
2889 : !> \param[in] successProb : The probability of success.
2890 : !> \param[in] maxNumTrial : The maximum number of trails possible.
2891 : !> After `maxNumTrial` tries, the Geometric distribution restarts from index `1`.
2892 : !> \param[in] numTrial : The length of the array `SuccessStep`. Note that `numTrial < maxNumTrial` must hold.
2893 : !> \param[in] SuccessStep : A vector of length `(1:numTrial)` of integers that represent
2894 : !> the steps at which the Bernoulli successes occur.
2895 : !>
2896 : !> \return
2897 : !> `LogProbGeoCyclic` : A real-valued vector of length `(1:numTrial)` whose values represent the probabilities
2898 : !> of having Bernoulli successes at the corresponding SuccessStep values.
2899 : !>
2900 : !> \warning
2901 : !> Any value of SuccessStep must be an integer numbers between `1` and `maxNumTrial`.
2902 : !> The onus is on the user to ensure this condition holds.
2903 : !>
2904 : !> \author
2905 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2906 586978 : pure function getLogProbGeoCyclic(successProb,maxNumTrial,numTrial,SuccessStep) result(LogProbGeoCyclic)
2907 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2908 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProbGeoCyclic
2909 : #endif
2910 3 : use Constants_mod, only: IK, RK, NEGLOGINF_RK
2911 : implicit none
2912 : real(RK) , intent(in) :: successProb
2913 : integer(IK) , intent(in) :: maxNumTrial
2914 : integer(IK) , intent(in) :: numTrial
2915 : integer(IK) , intent(in) :: SuccessStep(numTrial)
2916 : real(RK) :: LogProbGeoCyclic(numTrial)
2917 75529 : real(RK) :: failureProb, logProbSuccess, logProbFailure, logDenominator, exponentiation
2918 75529 : if (successProb>0._RK .and. successProb<1._RK) then ! tolerate log(0)
2919 75529 : failureProb = 1._RK - successProb
2920 75529 : logProbSuccess = log(successProb)
2921 75529 : logProbFailure = log(failureProb)
2922 75529 : exponentiation = maxNumTrial * logProbFailure
2923 75529 : if (exponentiation<NEGLOGINF_RK) then ! tolerate underflow
2924 0 : logDenominator = 0._RK
2925 : else
2926 75529 : exponentiation = exp(exponentiation)
2927 75529 : if (exponentiation<1._RK) then ! tolerate log(0)
2928 75529 : logDenominator = log(1._RK-exponentiation)
2929 : else
2930 0 : logDenominator = NEGLOGINF_RK
2931 : end if
2932 : end if
2933 435920 : LogProbGeoCyclic = logProbSuccess + (SuccessStep-1) * logProbFailure - logDenominator
2934 0 : elseif (successProb==0._RK) then
2935 0 : LogProbGeoCyclic = -log(real(maxNumTrial,kind=RK))
2936 0 : elseif (successProb==1._RK) then
2937 0 : LogProbGeoCyclic(1) = 0._RK
2938 0 : LogProbGeoCyclic(2:numTrial) = NEGLOGINF_RK
2939 : else
2940 0 : LogProbGeoCyclic = NEGLOGINF_RK
2941 : end if
2942 75529 : end function getLogProbGeoCyclic
2943 :
2944 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2945 :
2946 : !> \brief
2947 : !> Return the standard normal distribution PDF value.
2948 : !>
2949 : !> \param[in] z : The input value at which the PDF will be computed.
2950 : !>
2951 : !> \return
2952 : !> `snormPDF` : The standard normal distribution PDF value.
2953 : !>
2954 : !> \author
2955 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2956 3 : function getSNormPDF(z) result(snormPDF)
2957 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2958 : !DEC$ ATTRIBUTES DLLEXPORT :: getSNormPDF
2959 : #endif
2960 75529 : use Constants_mod, only: INVSQRT2PI
2961 : implicit none
2962 : real(RK), intent(in) :: z
2963 : real(RK) :: snormPDF
2964 3 : snormPDF = INVSQRT2PI * exp( -0.5_RK*z**2 )
2965 6 : end function getSNormPDF
2966 :
2967 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2968 :
2969 : !> \brief
2970 : !> Return the non-standard normal distribution PDF value.
2971 : !>
2972 : !> \param[in] avg : The mean of the Normal distribution.
2973 : !> \param[in] std : The standard deviation of the Normal distribution.
2974 : !> \param[in] var : The variance of the Normal distribution.
2975 : !> \param[in] x : The point at which the PDF will be computed.
2976 : !>
2977 : !> \return
2978 : !> `normPDF` : The normal distribution PDF value at the given input point.
2979 : !>
2980 : !> \author
2981 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
2982 3 : function getNormPDF(avg,std,var,x) result(normPDF)
2983 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
2984 : !DEC$ ATTRIBUTES DLLEXPORT :: getNormPDF
2985 : #endif
2986 3 : use Constants_mod, only: INVSQRT2PI
2987 : implicit none
2988 : real(RK), intent(in) :: avg,std,var,x
2989 : real(RK) :: normPDF
2990 3 : normPDF = INVSQRT2PI * exp( -(x-avg)**2/(2._RK*var) ) / std
2991 6 : end function getNormPDF
2992 :
2993 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2994 :
2995 : !> \brief
2996 : !> Return the non-standard normal distribution Cumulative Probability Density function (CDF) value.
2997 : !>
2998 : !> \param[in] avg : The mean of the Normal distribution.
2999 : !> \param[in] std : The standard deviation of the Normal distribution.
3000 : !> \param[in] x : The point at which the PDF will be computed.
3001 : !>
3002 : !> \return
3003 : !> `normCDF` : The normal distribution CDF value at the given input point.
3004 : !>
3005 : !> \author
3006 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
3007 3 : pure function getNormCDF(avg,std,x) result(normCDF)
3008 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3009 : !DEC$ ATTRIBUTES DLLEXPORT :: getNormCDF
3010 : #endif
3011 3 : use Constants_mod, only: RK, SQRT2
3012 : implicit none
3013 : real(RK), intent(in) :: avg,std,x
3014 : real(RK) :: normCDF
3015 3 : normCDF = 0.5_RK * ( 1._RK + erf( real((x-avg)/(SQRT2*std),kind=RK) ) )
3016 6 : end function getNormCDF
3017 :
3018 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3019 :
3020 : !> \brief
3021 : !> Return the standard normal distribution Cumulative Probability Density function (CDF) value.
3022 : !>
3023 : !> \param[in] x : The point at which the PDF will be computed.
3024 : !>
3025 : !> \return
3026 : !> `normCDF` : The normal distribution CDF value at the given input point.
3027 : !>
3028 : !> \remark
3029 : !> This procedure performs all calculations in `real32` real kind. If 64-bit accuracy matters more than performance,
3030 : !> then use the [getSNormCDF_DPR](@ref getsnormcdf_dpr) for a more-accurate double-precision but slower results.
3031 : !>
3032 : !> \author
3033 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
3034 3 : pure function getSNormCDF_SPR(x) result(normCDF)
3035 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3036 : !DEC$ ATTRIBUTES DLLEXPORT :: getSNormCDF_SPR
3037 : #endif
3038 : use iso_fortran_env, only: RK => real32
3039 : implicit none
3040 : real(RK) , intent(in) :: x
3041 : real(RK) :: normCDF
3042 : real(RK), parameter :: INVSQRT2 = 1._RK / sqrt(2._RK)
3043 3 : normCDF = 0.5_RK * ( 1._RK + erf( real(x*INVSQRT2,kind=RK) ) )
3044 3 : end function getSNormCDF_SPR
3045 :
3046 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3047 :
3048 : !> \brief
3049 : !> Return the standard normal distribution Cumulative Probability Density function (CDF) value.
3050 : !>
3051 : !> \param[in] x : The point at which the PDF will be computed.
3052 : !>
3053 : !> \return
3054 : !> `normCDF` : The normal distribution CDF value at the given input point.
3055 : !>
3056 : !> \remark
3057 : !> This procedure performs all calculations in `DPR` (`real64`) real kind. If performance matters more than 64-bit accuracy,
3058 : !> then use the [getSNormCDF_SPR](@ref getsnormcdf_spr) for a faster, but less-accurate single-precision results.
3059 : !>
3060 : !> \author
3061 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
3062 153 : pure function getSNormCDF_DPR(x) result(normCDF)
3063 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3064 : !DEC$ ATTRIBUTES DLLEXPORT :: getSNormCDF_DPR
3065 : #endif
3066 : use iso_fortran_env, only: RK => real64
3067 : implicit none
3068 : real(RK), intent(in) :: x
3069 : real(RK) :: normCDF
3070 : real(RK), parameter :: INVSQRT2 = 1._RK / sqrt(2._RK)
3071 153 : normCDF = 0.5_RK * ( 1._RK + erf( real(x*INVSQRT2,kind=RK) ) )
3072 3 : end function getSNormCDF_DPR
3073 :
3074 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3075 :
3076 : !> \brief
3077 : !> Return the Beta distribution Cumulative Probability Density function (CDF) value.
3078 : !>
3079 : !> \param[in] alpha : The first shape parameter of the Beta distribution.
3080 : !> \param[in] beta : The second shape parameter of the Beta distribution.
3081 : !> \param[in] x : The point at which the CDF will be computed.
3082 : !>
3083 : !> \return
3084 : !> `betaCDF` : The Beta distribution CDF value at the given input point.
3085 : !>
3086 : !> \warning
3087 : !> If `x` is not in the range `[0,1]`, a negative value for `betaCDF` will be returned to indicate the occurrence of error.
3088 : !>
3089 : !> \warning
3090 : !> The onus is on the user to ensure that the input (`alpha`, `beta`) shape parameters are positive.
3091 : !>
3092 : !> \remark
3093 : !> This procedure performs all calculations in `real32` real kind. If 64-bit accuracy matters more than performance,
3094 : !> then use the [getBetaCDF_DPR](@ref getbetacdf_dpr) for a more-accurate double-precision but slower results.
3095 : !>
3096 : !> \todo
3097 : !> The efficiency of this code can be improved by making `x` a vector on input.
3098 : !>
3099 : !> \author
3100 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
3101 9 : function getBetaCDF_SPR(alpha,beta,x) result(betaCDF)
3102 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3103 : !DEC$ ATTRIBUTES DLLEXPORT :: getBetaCDF_SPR
3104 : #endif
3105 : use iso_fortran_env, only: RK => real32
3106 : implicit none
3107 : real(RK), intent(in) :: alpha, beta, x
3108 9 : real(RK) :: bt
3109 : real(RK) :: betaCDF
3110 9 : if (x < 0._RK .or. x > 1._RK) then
3111 6 : betaCDF = -1._RK
3112 6 : return
3113 : end if
3114 3 : if (x == 0._RK .or. x == 1._RK) then
3115 0 : bt = 0._RK
3116 : else
3117 : bt = exp( log_gamma(real(alpha+beta,kind=RK)) &
3118 : - log_gamma(real(alpha,kind=RK)) - log_gamma(real(beta,kind=RK)) &
3119 3 : + alpha*log(x) + beta*log(1._RK-x) )
3120 : end if
3121 3 : if ( x < (alpha+1._RK) / (alpha+beta+2._RK) ) then
3122 0 : betaCDF = bt * getBetaContinuedFraction(alpha,beta,x) / alpha
3123 : else
3124 3 : betaCDF = 1._RK - bt * getBetaContinuedFraction(beta,alpha,1._RK-x) / beta
3125 : end if
3126 162 : end function getBetaCDF_SPR
3127 :
3128 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3129 :
3130 : !> \brief
3131 : !> Return the Beta distribution Cumulative Probability Density function (CDF) value.
3132 : !>
3133 : !> \param[in] alpha : The first shape parameter of the Beta distribution.
3134 : !> \param[in] beta : The second shape parameter of the Beta distribution.
3135 : !> \param[in] x : The point at which the CDF will be computed.
3136 : !>
3137 : !> \return
3138 : !> `betaCDF` : The Beta distribution CDF value at the given input point.
3139 : !>
3140 : !> \warning
3141 : !> If `x` is not in the range `[0,1]`, a negative value for `betaCDF` will be returned to indicate the occurrence of error.
3142 : !>
3143 : !> \warning
3144 : !> The onus is on the user to ensure that the input (`alpha`, `beta`) shape parameters are positive.
3145 : !>
3146 : !> \remark
3147 : !> This procedure performs all calculations in `DPR` (`real64`) real kind. If performance matters more than 64-bit accuracy,
3148 : !> then use the [getBetaCDF_SPR](@ref getbetacdf_spr) for a faster, but less-accurate single-precision results.
3149 : !>
3150 : !> \todo
3151 : !> The efficiency of this code can be improved by making `x` a vector on input.
3152 : !>
3153 : !> \author
3154 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
3155 15 : function getBetaCDF_DPR(alpha,beta,x) result(betaCDF)
3156 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3157 : !DEC$ ATTRIBUTES DLLEXPORT :: getBetaCDF_DPR
3158 : #endif
3159 : use iso_fortran_env, only: RK => real64
3160 : implicit none
3161 : real(RK), intent(in) :: alpha, beta, x
3162 15 : real(RK) :: bt
3163 : real(RK) :: betaCDF
3164 15 : if (x < 0._RK .or. x > 1._RK) then
3165 6 : betaCDF = -1._RK
3166 6 : return
3167 : end if
3168 9 : if (x == 0._RK .or. x == 1._RK) then
3169 0 : bt = 0._RK
3170 : else
3171 : bt = exp( log_gamma(real(alpha+beta,kind=RK)) &
3172 : - log_gamma(real(alpha,kind=RK)) - log_gamma(real(beta,kind=RK)) &
3173 9 : + alpha*log(x) + beta*log(1._RK-x) )
3174 : end if
3175 9 : if ( x < (alpha+1._RK) / (alpha+beta+2._RK) ) then
3176 6 : betaCDF = bt * getBetaContinuedFraction(alpha,beta,x) / alpha
3177 : else
3178 3 : betaCDF = 1._RK - bt * getBetaContinuedFraction(beta,alpha,1._RK-x) / beta
3179 : end if
3180 24 : end function getBetaCDF_DPR
3181 :
3182 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3183 :
3184 : !> \brief
3185 : !> Return the Beta Continued Fraction (BCF).
3186 : !>
3187 : !> \param[in] alpha : The first shape parameter of the Beta distribution.
3188 : !> \param[in] beta : The second shape parameter of the Beta distribution.
3189 : !> \param[in] x : The point at which the BCF will be computed.
3190 : !>
3191 : !> \return
3192 : !> `betaCDF` : The BCF value at the given input point.
3193 : !>
3194 : !> \author
3195 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
3196 3 : function getBetaContinuedFraction_SPR(alpha,beta,x) result(betaContinuedFraction)
3197 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3198 : !DEC$ ATTRIBUTES DLLEXPORT :: getBetaContinuedFraction_SPR
3199 : #endif
3200 : use iso_fortran_env, only: RK => real32
3201 : implicit none
3202 : real(RK) , intent(in) :: alpha,beta,x
3203 : real(RK) , parameter :: eps = epsilon(x), fpmin = tiny(x)/eps
3204 : integer(IK), parameter :: maxit = 100_IK
3205 3 : real(RK) :: aa,c,d,del,qab,qam,qap
3206 : real(RK) :: betaContinuedFraction
3207 : integer(IK) :: m,m2
3208 3 : qab = alpha+beta
3209 3 : qap = alpha+1._RK
3210 3 : qam = alpha-1._RK
3211 3 : c = 1._RK
3212 3 : d = 1._RK-qab*x/qap
3213 0 : if (abs(d) < fpmin) d = fpmin
3214 3 : d = 1._RK/d
3215 3 : betaContinuedFraction = d
3216 6 : do m = 1,maxit
3217 6 : m2 = 2*m
3218 6 : aa = m*(beta-m)*x/((qam+m2)*(alpha+m2))
3219 6 : d = 1._RK+aa*d
3220 6 : if (abs(d) < fpmin) d = fpmin
3221 6 : c = 1._RK+aa/c
3222 6 : if (abs(c) < fpmin) c = fpmin
3223 6 : d = 1._RK/d
3224 6 : betaContinuedFraction = betaContinuedFraction*d*c
3225 6 : aa = -(alpha+m)*(qab+m)*x/((alpha+m2)*(qap+m2))
3226 6 : d = 1._RK+aa*d
3227 6 : if (abs(d) < fpmin) d = fpmin
3228 6 : c = 1._RK+aa/c
3229 6 : if (abs(c) < fpmin) c = fpmin
3230 6 : d = 1._RK/d
3231 6 : del = d*c
3232 6 : betaContinuedFraction = betaContinuedFraction*del
3233 6 : if (abs(del-1._RK) <= eps) exit
3234 : end do
3235 3 : if (m > maxit) then
3236 : ! LCOV_EXCL_START
3237 : error stop
3238 : !call abortProgram( output_unit , 1 , 1 , &
3239 : !'Statitistics@getBetaContinuedFraction_SPR() failed: alpha or beta too big, or maxit too small.' )
3240 : end if
3241 : ! LCOV_EXCL_STOP
3242 18 : end function getBetaContinuedFraction_SPR
3243 :
3244 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3245 :
3246 : !> \brief
3247 : !> Return the Beta Continued Fraction (BCF).
3248 : !>
3249 : !> \param[in] alpha : The first shape parameter of the Beta distribution.
3250 : !> \param[in] beta : The second shape parameter of the Beta distribution.
3251 : !> \param[in] x : The point at which the BCF will be computed.
3252 : !>
3253 : !> \return
3254 : !> `betaCDF` : The BCF value at the given input point.
3255 : !>
3256 : !> \author
3257 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
3258 9 : function getBetaContinuedFraction_DPR(alpha,beta,x) result(betaContinuedFraction)
3259 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3260 : !DEC$ ATTRIBUTES DLLEXPORT :: getBetaContinuedFraction_DPR
3261 : #endif
3262 : use iso_fortran_env, only: RK => real64
3263 : implicit none
3264 : real(RK) , intent(in) :: alpha,beta,x
3265 : real(RK) , parameter :: eps = epsilon(x), fpmin = tiny(x)/eps
3266 : integer(IK), parameter :: maxit = 100_IK
3267 9 : real(RK) :: aa,c,d,del,qab,qam,qap
3268 : real(RK) :: betaContinuedFraction
3269 : integer(IK) :: m,m2
3270 9 : qab = alpha+beta
3271 9 : qap = alpha+1._RK
3272 9 : qam = alpha-1._RK
3273 9 : c = 1._RK
3274 9 : d = 1._RK-qab*x/qap
3275 0 : if (abs(d) < fpmin) d = fpmin
3276 9 : d = 1._RK/d
3277 9 : betaContinuedFraction = d
3278 90 : do m = 1,maxit
3279 90 : m2 = 2*m
3280 90 : aa = m*(beta-m)*x/((qam+m2)*(alpha+m2))
3281 90 : d = 1._RK+aa*d
3282 90 : if (abs(d) < fpmin) d = fpmin
3283 90 : c = 1._RK+aa/c
3284 90 : if (abs(c) < fpmin) c = fpmin
3285 90 : d = 1._RK/d
3286 90 : betaContinuedFraction = betaContinuedFraction*d*c
3287 90 : aa = -(alpha+m)*(qab+m)*x/((alpha+m2)*(qap+m2))
3288 90 : d = 1._RK+aa*d
3289 90 : if (abs(d) < fpmin) d = fpmin
3290 90 : c = 1._RK+aa/c
3291 90 : if (abs(c) < fpmin) c = fpmin
3292 90 : d = 1._RK/d
3293 90 : del = d*c
3294 90 : betaContinuedFraction = betaContinuedFraction*del
3295 90 : if (abs(del-1._RK) <= eps) exit
3296 : end do
3297 9 : if (m > maxit) then
3298 : ! LCOV_EXCL_START
3299 : error stop
3300 : !call abortProgram( output_unit , 1 , 1 , &
3301 : !'Statitistics@getBetaContinuedFraction_DPR() failed: alpha or beta too big, or maxit too small.' )
3302 : end if
3303 : ! LCOV_EXCL_STOP
3304 12 : end function getBetaContinuedFraction_DPR
3305 :
3306 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3307 :
3308 : !> \brief
3309 : !> Return the one-sample Kolmogorov–Smirnov (KS) test results for the null hypothesis that the data
3310 : !> in vector `Point` comes from a distribution whose CDF is specified by the input `getCDF()` function.
3311 : !>
3312 : !> \param[in] np : The number of points in the input vector `Point`.
3313 : !> \param[inout] Point : The sample. On return, this array will be sorted in Ascending order.
3314 : !> \param[in] getCDF : The function returning the Cumulative Distribution Function (CDF) of the sample.
3315 : !> \param[out] statKS : The KS test statistic.
3316 : !> \param[out] probKS : The `p`-value of the test, returned as a scalar value in the range `[0,1]`.
3317 : !> The output `probKS` is the probability of observing a test statistic as extreme as,
3318 : !> or more extreme than, the observed value under the null hypothesis.
3319 : !> Small values of `probKS` cast doubt on the validity of the null hypothesis.
3320 : !> \param[out] Err : An object of class [Err_type](@ref err_mod::err_type) containing information
3321 : !> about the occurrence of any error during the KS test computation.
3322 : !>
3323 : !> \author
3324 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
3325 3 : subroutine doKS1(np,Point,getCDF,statKS,probKS,Err)
3326 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3327 : !DEC$ ATTRIBUTES DLLEXPORT :: doKS1
3328 : #endif
3329 9 : use Sort_mod, only : sortAscending
3330 : use Err_mod, only: Err_type
3331 :
3332 : implicit none
3333 :
3334 : integer(IK) , intent(in) :: np
3335 : real(RK) , intent(out) :: statKS,probKS
3336 : real(RK) , intent(inout) :: Point(np)
3337 : type(Err_type) , intent(out) :: Err
3338 :
3339 : character(*) , parameter :: PROCEDURE_NAME = MODULE_NAME//"@doKS1"
3340 3 : real(RK) :: npSqrt
3341 3 : real(RK) :: cdf,cdfObserved,dt,frac
3342 : integer(IK) :: j
3343 :
3344 : interface
3345 : function getCDF(x)
3346 : use Constants_mod, only: RK
3347 : real(RK), intent(in) :: x
3348 : real(RK) :: getCDF
3349 : end function getCDF
3350 : end interface
3351 :
3352 3 : call sortAscending(np,Point,Err)
3353 3 : if (Err%occurred) then
3354 : ! LCOV_EXCL_START
3355 : Err%msg = PROCEDURE_NAME//Err%msg
3356 : return
3357 : end if
3358 : ! LCOV_EXCL_STOP
3359 :
3360 3 : statKS = 0._RK
3361 3 : cdfObserved = 0._RK
3362 3 : npSqrt = np
3363 153 : do j = 1,np
3364 150 : frac = j/npSqrt
3365 150 : cdf = getCDF(Point(j))
3366 150 : dt = max( abs(cdfObserved-cdf) , abs(frac-cdf) )
3367 150 : if( dt > statKS ) statKS = dt
3368 153 : cdfObserved = frac
3369 : end do
3370 3 : npSqrt = sqrt(npSqrt)
3371 3 : probKS = getProbKS( (npSqrt+0.12_RK+0.11_RK/npSqrt)*statKS )
3372 :
3373 3 : end subroutine doKS1
3374 :
3375 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3376 :
3377 : !> \brief
3378 : !> Return the one-sample Kolmogorov–Smirnov (KS) test results for the assumption that the
3379 : !> points originate from a uniform distribution in [0,1]. So, all Point must be in [0,1] on input.
3380 : !>
3381 : !> \param[in] np : The number of points in the input vector `Point`.
3382 : !> \param[inout] Point : The sample. On return, this array will be sorted in Ascending order.
3383 : !> \param[out] statKS : The KS test statistic.
3384 : !> \param[out] probKS : The `p`-value of the test, returned as a scalar value in the range `[0,1]`.
3385 : !> The output `probKS` is the probability of observing a test statistic as extreme as,
3386 : !> or more extreme than, the observed value under the null hypothesis.
3387 : !> Small values of `probKS` cast doubt on the validity of the null hypothesis.
3388 : !> \param[out] Err : An object of class [Err_type](@ref err_mod::err_type) containing information
3389 : !> about the occurrence of any error during the KS test computation.
3390 : !>
3391 : !> \author
3392 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
3393 3 : pure subroutine doUniformKS1(np,Point,statKS,probKS,Err)
3394 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3395 : !DEC$ ATTRIBUTES DLLEXPORT :: doUniformKS1
3396 : #endif
3397 3 : use Sort_mod, only : sortAscending
3398 : use Err_mod, only: Err_type
3399 :
3400 : implicit none
3401 :
3402 : integer(IK) , intent(in) :: np
3403 : real(RK) , intent(out) :: statKS,probKS
3404 : real(RK) , intent(inout) :: Point(np)
3405 : type(Err_type) , intent(out) :: Err
3406 :
3407 : character(*) , parameter :: PROCEDURE_NAME = MODULE_NAME//"@doUniformKS1"
3408 3 : real(RK) :: npSqrt
3409 3 : real(RK) :: cdf,cdfObserved,dt,frac
3410 : integer(IK) :: j
3411 :
3412 3 : call sortAscending(np,Point,Err)
3413 3 : if (Err%occurred) then
3414 : ! LCOV_EXCL_START
3415 : Err%msg = PROCEDURE_NAME//Err%msg
3416 : return
3417 : ! LCOV_EXCL_STOP
3418 : end if
3419 :
3420 3 : statKS = 0._RK
3421 3 : cdfObserved = 0._RK
3422 3 : npSqrt = np
3423 153 : do j = 1,np
3424 150 : frac = j/npSqrt
3425 150 : cdf = Point(j)
3426 150 : dt = max( abs(cdfObserved-cdf) , abs(frac-cdf) )
3427 150 : if( dt > statKS ) statKS = dt
3428 153 : cdfObserved = frac
3429 : end do
3430 3 : npSqrt = sqrt(npSqrt)
3431 3 : probKS = getProbKS( (npSqrt+0.12_RK+0.11_RK/npSqrt)*statKS )
3432 3 : end subroutine doUniformKS1
3433 :
3434 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3435 :
3436 : !> \brief
3437 : !> Return the two-sample Kolmogorov–Smirnov (KS) test results under the assumption that the
3438 : !> points originate from the same distribution. It is assumed that the two input arrays are sorted ascending on input.
3439 : !>
3440 : !> \param[in] np1 : The number of points in the input vector `SortedPoint1`.
3441 : !> \param[in] np2 : The number of points in the input vector `SortedPoint2`.
3442 : !> \param[inout] SortedPoint1 : The first input sorted sample. On input, it must be sorted in ascending-order.
3443 : !> \param[inout] SortedPoint2 : The second input sorted sample. On input, it must be sorted in ascending-order.
3444 : !> \param[out] statKS : The KS test statistic.
3445 : !> \param[out] probKS : The `p`-value of the test, returned as a scalar value in the range `[0,1]`.
3446 : !> The output `probKS` is the probability of observing a test statistic as extreme as,
3447 : !> or more extreme than, the observed value under the null hypothesis.
3448 : !> Small values of `probKS` cast doubt on the validity of the null hypothesis.
3449 : !>
3450 : !> \author
3451 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
3452 315 : subroutine doSortedKS2(np1,np2,SortedPoint1,SortedPoint2,statKS,probKS)
3453 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3454 : !DEC$ ATTRIBUTES DLLEXPORT :: doSortedKS2
3455 : #endif
3456 : integer(IK) , intent(in) :: np1, np2
3457 : real(RK) , intent(in) :: SortedPoint1(np1), SortedPoint2(np2)
3458 : real(RK) , intent(out) :: statKS, probKS
3459 315 : real(RK) :: dummy1,dummy2,dt,np1_RK,np2_RK,npEffective,cdf1,cdf2
3460 : integer(IK) :: j1,j2
3461 315 : np1_RK = np1
3462 315 : np2_RK = np2
3463 315 : j1 = 1_IK
3464 315 : j2 = 1_IK
3465 315 : cdf1 = 0._RK
3466 315 : cdf2 = 0._RK
3467 315 : statKS = 0._RK
3468 986343 : do
3469 986658 : if ( j1<=np1 .and. j2<=np2 )then
3470 986343 : dummy1 = SortedPoint1(j1)
3471 986343 : dummy2 = SortedPoint2(j2)
3472 986343 : if( dummy1 <= dummy2 ) then
3473 494247 : cdf1 = j1 / np1_RK
3474 494247 : j1 = j1 + 1
3475 : endif
3476 986343 : if( dummy2 <= dummy1 ) then
3477 494256 : cdf2 = j2 / np2_RK
3478 494256 : j2 = j2 + 1_IK
3479 : endif
3480 986343 : dt = abs(cdf2-cdf1)
3481 986343 : if (dt>statKS) statKS = dt
3482 986343 : cycle
3483 : end if
3484 315 : exit
3485 : end do
3486 315 : npEffective = sqrt( np1_RK * np2_RK / ( np1_RK + np2_RK ) )
3487 315 : probKS = getProbKS( ( npEffective + 0.12_RK + 0.11_RK / npEffective ) * statKS )
3488 3 : end subroutine doSortedKS2
3489 :
3490 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3491 :
3492 : !> \brief
3493 : !> Return the Kolmogorov–Smirnov (KS) probability.
3494 : !>
3495 : !> \author
3496 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
3497 321 : pure function getProbKS(lambda) result(probKS)
3498 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3499 : !DEC$ ATTRIBUTES DLLEXPORT :: getProbKS
3500 : #endif
3501 : implicit none
3502 : real(RK) , intent(in) :: lambda
3503 : real(RK) , parameter :: EPS1 = 0.001_RK, EPS2 = 1.e-8_RK
3504 : integer(IK), parameter :: NITER = 100
3505 : integer(IK) :: j
3506 321 : real(RK) :: a2,fac,term,termbf
3507 : real(RK) :: probKS
3508 321 : a2 = -2._RK*lambda**2
3509 321 : fac = 2._RK
3510 321 : probKS = 0._RK
3511 321 : termbf = 0._RK
3512 870 : do j = 1, NITER
3513 870 : term = fac*exp(a2*j**2)
3514 870 : probKS = probKS+term
3515 870 : if (abs(term) <= EPS1*termbf .or. abs(term) <= EPS2*probKS) return
3516 549 : fac = -fac
3517 549 : termbf = abs(term)
3518 : end do
3519 0 : probKS = 1._RK
3520 636 : end function getProbKS
3521 :
3522 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3523 :
3524 : !> \brief
3525 : !> Returns the uniform CDF on support [0,1). This is rather redundant, aint it? but sometimes, needed.
3526 : !>
3527 : !> \param[in] x : The point at which the CDF must be computed.
3528 : !>
3529 : !> \author
3530 : !> Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, ICES, The University of Texas at Austin.
3531 3 : pure function getUniformCDF(x)
3532 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3533 : !DEC$ ATTRIBUTES DLLEXPORT :: getUniformCDF
3534 : #endif
3535 : implicit none
3536 : real(RK), intent(in) :: x
3537 : real(RK) :: getUniformCDF
3538 3 : getUniformCDF = x
3539 321 : end function getUniformCDF
3540 :
3541 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3542 :
3543 : !> \brief
3544 : !> Return the 1-D histogram (Density plot) of the input vector `X`.
3545 : !> The number of bins in the `X` range (`nxbin`) is determined by the user.
3546 : !> The range of `X`, `[xmin, xmax]`, should be also given by the user.
3547 : !> The program returns two arrays of `Xbin` and `Density(x)` as output.
3548 : !>
3549 : !> \param[in] method : The method by which the hist count should be returned:
3550 : !> + `"pdf"` : Return the probability density function histogram.
3551 : !> + `"count"` : Return the count histogram.
3552 : !> \param[in] xmin : The minimum of the histogram binning.
3553 : !> \param[in] xmax : The maximum of the histogram binning.
3554 : !> \param[in] nxbin : The number of histogram bins.
3555 : !> \param[in] np : The length of input vector `X`.
3556 : !> \param[in] X : The vector of length `nxbin` of values to be binned.
3557 : !> \param[out] Xbin : The vector of length `nxbin` of values representing the bin left corners.
3558 : !> \param[out] Density : The vector of length `nxbin` of values representing the densities in each bin.
3559 : !> \param[out] errorOccurred : The logical output flag indicating whether error has occurred.
3560 : !>
3561 : !> \author
3562 : !> Amir Shahmoradi, Sep 1, 2017, 12:30 AM, ICES, UT Austin
3563 9 : subroutine getHist1D(method, xmin, xmax, nxbin, np, X, Xbin, Density, errorOccurred)
3564 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3565 : !DEC$ ATTRIBUTES DLLEXPORT :: getHist1D
3566 : #endif
3567 : implicit none
3568 : character(*), intent(in) :: method
3569 : integer(IK) , intent(in) :: np,nxbin
3570 : real(RK) , intent(in) :: xmin,xmax
3571 : real(RK) , intent(in) :: X(np)
3572 : real(RK) , intent(out) :: Xbin(nxbin), Density(nxbin)
3573 : logical , intent(out) :: errorOccurred
3574 9 : real(RK) :: xbinsize
3575 : integer(IK) :: i, ip, thisXbin, npEffective
3576 :
3577 9 : xbinsize = (xmax-xmin) / real(nxbin,kind=RK)
3578 189 : Xbin = [ (xmin + real(i-1,kind=RK)*xbinsize,i=1,nxbin) ]
3579 :
3580 99 : Density = 0._RK
3581 9 : npEffective = 0_IK
3582 459 : do ip = 1, np
3583 459 : if (X(ip)>=xmin .and. X(ip)<xmax) then
3584 450 : npEffective = npEffective + 1_IK
3585 450 : thisXbin = getBin(X(ip),xmin,nxbin,xbinsize)
3586 450 : Density(thisXbin) = Density(thisXbin) + 1._RK
3587 : end if
3588 : end do
3589 :
3590 99 : Xbin = Xbin + 0.5_RK * xbinsize
3591 :
3592 9 : if(method=="count") then
3593 3 : errorOccurred = .false.
3594 3 : return
3595 6 : elseif(method=="pdf") then
3596 33 : Density = Density / real(npEffective,kind=RK)
3597 3 : errorOccurred = .false.
3598 : else
3599 3 : errorOccurred = .true.
3600 : end if
3601 :
3602 12 : end subroutine getHist1D
3603 :
3604 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3605 :
3606 : !> \brief
3607 : !> Return the 2-D histogram (Density plot) of a set of data points with (X,Y) coordinates.
3608 : !> The number of bins in the `X` and `Y` directions (`[nxbin, nybin]`) are determined by the user.
3609 : !> The range of `X` and `Y` (`xmin`,`xmax`,`ymin`,`ymax`) should also be given by the user.
3610 : !> The program returns three arrays of `Xbin`, `Ybin`, and `Density(y,x)` as the output.
3611 : !>
3612 : !> \param[in] histType : The method by which the normalization of the histogram counts should be done:
3613 : !> + `"count"` : Return the count histogram.
3614 : !> + `"pdf"` : Return the probability density function (PDF) histogram.
3615 : !> + `"pdf(y|x)"` : Return the conditional PDF of `y` given `x` histogram.
3616 : !> + `"pdf(x|y)"` : Return the conditional PDF of `x` given `y` histogram.
3617 : !> \param[in] xmin : The minimum of the histogram binning along the x-axis.
3618 : !> \param[in] xmax : The maximum of the histogram binning along the x-axis.
3619 : !> \param[in] ymin : The minimum of the histogram binning along the y-axis.
3620 : !> \param[in] ymax : The maximum of the histogram binning along the y-axis.
3621 : !> \param[in] nxbin : The number of histogram bins along the x-axis.
3622 : !> \param[in] nybin : The number of histogram bins along the y-axis.
3623 : !> \param[in] np : The length of input vector `X`.
3624 : !> \param[in] X : The vector of length `nxbin` of values to be binned.
3625 : !> \param[in] Y : The vector of length `nybin` of values to be binned.
3626 : !> \param[out] Xbin : The vector of length `nxbin` of values representing the bin centers.
3627 : !> \param[out] Ybin : The vector of length `nybin` of values representing the bin centers.
3628 : !> \param[out] Density : The array of shape `(nybin,nxbin)` of values representing the densities per bin.
3629 : !> \param[out] errorOccurred : A logical output value indicating whether an error has occurred.
3630 : !>
3631 : !> \author
3632 : !> Amir Shahmoradi, Sep 1, 2017, 12:00 AM, ICES, UT Austin
3633 15 : subroutine getHist2D(histType,xmin,xmax,ymin,ymax,nxbin,nybin,np,X,Y,Xbin,Ybin,Density,errorOccurred)
3634 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3635 : !DEC$ ATTRIBUTES DLLEXPORT :: getHist2D
3636 : #endif
3637 : !use, intrinsic :: iso_fortran_env, only: output_unit
3638 9 : use String_mod, only: getLowerCase
3639 : implicit none
3640 : character(*), intent(in) :: histType
3641 : integer(IK) , intent(in) :: np,nxbin,nybin
3642 : real(RK) , intent(in) :: xmin,xmax,ymin,ymax
3643 : real(RK) , intent(in) :: X(np),Y(np)
3644 : real(RK) , intent(out) :: Xbin(nxbin),Ybin(nybin),Density(nybin,nxbin)
3645 : logical , intent(out) :: errorOccurred
3646 15 : character(:), allocatable :: method
3647 15 : real(RK) :: xbinsize,ybinsize
3648 : integer(IK) :: i, ip, thisXbin, thisYbin, npEffective
3649 :
3650 15 : errorOccurred = .false.
3651 :
3652 15 : xbinsize = (xmax-xmin) / real(nxbin,kind=RK)
3653 15 : ybinsize = (ymax-ymin) / real(nybin,kind=RK)
3654 255 : Xbin = [ (xmin+real(i-1,kind=RK)*xbinsize,i=1,nxbin) ]
3655 225 : Ybin = [ (ymin+real(i-1,kind=RK)*ybinsize,i=1,nybin) ]
3656 :
3657 975 : Density = 0._RK
3658 15 : npEffective = 0_IK
3659 765 : do ip = 1,np
3660 765 : if (X(ip)>=xmin .and. X(ip)<xmax .and. Y(ip)>=ymin .and. Y(ip)<ymax) then
3661 750 : npEffective = npEffective + 1_IK
3662 750 : thisXbin = getBin(X(ip),xmin,nxbin,xbinsize)
3663 750 : thisYbin = getBin(Y(ip),ymin,nybin,ybinsize)
3664 750 : Density(thisYbin,thisXbin) = Density(thisYbin,thisXbin) + 1._RK
3665 : end if
3666 : end do
3667 :
3668 135 : Xbin = Xbin + 0.5_RK * xbinsize
3669 120 : Ybin = Ybin + 0.5_RK * ybinsize
3670 15 : method = getLowerCase(trim(adjustl(histType)))
3671 15 : if(method=="pdf") then
3672 195 : Density = Density / real(npEffective,kind=RK)
3673 12 : elseif(method=="pdf(y|x)") then
3674 27 : do i = 1,nxbin
3675 363 : Density(1:nybin,i) = Density(1:nybin,i) / sum(Density(1:nybin,i))
3676 : end do
3677 9 : elseif(method=="pdf(x|y)") then
3678 24 : do i = 1,nybin
3679 360 : Density(i,1:nxbin) = Density(i,1:nxbin) / sum(Density(i,1:nxbin))
3680 : end do
3681 6 : elseif(method=="count") then
3682 3 : return
3683 : else
3684 3 : errorOccurred = .true.
3685 : end if
3686 :
3687 15 : end subroutine getHist2D
3688 :
3689 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3690 :
3691 : !> \brief
3692 : !> Given the range of the variable `x`, `xmin:xmin+binsize*nbin`, and the number of bins, `nbin`, with which
3693 : !> this range is divided, find which bin the input value `x` falls among `[1:nbin]` bins.
3694 : !> The output `ibin` is the number that identifies the bin.
3695 : !>
3696 : !> \param[in] x : The input value whose bin ID is to be found.
3697 : !> \param[in] lowerBound : The lower limit on the value of `x`.
3698 : !> \param[in] nbin : The number of bins to be considered starting from `lowerBound`.
3699 : !> \param[in] binsize : The size of the bins. It must be exactly `(xmax - xmin) / nbin`.
3700 : !>
3701 : !> \return
3702 : !> `ibin` : The ID of the bin to which the input value `x` belongs.
3703 : !>
3704 : !> \warning
3705 : !> If `x <= xmin` or `x xmin + nbin * binsize`, then `ibin = -1` will be returned to indicate error.
3706 : !>
3707 : !> \remark
3708 : !> If `bmin < x <= bmax` then `x` belongs to this bin.
3709 : !>
3710 : !> \todo
3711 : !> The performance and interface of this routine can be significantly improved.
3712 : !> It is more sensible to pass a contiguous array of bin edges as input instead of `lowerBound` and `binsize`.
3713 : !> If the input point is not within any bins, an index of zero should be returned.
3714 : !>
3715 : !> \author
3716 : !> Version 3.0, Sep 1, 2017, 11:12 AM, Amir Shahmoradi, ICES, The University of Texas at Austin.
3717 1950 : pure function getBin(x,lowerBound,nbin,binsize) result(ibin)
3718 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3719 : !DEC$ ATTRIBUTES DLLEXPORT :: getBin
3720 : #endif
3721 :
3722 : implicit none
3723 : integer(IK), intent(in) :: nbin
3724 : real(RK) , intent(in) :: x,lowerBound,binsize
3725 1950 : real(RK) :: xmin,xmid
3726 : integer(IK) :: ibin,minbin,midbin,maxbin
3727 :
3728 1950 : if (x<lowerBound .or. x>=lowerBound+nbin*binsize) then
3729 : ! LCOV_EXCL_START
3730 : ibin = -1_IK
3731 : return
3732 : ! LCOV_EXCL_STOP
3733 : end if
3734 :
3735 1950 : minbin = 1
3736 1950 : maxbin = nbin
3737 1950 : xmin = lowerBound
3738 5547 : loopFindBin: do
3739 7497 : midbin = (minbin+maxbin) / 2
3740 7497 : xmid = xmin + midbin*binsize
3741 7497 : if (x<xmid) then
3742 3021 : if (minbin==midbin) then
3743 54 : ibin = minbin
3744 54 : exit loopFindBin
3745 : end if
3746 2967 : maxbin = midbin
3747 2967 : cycle loopFindBin
3748 : else
3749 4476 : if (minbin==midbin) then
3750 1896 : ibin = maxbin
3751 1896 : exit loopFindBin
3752 : end if
3753 2580 : minbin = midbin
3754 2580 : cycle loopFindBin
3755 : end if
3756 : end do loopFindBin
3757 :
3758 1965 : end function getBin
3759 :
3760 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3761 :
3762 : !> \brief
3763 : !> Return the quantiles of an input sample of points, given the input quantile probabilities.
3764 : !>
3765 : !> \param[in] np : The number of points in the input sample.
3766 : !> \param[in] nq : The number of output quantiles.
3767 : !> \param[in] SortedQuantileProbability : A sorted ascending-order vector of probabilities at which the quantiles will be returned.
3768 : !> \param[in] Point : The vector of length `np` representing the input sample.
3769 : !> \param[in] Weight : The vector of length `np` representing the weights of the points in the input sample.
3770 : !> \param[in] sumWeight : The sum of the vector weights of the points: `sum(Weight)`.
3771 : !>
3772 : !> \return
3773 : !> `Quantile` : The output vector of length `nq`, representing the quantiles corresponding to the input `SortedQuantileProbability` probabilities.
3774 : !>
3775 : !> \author
3776 : !> Amir Shahmoradi, Sep 1, 2017, 12:00 AM, ICES, UT Austin
3777 10969 : pure function getQuantile(np,nq,SortedQuantileProbability,Point,Weight,sumWeight) result(Quantile)
3778 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
3779 : !DEC$ ATTRIBUTES DLLEXPORT :: getQuantile
3780 : #endif
3781 1950 : use Constants_mod, only: IK, RK, NEGINF_RK
3782 : use Sort_mod, only: indexArray
3783 : use Err_mod, only: Err_type
3784 : implicit none
3785 : integer(IK) , intent(in) :: np, nq
3786 : real(RK) , intent(in) :: SortedQuantileProbability(nq), Point(np)
3787 : integer(IK) , intent(in), optional :: Weight(np), sumWeight
3788 : real(RK) :: Quantile(nq)
3789 1688 : integer(IK) :: ip, iq, iw, weightCounter, Indx(np), SortedQuantileDensity(nq)
3790 844 : type(Err_type) :: Err
3791 844 : call indexArray(np,Point,Indx,Err)
3792 844 : if (Err%occurred) then
3793 : ! LCOV_EXCL_START
3794 : Quantile = NEGINF_RK
3795 : return
3796 : ! LCOV_EXCL_STOP
3797 : end if
3798 844 : iq = 1_IK
3799 844 : if (present(sumWeight)) then
3800 841 : weightCounter = 0_IK
3801 8410 : SortedQuantileDensity = nint( SortedQuantileProbability * sumWeight )
3802 235193 : loopWeighted: do ip = 1, np
3803 954201 : do iw = 1, Weight(Indx(ip))
3804 719008 : weightCounter = weightCounter + 1_IK
3805 953360 : if (weightCounter>=SortedQuantileDensity(iq)) then
3806 7453 : Quantile(iq) = Point(Indx(ip))
3807 7453 : iq = iq + 1_IK
3808 7453 : if (iq>nq) exit loopWeighted
3809 : end if
3810 : end do
3811 : end do loopWeighted
3812 : else
3813 30 : SortedQuantileDensity = nint( SortedQuantileProbability * np )
3814 150 : loopNonWeighted: do ip = 1, np
3815 153 : if (ip>=SortedQuantileDensity(iq)) then
3816 27 : Quantile(iq) = Point(Indx(ip))
3817 27 : iq = iq + 1_IK
3818 27 : if (iq>nq) exit loopNonWeighted
3819 : end if
3820 : end do loopNonWeighted
3821 : end if
3822 960 : if (iq<=nq) Quantile(iq:nq) = Point(Indx(np))
3823 2532 : end function getQuantile
3824 :
3825 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3826 :
3827 : end module Statistics_mod ! LCOV_EXCL_LINE
|
