ParaMonte Fortran 2.0.0
Parallel Monte Carlo and Machine Learning Library
See the latest version documentation.
test_pm_hypoTest.F90
Go to the documentation of this file.
1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3!!!! !!!!
4!!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5!!!! !!!!
6!!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7!!!! !!!!
8!!!! This file is part of the ParaMonte library. !!!!
9!!!! !!!!
10!!!! LICENSE !!!!
11!!!! !!!!
12!!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13!!!! !!!!
14!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16
19
20module test_pm_statest
21
22 use pm_statest
23 use pm_test, only: test_type, LK
24 use pm_kind, only: LK
25
26 implicit none
27
28 private
29 public :: setTest
30 type(test_type) :: test
31
32 integer(IK) , parameter :: lenRnd = 50_IK
33
34 real(RK) :: UnifRnd(lenRnd) = [ 0.0759666916908419_RK &
35 , 0.2399161535536580_RK &
36 , 0.1233189348351660_RK &
37 , 0.1839077882824170_RK &
38 , 0.2399525256649030_RK &
39 , 0.4172670690843700_RK &
40 , 0.0496544303257421_RK &
41 , 0.9027161099152810_RK &
42 , 0.9447871897216460_RK &
43 , 0.4908640924680800_RK &
44 , 0.4892526384000190_RK &
45 , 0.3377194098213770_RK &
46 , 0.9000538464176620_RK &
47 , 0.3692467811202150_RK &
48 , 0.1112027552937870_RK &
49 , 0.7802520683211380_RK &
50 , 0.3897388369612530_RK &
51 , 0.2416912859138330_RK &
52 , 0.4039121455881150_RK &
53 , 0.0964545251683886_RK &
54 , 0.1319732926063350_RK &
55 , 0.9420505907754850_RK &
56 , 0.9561345402298020_RK &
57 , 0.5752085950784660_RK &
58 , 0.0597795429471558_RK &
59 , 0.2347799133724060_RK &
60 , 0.3531585712220710_RK &
61 , 0.8211940401979590_RK &
62 , 0.0154034376515551_RK &
63 , 0.0430238016578078_RK &
64 , 0.1689900294627040_RK &
65 , 0.6491154749564520_RK &
66 , 0.7317223856586700_RK &
67 , 0.6477459631363070_RK &
68 , 0.4509237064309450_RK &
69 , 0.5470088922863450_RK &
70 , 0.2963208056077730_RK &
71 , 0.7446928070741560_RK &
72 , 0.1889550150325450_RK &
73 , 0.6867754333653150_RK &
74 , 0.1835111557372700_RK &
75 , 0.3684845964903370_RK &
76 , 0.6256185607296900_RK &
77 , 0.7802274351513770_RK &
78 , 0.0811257688657853_RK &
79 , 0.9293859709687300_RK &
80 , 0.7757126786084020_RK &
81 , 0.4867916324031720_RK &
82 , 0.4358585885809190_RK &
83 , 0.4467837494298060_RK ]
84
85 real(RK) :: StdNormRnd1(lenRnd) = [ 0.537667139546100_RK &
86 , 1.83388501459509_RK &
87 , -2.25884686100365_RK &
88 , 0.862173320368121_RK &
89 , 0.318765239858981_RK &
90 , -1.30768829630527_RK &
91 , -0.433592022305684_RK &
92 , 0.342624466538650_RK &
93 , 3.57839693972576_RK &
94 , 2.76943702988488_RK &
95 , -1.34988694015652_RK &
96 , 3.03492346633185_RK &
97 , 0.725404224946106_RK &
98 , -0.0630548731896562_RK &
99 , 0.714742903826096_RK &
100 , -0.204966058299775_RK &
101 , -0.124144348216312_RK &
102 , 1.48969760778547_RK &
103 , 1.40903448980048_RK &
104 , 1.41719241342961_RK &
105 , 0.671497133608081_RK &
106 , -1.20748692268504_RK &
107 , 0.717238651328839_RK &
108 , 1.63023528916473_RK &
109 , 0.488893770311789_RK &
110 , 1.03469300991786_RK &
111 , 0.726885133383238_RK &
112 , -0.303440924786016_RK &
113 , 0.293871467096658_RK &
114 , -0.787282803758638_RK &
115 , 0.888395631757642_RK &
116 , -1.14707010696915_RK &
117 , -1.06887045816803_RK &
118 , -0.809498694424876_RK &
119 , -2.94428416199490_RK &
120 , 1.43838029281510_RK &
121 , 0.325190539456198_RK &
122 , -0.754928319169703_RK &
123 , 1.37029854009523_RK &
124 , -1.71151641885370_RK &
125 , -0.102242446085491_RK &
126 , -0.241447041607358_RK &
127 , 0.319206739165502_RK &
128 , 0.312858596637428_RK &
129 , -0.864879917324457_RK &
130 , -0.0300512961962686_RK &
131 , -0.164879019209038_RK &
132 , 0.627707287528727_RK &
133 , 1.09326566903948_RK &
134 , 1.10927329761440_RK ]
135
136 real(RK) :: StdNormRnd2(lenRnd) = [ -0.863652821988714_RK &
137 , 0.0773590911304249_RK &
138 , -1.21411704361541_RK &
139 , -1.11350074148676_RK &
140 , -0.00684932810334806_RK &
141 , 1.53263030828475_RK &
142 , -0.769665913753682_RK &
143 , 0.371378812760058_RK &
144 , -0.225584402271252_RK &
145 , 1.11735613881447_RK &
146 , -1.08906429505224_RK &
147 , 0.0325574641649735_RK &
148 , 0.552527021112224_RK &
149 , 1.10061021788087_RK &
150 , 1.54421189550395_RK &
151 , 0.0859311331754255_RK &
152 , -1.49159031063761_RK &
153 , -0.742301837259857_RK &
154 , -1.06158173331999_RK &
155 , 2.35045722400204_RK &
156 , -0.615601881466894_RK &
157 , 0.748076783703985_RK &
158 , -0.192418510588264_RK &
159 , 0.888610425420721_RK &
160 , -0.764849236567874_RK &
161 , -1.40226896933876_RK &
162 , -1.42237592509150_RK &
163 , 0.488193909859941_RK &
164 , -0.177375156618825_RK &
165 , -0.196053487807333_RK &
166 , 1.41931015064255_RK &
167 , 0.291584373984183_RK &
168 , 0.197811053464361_RK &
169 , 1.58769908997406_RK &
170 , -0.804465956349547_RK &
171 , 0.696624415849607_RK &
172 , 0.835088165072682_RK &
173 , -0.243715140377952_RK &
174 , 0.215670086403744_RK &
175 , -1.16584393148205_RK &
176 , -1.14795277889859_RK &
177 , 0.104874716016494_RK &
178 , 0.722254032225002_RK &
179 , 2.58549125261624_RK &
180 , -0.666890670701386_RK &
181 , 0.187331024578940_RK &
182 , -0.0824944253709554_RK &
183 , -1.93302291785099_RK &
184 , -0.438966153934773_RK &
185 , -1.79467884145512_RK ]
186
187 integer(IK) , parameter :: lenDistSortedDiff = 200_IK
188 real(RK) , parameter :: DistSortedDiff(*) = [ 0.00137390044027030_RK & ! LCOV_EXCL_LINE
189 , 0.00106568311889532_RK & ! LCOV_EXCL_LINE
190 , 0.00261528537896705_RK & ! LCOV_EXCL_LINE
191 , 0.00418324458049812_RK & ! LCOV_EXCL_LINE
192 , 0.00173791427946435_RK & ! LCOV_EXCL_LINE
193 , 0.00683744981995738_RK & ! LCOV_EXCL_LINE
194 , 0.00122509129524828_RK & ! LCOV_EXCL_LINE
195 , 0.000963897105374811_RK & ! LCOV_EXCL_LINE
196 , 0.00659617936766466_RK & ! LCOV_EXCL_LINE
197 , 0.00121099678988135_RK & ! LCOV_EXCL_LINE
198 , 0.00152627239164915_RK & ! LCOV_EXCL_LINE
199 , 5.59784103162375e-05_RK & ! LCOV_EXCL_LINE
200 , 0.00396276701851184_RK & ! LCOV_EXCL_LINE
201 , 0.00585926454186925_RK & ! LCOV_EXCL_LINE
202 , 0.00183541725616831_RK & ! LCOV_EXCL_LINE
203 , 0.00213464221769843_RK & ! LCOV_EXCL_LINE
204 , 0.00456380640221332_RK & ! LCOV_EXCL_LINE
205 , 0.000300222645855053_RK & ! LCOV_EXCL_LINE
206 , 0.00113346207247378_RK & ! LCOV_EXCL_LINE
207 , 0.00403412924770930_RK & ! LCOV_EXCL_LINE
208 , 0.00876326613184120_RK & ! LCOV_EXCL_LINE
209 , 0.00240591898034481_RK & ! LCOV_EXCL_LINE
210 , 0.00111791498744240_RK & ! LCOV_EXCL_LINE
211 , 0.000314496869014858_RK & ! LCOV_EXCL_LINE
212 , 0.00166278866932723_RK & ! LCOV_EXCL_LINE
213 , 0.00123356868877644_RK & ! LCOV_EXCL_LINE
214 , 0.00136174700279745_RK & ! LCOV_EXCL_LINE
215 , 0.00180033865772067_RK & ! LCOV_EXCL_LINE
216 , 0.000188696505842079_RK & ! LCOV_EXCL_LINE
217 , 0.00447901518516258_RK & ! LCOV_EXCL_LINE
218 , 0.000649229339662938_RK & ! LCOV_EXCL_LINE
219 , 0.00665965285198666_RK & ! LCOV_EXCL_LINE
220 , 0.000639392345076595_RK & ! LCOV_EXCL_LINE
221 , 0.00439231287788411_RK & ! LCOV_EXCL_LINE
222 , 0.00514580427373124_RK & ! LCOV_EXCL_LINE
223 , 0.00325615977507276_RK & ! LCOV_EXCL_LINE
224 , 0.00262444209665191_RK & ! LCOV_EXCL_LINE
225 , 0.00287172623926713_RK & ! LCOV_EXCL_LINE
226 , 0.00726524941509710_RK & ! LCOV_EXCL_LINE
227 , 0.000679801673850622_RK & ! LCOV_EXCL_LINE
228 , 0.000771848959060684_RK & ! LCOV_EXCL_LINE
229 , 0.00236687990071305_RK & ! LCOV_EXCL_LINE
230 , 0.00279738279147568_RK & ! LCOV_EXCL_LINE
231 , 0.00134507052319155_RK & ! LCOV_EXCL_LINE
232 , 0.00192351678728353_RK & ! LCOV_EXCL_LINE
233 , 0.000403579255582320_RK & ! LCOV_EXCL_LINE
234 , 0.00466075625864260_RK & ! LCOV_EXCL_LINE
235 , 0.00128562661588161_RK & ! LCOV_EXCL_LINE
236 , 0.00240934864593989_RK & ! LCOV_EXCL_LINE
237 , 0.00379834817045643_RK & ! LCOV_EXCL_LINE
238 , 0.00626106613651900_RK & ! LCOV_EXCL_LINE
239 , 0.00155470502430244_RK & ! LCOV_EXCL_LINE
240 , 0.00386451831278178_RK & ! LCOV_EXCL_LINE
241 , 0.00633652247623751_RK & ! LCOV_EXCL_LINE
242 , 0.00212716547361247_RK & ! LCOV_EXCL_LINE
243 , 0.000117519040001346_RK & ! LCOV_EXCL_LINE
244 , 0.00268124387647972_RK & ! LCOV_EXCL_LINE
245 , 0.00330190345409753_RK & ! LCOV_EXCL_LINE
246 , 0.00181556281964235_RK & ! LCOV_EXCL_LINE
247 , 0.00680782445764760_RK & ! LCOV_EXCL_LINE
248 , 0.00344385667182678_RK & ! LCOV_EXCL_LINE
249 , 0.00341876847109601_RK & ! LCOV_EXCL_LINE
250 , 0.00156480476836285_RK & ! LCOV_EXCL_LINE
251 , 0.00392916582362890_RK & ! LCOV_EXCL_LINE
252 , 0.000536610168473173_RK & ! LCOV_EXCL_LINE
253 , 0.00308379345606979_RK & ! LCOV_EXCL_LINE
254 , 0.00382371993615627_RK & ! LCOV_EXCL_LINE
255 , 0.00143530941667058_RK & ! LCOV_EXCL_LINE
256 , 0.000273761654861704_RK & ! LCOV_EXCL_LINE
257 , 0.000218939488190961_RK & ! LCOV_EXCL_LINE
258 , 0.00395298813907463_RK & ! LCOV_EXCL_LINE
259 , 0.00461683348784969_RK & ! LCOV_EXCL_LINE
260 , 0.000587765974902510_RK & ! LCOV_EXCL_LINE
261 , 0.00130109023546354_RK & ! LCOV_EXCL_LINE
262 , 0.000446086978761473_RK & ! LCOV_EXCL_LINE
263 , 0.00102849200817856_RK & ! LCOV_EXCL_LINE
264 , 0.00475349511066359_RK & ! LCOV_EXCL_LINE
265 , 0.00177004768590394_RK & ! LCOV_EXCL_LINE
266 , 0.00118757227387534_RK & ! LCOV_EXCL_LINE
267 , 0.00511260207997721_RK & ! LCOV_EXCL_LINE
268 , 0.000966468749441840_RK & ! LCOV_EXCL_LINE
269 , 0.00153506873793674_RK & ! LCOV_EXCL_LINE
270 , 0.000163947408605591_RK & ! LCOV_EXCL_LINE
271 , 0.00159865940856585_RK & ! LCOV_EXCL_LINE
272 , 0.00319324233871277_RK & ! LCOV_EXCL_LINE
273 , 0.00182318600045628_RK & ! LCOV_EXCL_LINE
274 , 0.00270039749021500_RK & ! LCOV_EXCL_LINE
275 , 0.00233393067978271_RK & ! LCOV_EXCL_LINE
276 , 0.000304246250464657_RK & ! LCOV_EXCL_LINE
277 , 0.00107643657879242_RK & ! LCOV_EXCL_LINE
278 , 0.00149167058888822_RK & ! LCOV_EXCL_LINE
279 , 3.22029433108551e-05_RK & ! LCOV_EXCL_LINE
280 , 0.000507552164931036_RK & ! LCOV_EXCL_LINE
281 , 0.00284976448963692_RK & ! LCOV_EXCL_LINE
282 , 0.00261297396275373_RK & ! LCOV_EXCL_LINE
283 , 0.00115230561975765_RK & ! LCOV_EXCL_LINE
284 , 0.000589761776301434_RK & ! LCOV_EXCL_LINE
285 , 0.00254655116972891_RK & ! LCOV_EXCL_LINE
286 , 0.000736372303162147_RK & ! LCOV_EXCL_LINE
287 , 0.00298938200986687_RK & ! LCOV_EXCL_LINE
288 , 0.00115908892632211_RK & ! LCOV_EXCL_LINE
289 , 0.0118304957984586_RK & ! LCOV_EXCL_LINE
290 , 0.00427409203243478_RK & ! LCOV_EXCL_LINE
291 , 0.00408018550920808_RK & ! LCOV_EXCL_LINE
292 , 0.00137769679358790_RK & ! LCOV_EXCL_LINE
293 , 0.00172460242999972_RK & ! LCOV_EXCL_LINE
294 , 0.000397268483832813_RK & ! LCOV_EXCL_LINE
295 , 0.00640511077756611_RK & ! LCOV_EXCL_LINE
296 , 0.00160078461814450_RK & ! LCOV_EXCL_LINE
297 , 0.00116751842007934_RK & ! LCOV_EXCL_LINE
298 , 0.00689243594329803_RK & ! LCOV_EXCL_LINE
299 , 5.96420573087952e-05_RK & ! LCOV_EXCL_LINE
300 , 0.00120353557763264_RK & ! LCOV_EXCL_LINE
301 , 0.00542771429675204_RK & ! LCOV_EXCL_LINE
302 , 0.00610071142686630_RK & ! LCOV_EXCL_LINE
303 , 0.00213356732282444_RK & ! LCOV_EXCL_LINE
304 , 0.000281572449424172_RK & ! LCOV_EXCL_LINE
305 , 0.00129330702548425_RK & ! LCOV_EXCL_LINE
306 , 0.000169205982398446_RK & ! LCOV_EXCL_LINE
307 , 0.000101693454083507_RK & ! LCOV_EXCL_LINE
308 , 0.000396590512815487_RK & ! LCOV_EXCL_LINE
309 , 0.000271783190621933_RK & ! LCOV_EXCL_LINE
310 , 0.00207176434429135_RK & ! LCOV_EXCL_LINE
311 , 0.00195525703558364_RK & ! LCOV_EXCL_LINE
312 , 0.00133407286142151_RK & ! LCOV_EXCL_LINE
313 , 0.00146450204323612_RK & ! LCOV_EXCL_LINE
314 , 0.00172854780161391_RK & ! LCOV_EXCL_LINE
315 , 0.00175309009682001_RK & ! LCOV_EXCL_LINE
316 , 0.00591372011950542_RK & ! LCOV_EXCL_LINE
317 , 0.00383246845630059_RK & ! LCOV_EXCL_LINE
318 , 0.00478224831434315_RK & ! LCOV_EXCL_LINE
319 , 0.00140155544663756_RK & ! LCOV_EXCL_LINE
320 , 0.00311278264609216_RK & ! LCOV_EXCL_LINE
321 , 0.00117361727622267_RK & ! LCOV_EXCL_LINE
322 , 0.00197417012646872_RK & ! LCOV_EXCL_LINE
323 , 0.00124210536927416_RK & ! LCOV_EXCL_LINE
324 , 0.000362752683200629_RK & ! LCOV_EXCL_LINE
325 , 0.000939106199745576_RK & ! LCOV_EXCL_LINE
326 , 0.000794249402029323_RK & ! LCOV_EXCL_LINE
327 , 0.000913827979191928_RK & ! LCOV_EXCL_LINE
328 , 0.00720437268872554_RK & ! LCOV_EXCL_LINE
329 , 0.000205371175825086_RK & ! LCOV_EXCL_LINE
330 , 0.00715773811743015_RK & ! LCOV_EXCL_LINE
331 , 0.00251452853672052_RK & ! LCOV_EXCL_LINE
332 , 0.000861445629684599_RK & ! LCOV_EXCL_LINE
333 , 0.00418790413483983_RK & ! LCOV_EXCL_LINE
334 , 0.000370461989028015_RK & ! LCOV_EXCL_LINE
335 , 0.00125679370708720_RK & ! LCOV_EXCL_LINE
336 , 0.00349303438590243_RK & ! LCOV_EXCL_LINE
337 , 0.000724203659035916_RK & ! LCOV_EXCL_LINE
338 , 0.000491862014471267_RK & ! LCOV_EXCL_LINE
339 , 0.00287196585168448_RK & ! LCOV_EXCL_LINE
340 , 0.00185611717084277_RK & ! LCOV_EXCL_LINE
341 , 0.000150069773691919_RK & ! LCOV_EXCL_LINE
342 , 0.000848534416059033_RK & ! LCOV_EXCL_LINE
343 , 0.00123442270407526_RK & ! LCOV_EXCL_LINE
344 , 0.00103148223612315_RK & ! LCOV_EXCL_LINE
345 , 0.00132400045208614_RK & ! LCOV_EXCL_LINE
346 , 0.00435456041987981_RK & ! LCOV_EXCL_LINE
347 , 0.00117222896204028_RK & ! LCOV_EXCL_LINE
348 , 0.000873409715677065_RK & ! LCOV_EXCL_LINE
349 , 0.00512459761807316_RK & ! LCOV_EXCL_LINE
350 , 0.000423255105326148_RK & ! LCOV_EXCL_LINE
351 , 0.00801979475421000_RK & ! LCOV_EXCL_LINE
352 , 0.00452332834787750_RK & ! LCOV_EXCL_LINE
353 , 0.000662717679231761_RK & ! LCOV_EXCL_LINE
354 , 0.00612674945114589_RK & ! LCOV_EXCL_LINE
355 , 0.000241023672676199_RK & ! LCOV_EXCL_LINE
356 , 0.00241368458759861_RK & ! LCOV_EXCL_LINE
357 , 0.00721401234545482_RK & ! LCOV_EXCL_LINE
358 , 0.000617822886748609_RK & ! LCOV_EXCL_LINE
359 , 4.20361472183162e-05_RK & ! LCOV_EXCL_LINE
360 , 0.000176033863938718_RK & ! LCOV_EXCL_LINE
361 , 0.00261140736592547_RK & ! LCOV_EXCL_LINE
362 , 0.00161441323061740_RK & ! LCOV_EXCL_LINE
363 , 0.00119563944997603_RK & ! LCOV_EXCL_LINE
364 , 0.000161186670205815_RK & ! LCOV_EXCL_LINE
365 , 0.00186893825670897_RK & ! LCOV_EXCL_LINE
366 , 0.0156568954429827_RK & ! LCOV_EXCL_LINE
367 , 0.00333188871092183_RK & ! LCOV_EXCL_LINE
368 , 0.00460140509714901_RK & ! LCOV_EXCL_LINE
369 , 0.00106350049875970_RK & ! LCOV_EXCL_LINE
370 , 0.000278484374323207_RK & ! LCOV_EXCL_LINE
371 , 0.000109843896551887_RK & ! LCOV_EXCL_LINE
372 , 0.00517794656024040_RK & ! LCOV_EXCL_LINE
373 , 0.000764568679501587_RK & ! LCOV_EXCL_LINE
374 , 0.00106947825226744_RK & ! LCOV_EXCL_LINE
375 , 0.00475394751147906_RK & ! LCOV_EXCL_LINE
376 , 0.00700310227105183_RK & ! LCOV_EXCL_LINE
377 , 0.00280107346803815_RK & ! LCOV_EXCL_LINE
378 , 0.00594747756623304_RK & ! LCOV_EXCL_LINE
379 , 0.00268586958531547_RK & ! LCOV_EXCL_LINE
380 , 0.00120417364527792_RK & ! LCOV_EXCL_LINE
381 , 0.00344542024642325_RK & ! LCOV_EXCL_LINE
382 , 0.00217032516670990_RK & ! LCOV_EXCL_LINE
383 , 0.00170145556378476_RK & ! LCOV_EXCL_LINE
384 , 0.00159968663729304_RK & ! LCOV_EXCL_LINE
385 , 0.00153233985478085_RK & ! LCOV_EXCL_LINE
386 , 0.00219750724001799_RK & ! LCOV_EXCL_LINE
387 , 0.00230064745992042_RK & ! LCOV_EXCL_LINE
388 ]
389 integer(IK) , parameter :: SampleSize(*) = [ 1_IK &
390 , 101_IK &
391 , 102_IK &
392 , 103_IK &
393 , 104_IK &
394 , 105_IK &
395 , 106_IK &
396 , 107_IK &
397 , 108_IK &
398 , 109_IK &
399 , 110_IK &
400 ]
401 integer(IK) , parameter :: lenSampleSize = size(SampleSize,kind=IK)
402 real(RK) , parameter :: ProbKS_REF(*) = [ 1._RK &
403 , 0.482474176998994_RK &
404 , 0.370817196700337_RK &
405 , 0.493146633387423_RK &
406 , 0.702120402082310_RK &
407 , 0.671010642014021_RK &
408 , 0.713612126474551_RK &
409 , 0.683487750289677_RK &
410 , 0.559321682773637_RK &
411 , 0.530113153332594_RK &
412 , 0.571577075807556_RK &
413 ] ! The KS two-sample test probabilities.
414 real(RK) , parameter :: StatKS_REF(*) = [ 0._RK &
415 , 0.162352941176471_RK &
416 , 0.176470588235294_RK &
417 , 0.159502262443439_RK &
418 , 0.134615384615385_RK &
419 , 0.137518142235123_RK &
420 , 0.132075471698113_RK &
421 , 0.134870719776380_RK &
422 , 0.148148148148148_RK &
423 , 0.150841750841751_RK &
424 , 0.145454545454545_RK &
425 ] ! The KS two-sample test statistics.
426
427!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
428
429contains
430
431!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
432
433 subroutine setTest()
434
435 test = test_type(MODULE_NAME)
436 call test_Uniformity()
437 return
438 call test%run(test_doKS1_1, SK_"test_doKS1_1")
439 call test%run(test_doSortedKS2_1, SK_"test_doSortedKS2_1")
440 call test%run(test_doUniformKS1_1, SK_"test_doUniformKS1_1")
441 call test%run(test_getCumDenComKS_1, SK_"test_getCumDenComKS_1")
442 call test%run(test_getSampleDisparity_1, SK_"test_getSampleDisparity_1")
443 call test%run(test_getSampleDisparity_2, SK_"test_getSampleDisparity_2")
444 !call test%run(test_performance_1, SK_"test_performance_1")
445 call test%summarize()
446
447 end subroutine setTest
448
449!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
450
451 function test_doKS1_1() result(assertion)
452
453 use pm_distNorm, only: getNormCDF
454 use pm_kind, only: RK, IK
455
456 implicit none
457
458 logical(LK) :: assertion
459 real(RK) , parameter :: TOLERANCE = 1.e-7_RK
460
461 real(RK), parameter :: probKS_ref = .3763758622852317E-01_RK
462 real(RK), parameter :: statKS_ref = .1955719390701096_RK
463 real(RK) :: statKS
464 real(RK) :: probKS
465 real(RK) :: difference
466
467 call doKS1( np = lenRnd &
468 , Point = StdNormRnd1 &
469 , getCDF = getCDF &
470 , statKS = statKS &
471 , probKS = probKS &
472 )
473
474 difference = abs( (probKS - probKS_ref) / probKS_ref )
475 assertion = difference < TOLERANCE
476
477 if (test%traceable .and. .not. assertion) then
478 ! LCOV_EXCL_START
479 write(test%disp%unit,"(*(g0,:,', '))")
480 write(test%disp%unit,"(*(g0,:,', '))") "probKS_ref :", probKS_ref
481 write(test%disp%unit,"(*(g0,:,', '))") "probKS :", probKS
482 write(test%disp%unit,"(*(g0,:,', '))") "difference :", difference
483 write(test%disp%unit,"(*(g0,:,', '))")
484 ! LCOV_EXCL_STOP
485 end if
486
487 difference = abs( (statKS - statKS_ref) / statKS_ref )
488 assertion = assertion .and. difference < TOLERANCE
489 call test%assert(assertion)
490
491 if (test%traceable .and. .not. assertion) then
492 ! LCOV_EXCL_START
493 write(test%disp%unit,"(*(g0,:,', '))")
494 write(test%disp%unit,"(*(g0,:,', '))") "statKS_ref :", statKS_ref
495 write(test%disp%unit,"(*(g0,:,', '))") "statKS :", statKS
496 write(test%disp%unit,"(*(g0,:,', '))") "difference :", difference
497 write(test%disp%unit,"(*(g0,:,', '))")
498 ! LCOV_EXCL_STOP
499 end if
500
501 contains
502
503 PURE function getCDF(x) result(cdf)
504 implicit none
505 real(RK) , intent(in) :: x
506 real(RK) :: cdf
507 cdf = getNormCDF(x)
508 end function
509
510 end function test_doKS1_1
511
512!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
513
514 function test_getCumDenComKS_1() result(assertion)
515 use pm_kind, only: RK, IK
516 implicit none
517 integer(IK) :: i
518 logical(LK) :: assertion
519 real(RK) , parameter :: TOLERANCE = 1.e-11_RK
520 real(RK) , parameter :: Score(*) = [(0.01_RK * real(i,RK), i = -1, 100)]
521 real(RK) :: Difference(size(Score))
522 real(RK) :: ProbKS_ref(size(Score))
523 real(RK) :: ProbKS(size(Score))
524 do i = 1, size(Score)
525 ProbKS_ref(i) = getProbKS(Score(i))
526 ProbKS(i) = getCumDenComKS(Score(i))
527 Difference(i) = abs(ProbKS(i) - ProbKS_ref(i))
528 end do
529 assertion = all(Difference < TOLERANCE)
530 if (test%traceable .and. .not. assertion) then
531 ! LCOV_EXCL_START
532 write(test%disp%unit,"(*(g0,:,', '))")
533 write(test%disp%unit,"(*(g0,:,', '))") "ProbKS_ref :", ProbKS_ref
534 write(test%disp%unit,"(*(g0,:,', '))") "ProbKS :", ProbKS
535 write(test%disp%unit,"(*(g0,:,', '))") "Difference :", Difference
536 write(test%disp%unit,"(*(g0,:,', '))") "maxval(Difference) :", maxval(Difference)
537 write(test%disp%unit,"(*(g0,:,', '))")
538 ! LCOV_EXCL_STOP
539 end if
540 end function test_getCumDenComKS_1
541
542!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
543
547 function test_performance_1() result(assertion)
548 use pm_kind, only: RK, IK
549 use pm_timer, only: timerSYS_type
550 use pm_err, only: err_type
551 implicit none
552 integer(IK) :: i
553 integer(IK) , parameter :: NSCORE = 10**6_IK
554 logical(LK) :: assertion
555 real(RK) :: Score(NSCORE)
556 real(RK) :: ProbKS(NSCORE), dummy
557 type(timerSYS_type) :: timer
558 timer = timerSYS_type()
559 assertion = .true.
560 do i = -1, NSCORE - 2
561 Score(i+2) = 0.01_RK * real(i,RK)
562 end do
563 timer%start = timer%time()
564 do i = 1, NSCORE
565 ProbKS(i) = getCumDenComKS(Score(i))
566 end do
567 timer%delta = timer%time(since = timer%start)
568 write(test%disp%unit,"(*(g0,:,', '))")
569 write(test%disp%unit,"(*(g0,:,', '))") "time(getCumDenComKS):", timer%delta
570 dummy = timer%delta
571 timer%start = timer%time()
572 do i = 1, NSCORE
573 ProbKS(i) = getProbKS(Score(i))
574 end do
575 timer%delta = timer%time(since = timer%start)
576 write(test%disp%unit,"(*(g0,:,', '))") "time(getProbKS):", timer%delta
577 write(test%disp%unit,"(*(g0,:,', '))") "ratio(getProbKS/getCumDenComKS):", timer%delta / dummy
578 write(test%disp%unit,"(*(g0,:,', '))")
579 end function test_performance_1
580
581!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
582
583 function test_doUniformKS1_1() result(assertion)
584
585 use pm_kind, only: RK, IK
586
587 implicit none
588
589 logical(LK) :: assertion
590 real(RK) , parameter :: TOLERANCE = 1.e-7_RK
591
592 real(RK), parameter :: probKS_ref = .1982797523608350_RK
593 real(RK), parameter :: statKS_ref = .1491359075319200_RK
594 real(RK) :: statKS
595 real(RK) :: probKS
596 real(RK) :: difference
597
598 call doUniformKS1( np = lenRnd &
599 , Point = UnifRnd &
600 , statKS = statKS &
601 , probKS = probKS &
602 )
603
604 difference = abs( (probKS - probKS_ref) / probKS_ref )
605 assertion = difference < TOLERANCE
606 call test%assert(assertion)
607
608 if (test%traceable .and. .not. assertion) then
609 ! LCOV_EXCL_START
610 write(test%disp%unit,"(*(g0,:,', '))")
611 write(test%disp%unit,"(*(g0,:,', '))") "probKS_ref :", probKS_ref
612 write(test%disp%unit,"(*(g0,:,', '))") "probKS :", probKS
613 write(test%disp%unit,"(*(g0,:,', '))") "difference :", difference
614 write(test%disp%unit,"(*(g0,:,', '))")
615 ! LCOV_EXCL_STOP
616 end if
617
618 difference = abs( (statKS - statKS_ref) / statKS_ref )
619 assertion = assertion .and. difference < TOLERANCE
620 call test%assert(assertion)
621
622 if (test%traceable .and. .not. assertion) then
623 ! LCOV_EXCL_START
624 write(test%disp%unit,"(*(g0,:,', '))")
625 write(test%disp%unit,"(*(g0,:,', '))") "statKS_ref :", statKS_ref
626 write(test%disp%unit,"(*(g0,:,', '))") "statKS :", statKS
627 write(test%disp%unit,"(*(g0,:,', '))") "difference :", difference
628 write(test%disp%unit,"(*(g0,:,', '))")
629 ! LCOV_EXCL_STOP
630 end if
631
632 end function test_doUniformKS1_1
633
634!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
635
636 function test_doSortedKS2_1() result(assertion)
637
638 use pm_kind, only: RK, IK
639 use pm_arraySort, only: setSorted
640
641 implicit none
642
643 logical(LK) :: assertion
644 integer(IK) , parameter :: lenRnd = 50_IK
645 real(RK) , parameter :: TOLERANCE = 1.e-12_RK
646
647 real(RK), parameter :: probKS_ref = 0.056045859714425_RK
648 real(RK), parameter :: statKS_ref = 0.260000000000000_RK
649 real(RK) :: difference, probKS, statKS, lambda
650
651 assertion = .true._LK
652
653 call setSorted(StdNormRnd1)
654 call setSorted(StdNormRnd2)
655
656 call doSortedKS2( sortedSample1 = StdNormRnd1 &
657 , SortedSample2 = StdNormRnd2 &
658 , probKS = probKS &
659 , statKS = statKS &
660 , lambda = lambda &
661 )
662
663 difference = 2 * abs(probKS - probKS_ref) / (probKS_ref + probKS)
664 assertion = assertion .and. difference < TOLERANCE
665
666 !write(*,*) "statKS, diff", statKS, abs(statKS - statKS_ref)
667 !write(*,*) "probKS, diff", probKS, abs(probKS - probKS_ref)
668 if (test%traceable .and. .not. assertion) then
669 ! LCOV_EXCL_START
670 write(test%disp%unit,"(*(g0,:,', '))")
671 write(test%disp%unit,"(*(g0,:,', '))") "probKS_ref :", probKS_ref
672 write(test%disp%unit,"(*(g0,:,', '))") "probKS :", probKS
673 write(test%disp%unit,"(*(g0,:,', '))") "difference :", difference
674 write(test%disp%unit,"(*(g0,:,', '))") "TOLERANCE :", TOLERANCE
675 write(test%disp%unit,"(*(g0,:,', '))")
676 ! LCOV_EXCL_STOP
677 end if
678
679 call test%assert(assertion)
680
681 difference = 2 * abs(statKS - statKS_ref) / (statKS_ref + statKS)
682 assertion = assertion .and. difference < TOLERANCE
683 call test%assert(assertion)
684
685 if (test%traceable .and. .not. assertion) then
686 ! LCOV_EXCL_START
687 write(test%disp%unit,"(*(g0,:,', '))")
688 write(test%disp%unit,"(*(g0,:,', '))") "statKS_ref :", statKS_ref
689 write(test%disp%unit,"(*(g0,:,', '))") "statKS :", statKS
690 write(test%disp%unit,"(*(g0,:,', '))") "difference :", difference
691 write(test%disp%unit,"(*(g0,:,', '))") "TOLERANCE :", TOLERANCE
692 write(test%disp%unit,"(*(g0,:,', '))")
693 ! LCOV_EXCL_STOP
694 end if
695
696 end function test_doSortedKS2_1
697
698!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
699
702 function test_getSampleDisparity_1() result(assertion)
703
704 use pm_kind, only: RK, IK
705
706 implicit none
707
708 logical(LK) :: assertion
709 real(RK) :: ProbKS(2)
710 real(RK) :: StatKS(2)
711 real(RK) :: Lambda(2)
712
713 assertion = .true._LK
714 call setSampleDisparityAuto( SortedSample = DistSortedDiff & ! LCOV_EXCL_LINE
715 , SampleSize = [3, 4] & ! [0_IK, 1_IK] & ! LCOV_EXCL_LINE
716 , StatKS = StatKS & ! LCOV_EXCL_LINE
717 , ProbKS = ProbKS & ! LCOV_EXCL_LINE
718 , Lambda = Lambda & ! LCOV_EXCL_LINE
719 )
720 call test%assert(assertion)
721
722 call setSampleDisparityAuto( SortedSample = DistSortedDiff & ! LCOV_EXCL_LINE
723 , SampleSize = [ 1, lenDistSortedDiff] & ! [1_IK, lenDistSortedDiff + 1_IK] & ! LCOV_EXCL_LINE
724 , StatKS = StatKS & ! LCOV_EXCL_LINE
725 , ProbKS = ProbKS & ! LCOV_EXCL_LINE
726 , Lambda = Lambda & ! LCOV_EXCL_LINE
727 )
728 call test%assert(assertion)
729
730 end function test_getSampleDisparity_1
731
732!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
733
736 function test_getSampleDisparity_2() result(assertion)
737
738 use pm_kind, only: RK, IK
739
740 implicit none
741
742 logical(LK) :: assertion
743 real(RK) , allocatable :: Difference(:)
744 real(RK) , parameter :: TOLERANCE = 1.e-11_RK
745 real(RK) :: ProbKS(lenSampleSize)
746 real(RK) :: StatKS(lenSampleSize)
747 real(RK) :: Lambda(lenSampleSize)
748
749 assertion = .true.
750 call setSampleDisparityAuto( SortedSample = DistSortedDiff & ! LCOV_EXCL_LINE
751 , SampleSize = SampleSize & ! LCOV_EXCL_LINE
752 , ProbKS = ProbKS & ! LCOV_EXCL_LINE
753 , StatKS = StatKS & ! LCOV_EXCL_LINE
754 , Lambda = Lambda & ! LCOV_EXCL_LINE
755 )
756 call test%assert(assertion)
757 Difference = 2 * abs(ProbKS - ProbKS_REF) / (ProbKS_REF + ProbKS)
758 assertion = assertion .and. all(Difference < TOLERANCE)
759 call test%assert(assertion)
760
761 !write(*,*) "ProbKS, diff", ProbKS, abs(ProbKS - ProbKS_REF)
762 !write(*,*) "StatKS, diff", StatKS, abs(StatKS - StatKS_REF)
763 !write(*,*) "Lambda, diff", Lambda, abs(Lambda - StatKS_REF)
764
765 if (test%traceable .and. .not. assertion) then
766 ! LCOV_EXCL_START
767 write(test%disp%unit,"(*(g0,:,', '))")
768 write(test%disp%unit,"(*(g0,:,', '))") "ProbKS_REF:", ProbKS_REF
769 write(test%disp%unit,"(*(g0,:,', '))") "ProbKS :", ProbKS
770 write(test%disp%unit,"(*(g0,:,', '))") "Difference :", Difference
771 write(test%disp%unit,"(*(g0,:,', '))")
772 end if
773 ! LCOV_EXCL_STOP
774
775 Difference = abs(StatKS - StatKS_REF) !* 2 / (StatKS_REF + StatKS)
776 assertion = assertion .and. all(Difference < TOLERANCE)
777 call test%assert(assertion)
778
779 if (test%traceable .and. .not. assertion) then
780 ! LCOV_EXCL_START
781 write(test%disp%unit,"(*(g0,:,', '))")
782 write(test%disp%unit,"(*(g0,:,', '))") "StatKS_REF :", StatKS_REF
783 write(test%disp%unit,"(*(g0,:,', '))") "StatKS :", StatKS
784 write(test%disp%unit,"(*(g0,:,', '))") "Difference :", Difference
785 write(test%disp%unit,"(*(g0,:,', '))")
786 end if
787 ! LCOV_EXCL_STOP
788
789 end function test_getSampleDisparity_2
790
791!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
792
793 subroutine test_Uniformity()
794
795 use pm_kind, only: IK, LK, RK, RKS
796 use pm_sysPath, only: getDirCurrent
797 use pm_val2str, only: getStr
798 use pm_matrixInit, only: getMatInit
799 use pm_matrixInit, only: uppLowDia
800 use pm_matrixDet, only: setMatDetSqrtLog, uppDia, lowDia, transHerm
801 use pm_distUnifEll, only: setUnifEllRand
802 use pm_statest, only: setSampleDisparityAuto
803 !use pm_statest, only: doSortedKS2
804 use pm_arrayRange, only: getRange
805 use pm_mathCumSum, only: setCumSum
806 use pm_distanceEuclid, only: getDisEuclid, euclidsq
807 use pm_knn, only: setKnnSorted
808 use pm_knn, only: setDisSortedExpDiff
809 use pm_distanceEuclid, only: setDisMatEuclid
810 use pm_distanceEuclid, only: rdpack, uppLowDia, euclidu
811 use pm_ellipsoid, only: getLogVolUnitBall
812 use pm_timer, only: timerSYS_type
813
814 type(timerSYS_type) :: timer
815 integer(IK) :: ipnt, idim, iref
816 integer(IK) :: npntMinusOne
817 integer(IK) :: npntHalf
818 integer(IK) :: npnt
819 integer(IK) :: ndim
820 real(RK) :: minDistanceSq
821 real(RK) :: sqrtDetInvCovMat
822 real(RK) :: rho
823 real(RK) , allocatable :: refDistMat(:,:)
824 !real(RK) , allocatable :: PairDistMat(:,:)
825 !real(RK) , allocatable :: PairDistVec(:)
826 real(RK) , allocatable :: Reference(:)
827 real(RK) , allocatable :: Center(:)
828 real(RK) , allocatable :: Volume(:)
829 real(RK) , allocatable :: VolumeWeighted(:)
830 real(RK) , allocatable :: disSortedExpDiff(:)
831 real(RK) , allocatable :: sample(:,:), SampleDum(:,:) ! ndim, npnt
832 real(RK) , allocatable :: chocov(:,:)
833 real(RK) , allocatable :: DistanceSq(:)
834 real(RK) , allocatable :: ProbKS(:), StatKS(:), LambdaKS(:)
835 real(RK) , allocatable :: CumSumProbKS(:)
836 integer(IK) , allocatable :: ProxyIndex(:,:)
837 integer(IK) , allocatable :: IndexNN(:)
838 integer(IK) , allocatable :: IndexKS(:)
839 integer(IK) :: fileUnit
840 integer(IK) :: counter
841 integer(IK) :: info
842 character(*, SK), parameter :: CWD = SK_"/mnt/d/Dropbox/Projects/20220601_nsfCareer/codes/logz"
843 character(:, SK), allocatable :: filePath
844
845 rho = 0.9_RK
846 ndim = 2_IK
847 npnt = 2000_IK
848 Center = [(0._RK, idim = 1, ndim)]
849
850 ! Get the Covariance.
851
852 chocov = getMatInit([ndim, ndim + 1_IK], uppLowDia, rho, rho, 1._RK, doff = 1_IK)
853
854 ! Get the Cholesky lower factor.
855
856 call setMatDetSqrtLog(chocov(:, 2 : ndim + 1), uppDia, sqrtDetInvCovMat, info, chocov(:, 1 : ndim), transHerm)
857 if (info /= 0_IK) error stop "setMatChol() failed."
858
859 ! Normalize the volume to 1.
860
861 sqrtDetInvCovMat = exp(-sqrtDetInvCovMat / ndim)
862 do idim = 1, ndim
863 chocov(idim + 1 : ndim, idim) = chocov(idim + 1 : ndim, idim) * sqrtDetInvCovMat! / exp(getLogVolUnitBall(ndim)/ndim)
864 end do
865
866 ! Get the Cholesky lower factor.
867
868 allocate(sample(ndim, npnt))
869 call setUnifEllRand(sample, chocov, lowDia)
870
871 ! Find the closest point to the center.
872
873 minDistanceSq = +huge(minDistanceSq)
874 allocate(DistanceSq(npnt))
875 do ipnt = 1_IK, npnt
876 DistanceSq(ipnt) = getDisEuclid(center, sample(1:ndim, ipnt), euclidsq)
877 if (minDistanceSq > DistanceSq(ipnt)) then
878 minDistanceSq = DistanceSq(ipnt)
879 iref = ipnt
880 end if
881 end do
882 Reference = sample(1:ndim, iref)
883
884 ! Create a donut. \warning this may not be consistent with the above center point.
885
886 if (.false.) then
887 allocate(SampleDum, mold = sample)
888 counter = 0_IK
889 !DistanceSq = sum(sample**2, dim = 1)
890 !call setDisEuclid(distancesq, sample, spread(center, 2_ik, npnt))
891 do ipnt = 1_IK, npnt
892 !DistanceSq(ipnt) = getDisEuclid([(0._RK, idim = 1, ndim)], sample(:, ipnt), euclidsq)
893 if (DistanceSq(ipnt) < 0.4_RK**2 .or. DistanceSq(ipnt) > 0.6_RK**2) then
894 counter = counter + 1_IK
895 SampleDum(1:ndim, counter) = sample(1:ndim, ipnt)
896 end if
897 end do
898 npnt = counter
899 sample = SampleDum(1:ndim,1:counter)
900 deallocate(SampleDum)
901 end if
902
903 ! Define the sample sizes for which the KS test must be run.
904
905 npntMinusOne = npnt - 1
906 npntHalf = npntMinusOne / 2
907 allocate(disSortedExpDiff(0:npnt-1), IndexNN(npnt))
908
909 ! Write the points to the output file
910
911 filePath = CWD//SK_"/sample.txt"
912 open(newunit = fileUnit, file = filePath, status = "replace")
913 !print *, filePath
914 do ipnt = 1, npnt
915 !write(fileUnit, "("//getStr(ndim)//"(g0,:,','))") sample
916 write(fileUnit, "(*(g0,:,','))") sample(:, ipnt)
917 end do
918 close(fileUnit)
919
920 ! Sort the sample from the closest to the farthest with respect to the best reference point.
921
922 !print *, maxval(Reference)
923 !Center = [(0._RK, idim = 1, ndim)]
924 call setDisSortedExpDiff(sample, Center, disSortedExpDiff(0:npnt-1), IndexNN)
925 sample(1:ndim, 1:npnt) = sample(1:ndim, IndexNN)
926
927 ! Compute the pairwise distance matrix.
928
929 !!allocate(PairDistVec((npnt**2 - npnt) / 2))
930 !!call setPairDistSqVec(PairDistVec, sample)
931 !!PairDistVec = sqrt(PairDistVec)
932 !allocate(PairDistMat(npnt, npnt))
933 !call setPairDistSqMat(PairDistMat, sample)
934 !PairDistMat = sqrt(PairDistMat)
935
936 ! Compute the pairwise distance differences matrix.
937
938 allocate(refDistMat(npnt, npnt), ProxyIndex(npnt, npnt))
939 call setDisMatEuclid(refDistMat, rdpack, uppLowDia, sample, euclidu)
940 call setKnnSorted(refDistMat, ProxyIndex)
941
942 ! Estimate the volume.
943
944 call doCrossKS()
945 call doAutoKS()
946
947 contains
948
949 subroutine doCrossKS()
950
951 use pm_statest, only: doKS2
952 integer(IK) :: counter
953 integer(IK) :: lenIndexKS
954 integer(IK) :: maxIndexKS
955 integer(IK) :: oldIndexKS
956 integer(IK) :: index, knn
957 integer(IK) :: skip, jpnt, kpnt
958 integer(IK) :: jpntIndexInRefDistMat
959 real(RK) :: VolSortedDiff(npnt, npnt)
960 real(RK) :: VolSortedDiffAlt((npntMinusOne**2 - npntMinusOne))
961 real(RK) :: lambdaKS, statKS
962 real(RK) :: margin
963 real(RK) , allocatable :: ProbKS(:), Volume(:), VolumeWeighted(:), CumSumProbKS(:)
964
965 ! Compute the distance difference matrix.
966
967 do ipnt = 1_IK, npnt
968 VolSortedDiff(1:npnt, ipnt) = refDistMat(ProxyIndex(1:npnt, ipnt), ipnt)**ndim
969 do jpnt = npnt, 2_IK, -1_IK
970 VolSortedDiff(jpnt, ipnt) = VolSortedDiff(jpnt, ipnt) - VolSortedDiff(jpnt - 1, ipnt)
971 end do
972 end do
973 !VolSortedDiff(1:npnt, 1) = disSortedExpDiff
974
975 skip = 1_IK
976 IndexKS = [(idim, idim = 8_IK, npnt, skip)]
977 maxIndexKS = maxval(IndexKS, dim = 1)
978 lenIndexKS = size(IndexKS,1,IK)
979 allocate( ProbKS(maxIndexKS) &
980 , Volume(maxIndexKS) &
981 , VolumeWeighted(maxIndexKS) &
982 )
983
984 ! Get the KS tests
985
986 timer%start = timer%time()
987
988 ! First find all bounded neighbors within the specified index.
989
990 loopOverReference: do ipnt = 1_IK, 1!npnt
991
992 oldIndexKS = 1_IK
993 do index = 1_IK, lenIndexKS
994 knn = IndexKS(index)
995 !print *, ipnt, knn
996 counter = 0_IK
997 do jpnt = 1_IK, knn - 1_IK
998 ! find the farthest neighbor for subcircle jpnt within reference and knn.
999 ! The last point on the border has no neighbor within the super-circle, so do not loop over it.
1000 ! The jpnt'th neighbor with respect to the ipnt'th reference point has
1001 ! distance refDistMat(ProxyIndex(jpnt, ipnt), ipnt) from the reference.
1002 ! ProxyIndex(jpnt, ipnt) is the column of jpnt'th neighbor of ipnt reference in refDistMat.
1003 jpntIndexInRefDistMat = ProxyIndex(jpnt, ipnt) ! Note that `refDistMat` is symmetric.
1004 margin = refDistMat(ProxyIndex(knn, ipnt), ipnt) - refDistMat(jpntIndexInRefDistMat, ipnt) ! distance to border.
1005 !print *, margin
1006 !if (ProxyIndex(knn, ipnt) < jpntIndexInRefDistMat) error stop "nothing else matters."
1007 if (margin <= 0._RK) error stop "negative margin = "//getStr(margin)
1008 ! Loop over the neighbors of jpnt'th reference within the neighborhood of ipnt.
1009 ! the first point is the reference jpnt, so ignore it.
1010 ! But we do not know which point is the last neighbor of jpnt that is closer than `margin` to the jpnt'th reference.
1011 ! So loop over all points in the column. This is effectively a linear search,
1012 ! but can be done much faster using binary search.
1013 do kpnt = 2_IK, npnt
1014 if (refDistMat(ProxyIndex(kpnt, jpntIndexInRefDistMat), jpntIndexInRefDistMat) > margin) exit
1015 counter = counter + 1_IK
1016 VolSortedDiffAlt(counter) = VolSortedDiff(kpnt, jpnt)
1017 end do
1018 !print *, counter, refDistMat(ProxyIndex(knn, ipnt), ipnt), refDistMat(ProxyIndex(jpnt, ipnt), ipnt), margin
1019 end do
1020 if (counter == 0_IK) error stop "counter == 0_IK"
1021 call doKS2(VolSortedDiff(1:knn, ipnt), VolSortedDiffAlt(1:counter), ProbKS(knn), statKS, lambdaKS)
1022 ProbKS(oldIndexKS : knn - 1) = ProbKS(knn)
1023 oldIndexKS = knn + 1_IK
1024 end do
1025 write(*,*) "time = ", timer%time(since = timer%start)
1026
1027 ! Write the KS probabilities to the output file
1028
1029 filePath = CWD//SK_"/CrossProbKS.txt"
1030 open(newunit = fileUnit, file = filePath, status = "replace")
1031 do index = 1, maxIndexKS
1032 write(fileUnit, "(*(g0,:,','))") index, ProbKS(index)
1033 end do
1034 close(fileUnit)
1035
1036 ! Compute the mean volume
1037
1038 call setCumSum(Volume, VolSortedDiff(1:maxIndexKS,1))
1039 Volume = npnt * Volume / getRange(1_IK, maxIndexKS)
1040
1041 allocate(CumSumProbKS, mold = ProbKS)
1042 call setCumSum(CumSumProbKS, ProbKS)
1043 call setCumSum(VolumeWeighted, VolSortedDiff(1:maxIndexKS,1) * ProbKS)
1044 VolumeWeighted = npnt * VolumeWeighted / CumSumProbKS
1045
1046 filePath = CWD//SK_"/CrossVolumeKS.txt"
1047 open(newunit = fileUnit, file = filePath, status = "replace")
1048 do index = 1, maxIndexKS
1049 write(fileUnit, "(*(g0,:,','))") Volume(index), VolumeWeighted(index)
1050 end do
1051 close(fileUnit)
1052
1053 end do loopOverReference
1054
1055 end subroutine
1056
1057 subroutine doAutoKS()
1058
1059 allocate(Volume(npntHalf), VolumeWeighted(npntHalf), ProbKS(npntHalf), StatKS(npntHalf), LambdaKS(npntHalf))
1060 IndexKS = [(idim, idim = 1_IK, npntHalf, 1_IK)]
1061
1062 ! Get the KS tests
1063
1064 timer%start = timer%time()
1065 iref = 1_IK
1066 disSortedExpDiff(0:npntMinusOne) = refDistMat(ProxyIndex(:,iref), iref)**ndim
1067 disSortedExpDiff(1:npntMinusOne) = disSortedExpDiff(1:npntMinusOne) - disSortedExpDiff(0:npntMinusOne-1)
1068 call setSampleDisparityAuto(disSortedExpDiff(1:npntMinusOne), IndexKS, ProbKS, StatKS, LambdaKS)
1069 write(*,*) "time = ", timer%time(since = timer%start)
1070
1071 ! Write the KS probabilities to the output file
1072
1073 filePath = CWD//SK_"/AutoProbKS.txt"
1074 open(newunit = fileUnit, file = filePath, status = "replace")
1075 do ipnt = 1, size(IndexKS, 1, IK)
1076 write(fileUnit, "(*(g0,:,','))") 2*IndexKS(ipnt), ProbKS(ipnt)
1077 end do
1078 close(fileUnit)
1079
1080 ! Compute the mean volume
1081
1082 call setCumSum(Volume, disSortedExpDiff(1:npntHalf))
1083 Volume = npntMinusOne * Volume / getRange(1_IK, size(Volume, kind = IK))
1084
1085 allocate(CumSumProbKS, mold = ProbKS)
1086 call setCumSum(CumSumProbKS, ProbKS)
1087 call setCumSum(VolumeWeighted, disSortedExpDiff(1:npntHalf) * ProbKS)
1088 VolumeWeighted = npntMinusOne * VolumeWeighted / CumSumProbKS
1089
1090 filePath = CWD//SK_"/AutoVolumeKS.txt"
1091 open(newunit = fileUnit, file = filePath, status = "replace")
1092 do ipnt = 1, size(Volume, kind = IK)
1093 write(fileUnit, "(*(g0,:,','))") Volume(ipnt), VolumeWeighted(ipnt)
1094 end do
1095 close(fileUnit)
1096
1097 end subroutine
1098
1099 end subroutine
1100
1101!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1102
1103end module test_pm_statest
pm_arrayRange::getRange
Generate minimally-spaced character, integer, real sequences or sequences at fixed intervals of size ...
Definition: pm_arrayRange.F90:135
pm_arraySort::setSorted
Sort the input scalar string or contiguous vector in ascending order, or return the sorted indices of...
Definition: pm_arraySort.F90:5540
pm_distNorm::getNormCDF
Generate and return the Cumulative Distribution Function (CDF) of the univariate Normal distribution.
Definition: pm_distNorm.F90:754
pm_distUnifEll::setUnifEllRand
Return a (collection) of random vector(s) of size ndim from the ndim-dimensional MultiVariate Uniform...
Definition: pm_distUnifEll.F90:2714
pm_distanceEuclid::getDisEuclid
Generate and return the (squared) Euclidean distance of a (set of) point(s) in ndim-dimensions from a...
Definition: pm_distanceEuclid.F90:401
pm_distanceEuclid::setDisMatEuclid
Return the full or a subset of the Euclidean (squared) distance matrix of the input set of npnt point...
Definition: pm_distanceEuclid.F90:2663
pm_ellipsoid::getLogVolUnitBall
Generate and return the natural logarithm of the volume of an -dimensional ball of unit-radius.
Definition: pm_ellipsoid.F90:169
pm_knn::setKnnSorted
Return the input distance matrix whose columns are sorted in ascending order on output,...
Definition: pm_knn.F90:181
pm_mathCumSum::setCumSum
Return the cumulative sum of the input array, optionally in the backward direction and optionally,...
Definition: pm_mathCumSum.F90:423
pm_matrixDet::setMatDetSqrtLog
Return the natural logarithm of the square-root of the determinant of the input positive-definite squ...
Definition: pm_matrixDet.F90:1747
pm_matrixInit::getMatInit
Generate and return a matrix of shape (shape(1), shape(2)) with the upper/lower triangle and the diag...
Definition: pm_matrixInit.F90:289
pm_statest::getProbKS
Generate and return the probability of the null-hypothesis that sample1 of size nsam1 originates from...
Definition: pm_statest.F90:196
pm_sysPath::getDirCurrent
Generate and return the Current Working Directory (CWD) of the runtime shell.
Definition: pm_sysPath.F90:2759
pm_val2str::getStr
Generate and return the conversion of the input value to an output Fortran string,...
Definition: pm_val2str.F90:167
pm_arrayRange
This module contains procedures and generic interfaces for generating ranges of discrete character,...
Definition: pm_arrayRange.F90:36
pm_arraySort
This module contains procedures and generic interfaces for various sorting tasks.
Definition: pm_arraySort.F90:89
pm_distNorm
This module contains classes and procedures for computing various statistical quantities related to t...
Definition: pm_distNorm.F90:179
pm_distUnifEll
This module contains classes and procedures for computing various statistical quantities related to t...
Definition: pm_distUnifEll.F90:144
pm_distanceEuclid
This module contains procedures and generic interfaces for computing the Euclidean norm of a single p...
Definition: pm_distanceEuclid.F90:71
pm_distanceEuclid::euclidu
type(euclidu_type), parameter euclidu
This is a scalar parameter object of type euclidu_typethat is exclusively used to request unsafe meth...
Definition: pm_distanceEuclid.F90:207
pm_distanceEuclid::euclidsq
type(euclidsq_type), parameter euclidsq
This is a scalar parameter object of type euclidsq_typethat is exclusively used to request computing ...
Definition: pm_distanceEuclid.F90:271
pm_ellipsoid
This module contains classes and procedures for setting up and computing the properties of the hyper-...
Definition: pm_ellipsoid.F90:105
pm_err
This module contains classes and procedures for reporting and handling errors.
Definition: pm_err.F90:52
pm_kind
This module defines the relevant Fortran kind type-parameters frequently used in the ParaMonte librar...
Definition: pm_kind.F90:268
pm_kind::RK
integer, parameter RK
The default real kind in the ParaMonte library: real64 in Fortran, c_double in C-Fortran Interoperati...
Definition: pm_kind.F90:543
pm_kind::LK
integer, parameter LK
The default logical kind in the ParaMonte library: kind(.true.) in Fortran, kind(....
Definition: pm_kind.F90:541
pm_kind::IK
integer, parameter IK
The default integer kind in the ParaMonte library: int32 in Fortran, c_int32_t in C-Fortran Interoper...
Definition: pm_kind.F90:540
pm_kind::RKS
integer, parameter RKS
The single-precision real kind in Fortran mode. On most platforms, this is an 32-bit real kind.
Definition: pm_kind.F90:567
pm_knn
This module contains procedures and generic interfaces for computing the nearest neighbor statistics ...
Definition: pm_knn.F90:82
pm_mathCumSum
This module contains the procedures and interfaces for computing the cumulative sum of an array.
Definition: pm_mathCumSum.F90:48
pm_matrixDet
This module contains procedures and generic interfaces relevant to the computation of the determinant...
Definition: pm_matrixDet.F90:95
pm_matrixInit
This module contains procedures and generic interfaces relevant to generating and initializing matric...
Definition: pm_matrixInit.F90:116
pm_statest
This module contains classes and procedures for performing various statistical tests.
Definition: pm_statest.F90:77
pm_statest::MODULE_NAME
character(*, SK), parameter MODULE_NAME
Definition: pm_statest.F90:85
pm_sysPath
This module contains classes and procedures for manipulating system file/folder paths.
Definition: pm_sysPath.F90:274
pm_test
This module contains a simple unit-testing framework for the Fortran libraries, including the ParaMon...
Definition: pm_test.F90:42
pm_timer
This module contains the timer procedures and derived types to facilitate timing applications at runt...
Definition: pm_timer.F90:99
pm_val2str
This module contains the generic procedures for converting values of different types and kinds to For...
Definition: pm_val2str.F90:58
test_pm_statest
This module contains tests of the module pm_statest.
Definition: test_pm_hypoTest.F90:20
test_pm_statest::SampleSize
integer(IK), dimension(*), parameter SampleSize
Definition: test_pm_hypoTest.F90:389
test_pm_statest::lenSampleSize
integer(IK), parameter lenSampleSize
Definition: test_pm_hypoTest.F90:401
test_pm_statest::test_getSampleDisparity_2
logical(LK) function test_getSampleDisparity_2()
Any value of SampleSize smaller than one or larger than the size of the input Point should lead to an...
Definition: test_pm_hypoTest.F90:737
test_pm_statest::test_performance_1
logical(LK) function test_performance_1()
Compare performance of getCumDenComKS()` with pm_statest::getProbKS(). Results indicated that pm_stat...
Definition: test_pm_hypoTest.F90:548
test_pm_statest::test_doSortedKS2_1
logical(LK) function test_doSortedKS2_1()
Definition: test_pm_hypoTest.F90:637
test_pm_statest::StdNormRnd2
real(RK), dimension(lenRnd) StdNormRnd2
Definition: test_pm_hypoTest.F90:136
test_pm_statest::lenRnd
integer(IK), parameter lenRnd
Definition: test_pm_hypoTest.F90:32
test_pm_statest::setTest
subroutine setTest()
Definition: test_pm_hypoTest.F90:434
test_pm_statest::test_getCumDenComKS_1
logical(LK) function test_getCumDenComKS_1()
Definition: test_pm_hypoTest.F90:515
test_pm_statest::lenDistSortedDiff
integer(IK), parameter lenDistSortedDiff
Definition: test_pm_hypoTest.F90:187
test_pm_statest::StatKS_REF
real(RK), dimension(*), parameter StatKS_REF
Definition: test_pm_hypoTest.F90:414
test_pm_statest::test_Uniformity
subroutine test_Uniformity()
Definition: test_pm_hypoTest.F90:794
test_pm_statest::test_doKS1_1
logical(LK) function test_doKS1_1()
Definition: test_pm_hypoTest.F90:452
test_pm_statest::StdNormRnd1
real(RK), dimension(lenRnd) StdNormRnd1
Definition: test_pm_hypoTest.F90:85
test_pm_statest::test
type(test_type) test
Definition: test_pm_hypoTest.F90:30
test_pm_statest::UnifRnd
real(RK), dimension(lenRnd) UnifRnd
Definition: test_pm_hypoTest.F90:34
test_pm_statest::test_getSampleDisparity_1
logical(LK) function test_getSampleDisparity_1()
Any value of SampleSize smaller than one or larger than the size of the input Point should lead to an...
Definition: test_pm_hypoTest.F90:703
test_pm_statest::test_doUniformKS1_1
logical(LK) function test_doUniformKS1_1()
Definition: test_pm_hypoTest.F90:584
test_pm_statest::DistSortedDiff
real(RK), dimension(*), parameter DistSortedDiff
Definition: test_pm_hypoTest.F90:188
test_pm_statest::ProbKS_REF
real(RK), dimension(*), parameter ProbKS_REF
Definition: test_pm_hypoTest.F90:402
pm_err::err_type
This is the derived type for generating objects to gracefully and verbosely handle runtime unexpected...
Definition: pm_err.F90:157
pm_test::test_type
This is the derived type test_type for generating objects that facilitate testing of a series of proc...
Definition: pm_test.F90:209
pm_timer::timerSYS_type
This is the timerSYS_type class, containing attributes and static methods for setting up a timer base...
Definition: pm_timer.F90:502
doCrossKS
subroutine doCrossKS()
Definition: test_pm_hypoTest.F90:950
getCDF
PURE real(RK) function getCDF(x)
Definition: test_pm_hypoTest.F90:504
doAutoKS
subroutine doAutoKS()
Definition: test_pm_hypoTest.F90:1058