statistics.f
1 ! @section Synopsis
2 !
3 ! Statistical functions for VIF.
4 !
5 ! @author J. Marcel van der Veer
6 !
7 ! @section copyright
8 !
9 ! This file is part of VIF - vintage fortran compiler.
10 ! Copyright 2020-2025 J. Marcel van der Veer <algol68g@xs4all.nl>.
11 !
12 ! @section license
13 !
14 ! This program is free software; you can redistribute it and/or modify it
15 ! under the terms of the gnu general public license as published by the
16 ! free software foundation; either version 3 of the license, or
17 ! (at your option) any later version.
18 !
19 ! This program is distributed in the hope that it will be useful, but
20 ! without any warranty; without even the implied warranty of merchantability
21 ! or fitness for a particular purpose. See the GNU general public license for
22 ! more details. You should have received a copy of the GNU general public
23 ! license along with this program. If not, see <http://www.gnu.org/licenses/>.
24 !
25
26 REAL*8 FUNCTION dlnorm (x, upper)
27 !
28 ! DLNORM computes the cumulative density of the standard normal distribution.
29 !
30 ! Modified:
31 !
32 ! 28 March 1999
33 !
34 ! Author:
35 !
36 ! David Hill
37 ! Modifications by John Burkardt
38 !
39 ! Reference:
40 !
41 ! David Hill,
42 ! Algorithm AS 66:
43 ! The Normal Integral,
44 ! Applied Statistics,
45 ! Volume 22, Number 3, 1973, pages 424-427.
46 !
47 ! Parameters:
48 !
49 ! Input, real*8 X, is one endpoint of the semi-infinite interval
50 ! over which the integration takes place.
51 !
52 ! Input, logical UPPER, determines whether the upper or lower
53 ! interval is to be integrated:
54 ! .TRUE. => integrate from X to + Infinity;
55 ! .FALSE. => integrate from - Infinity to X.
56 !
57 ! Output, real*8 ALNORM, the integral of the standard normal
58 ! distribution over the desired interval.
59 !
60 IMPLICIT none
61
62 REAL*8 a1
63 PARAMETER (a1 = 5.75885480458d+00)
64 REAL*8 a2
65 PARAMETER (a2 = 2.62433121679d+00)
66 REAL*8 a3
67 PARAMETER (a3 = 5.92885724438d+00)
68 REAL*8 b1
69 PARAMETER (b1 = -29.8213557807d+00)
70 REAL*8 b2
71 PARAMETER (b2 = 48.6959930692d+00)
72 REAL*8 c1
73 PARAMETER (c1 = -0.000000038052d+00)
74 REAL*8 c2
75 PARAMETER (c2 = 0.000398064794d+00)
76 REAL*8 c3
77 PARAMETER (c3 = -0.151679116635d+00)
78 REAL*8 c4
79 PARAMETER (c4 = 4.8385912808d+00)
80 REAL*8 c5
81 PARAMETER (c5 = 0.742380924027d+00)
82 REAL*8 c6
83 PARAMETER (c6 = 3.99019417011d+00)
84 REAL*8 con
85 PARAMETER (con = 1.28d+00)
86 REAL*8 d1
87 PARAMETER (d1 = 1.00000615302d+00)
88 REAL*8 d2
89 PARAMETER (d2 = 1.98615381364d+00)
90 REAL*8 d3
91 PARAMETER (d3 = 5.29330324926d+00)
92 REAL*8 d4
93 PARAMETER (d4 = -15.1508972451d+00)
94 REAL*8 d5
95 PARAMETER (d5 = 30.789933034d+00)
96 REAL*8 ltone
97 PARAMETER (ltone = 7.0d+00)
98 REAL*8 p
99 PARAMETER (p = 0.398942280444d+00)
100 REAL*8 q
101 PARAMETER (q = 0.39990348504d+00)
102 REAL*8 r
103 PARAMETER (r = 0.398942280385d+00)
104 LOGICAL up
105 LOGICAL upper
106 REAL*8 utzero
107 PARAMETER (utzero = 18.66d+00)
108 REAL*8 x
109 REAL*8 y
110 REAL*8 z
111
112 up = upper
113 z = x
114
115 IF (z .lt. 0.0d+00) THEN
116 up = .not. up
117 z = - z
118 END IF
119
120 IF (ltone .lt. z .and.
121 & ((.not. up) .or. utzero .lt. z)) THEN
122
123 IF (up) THEN
124 dlnorm = 0.0d+00
125 ELSE
126 dlnorm = 1.0d+00
127 END IF
128
129 RETURN
130
131 END IF
132
133 y = 0.5d+00 * z * z
134
135 IF (z .le. con) THEN
136
137 dlnorm = 0.5d+00 - z * (p - q * y
138 & / (y + a1 + b1
139 & / (y + a2 + b2
140 & / (y + a3))))
141
142 ELSE
143
144 dlnorm = r * dexp (- y)
145 & / (z + c1 + d1
146 & / (z + c2 + d2
147 & / (z + c3 + d3
148 & / (z + c4 + d4
149 & / (z + c5 + d5
150 & / (z + c6))))))
151
152 END IF
153
154 IF (.not. up) THEN
155 dlnorm = 1.0d+00 - dlnorm
156 END IF
157
158 RETURN
159 END
160
161 REAL*8 FUNCTION prncst (st, idf, d, ifault)
162
163 !
164 ! PRNCST computes the lower tail of noncentral T distribution.
165 !
166 ! Modified:
167 !
168 ! 06 January 2008
169 !
170 ! Author:
171 !
172 ! BE Cooper
173 ! Modifications by John Burkardt
174 !
175 ! Reference:
176 !
177 ! BE Cooper,
178 ! Algorithm AS 5:
179 ! The Integral of the Non-Central T-Distribution,
180 ! Applied Statistics,
181 ! Volume 17, Number 2, 1968, page 193.
182 !
183 ! Parameters:
184 !
185 ! Input, real*8 ST, the argument.
186 !
187 ! Input, integer IDF, the number of degrees of freedom.
188 !
189 ! Input, real*8 D, the noncentrality parameter.
190 !
191 ! Output, integer IFAULT, error flag.
192 ! 0, no error occurred.
193 ! nonzero, an error occurred.
194 !
195 ! Output, real*8 PRNCST, the value of the lower tail of
196 ! the noncentral T distribution.
197 !
198 ! Local Parameters:
199 !
200 ! Local, real*8 G1, 1.0 / sqrt(2.0 * pi)
201 !
202 ! Local, real*8 G2, 1.0 / (2.0 * pi)
203 !
204 ! Local, real*8 G3, sqrt(2.0 * pi)
205 !
206 IMPLICIT none
207
208 REAL*8 a
209 REAL*8 ak
210 REAL*8 dlnorm
211 REAL*8 b
212 REAL*8 d
213 REAL*8 da
214 REAL*8 drb
215 REAL*8 emin
216 PARAMETER (emin = 12.5d+00)
217 REAL*8 f
218 REAL*8 fk
219 REAL*8 fkm1
220 REAL*8 fmkm1
221 REAL*8 fmkm2
222 REAL*8 g1
223 PARAMETER (g1 = 0.3989422804d+00)
224 REAL*8 g2
225 PARAMETER (g2 = 0.1591549431d+00)
226 REAL*8 g3
227 PARAMETER (g3 = 2.5066282746d+00)
228 INTEGER idf
229 INTEGER ifault
230 INTEGER ioe
231 INTEGER k
232 REAL*8 rb
233 REAL*8 st
234 REAL*8 sum
235 REAL*8 tfn
236
237 f = dble (idf)
238 !
239 ! For very large IDF, use the normal approximation.
240 !
241 IF (100 .lt. idf) THEN
242
243 ifault = 1
244
245 a = dsqrt (0.5d+00 * f)
246 & * dexp (dlngam (0.5d+00 * (f - 1.0d+00))
247 & - dlngam (0.5d+00 * f)) * d
248
249 prncst = dlnorm ((st - a) / dsqrt (f * (1.0d+00 + d * d)
250 & / (f - 2.0d+00) - a * a), .false.)
251
252 RETURN
253 END IF
254
255 ifault = 0
256 ioe = mod (idf, 2)
257 a = st / dsqrt (f)
258 b = f / (f + st * st)
259 rb = dsqrt (b)
260 da = d * a
261 drb = d * rb
262
263 IF (idf .eq. 1) THEN
264 prncst = dlnorm (drb, .true.) + 2.0d+00 * tfn (drb, a)
265 RETURN
266 END IF
267
268 sum = 0.0d+00
269
270 IF (dabs (drb) .lt. emin) THEN
271 fmkm2 = a * rb * dexp (- 0.5d+00 * drb * drb)
272 & * dlnorm (a * drb, .false.) * g1
273 ELSE
274 fmkm2 = 0.0d+00
275 END IF
276
277 fmkm1 = b * da * fmkm2
278
279 IF (dabs (d) .lt. emin) THEN
280 fmkm1 = fmkm1 + b * a * g2 * dexp (- 0.5d+00 * d * d)
281 END IF
282
283 IF (ioe .eq. 0) THEN
284 sum = fmkm2
285 ELSE
286 sum = fmkm1
287 END IF
288
289 ak = 1.0d+00
290 fk = 2.0d+00
291
292 DO k = 2, idf - 2, 2
293
294 fkm1 = fk - 1.0d+00
295 fmkm2 = b * (da * ak * fmkm1 + fmkm2) * fkm1 / fk
296 ak = 1.0d+00 / (ak * fkm1)
297 fmkm1 = b * (da * ak * fmkm2 + fmkm1) * fk / (fk + 1.0d+00)
298
299 IF (ioe .eq. 0) THEN
300 sum = sum + fmkm2
301 ELSE
302 sum = sum + fmkm1
303 END IF
304
305 ak = 1.0d+00 / (ak * fk)
306 fk = fk + 2.0d+00
307
308 END DO
309
310 IF (ioe .eq. 0) THEN
311 prncst = dlnorm (d, .true.) + sum * g3
312 ELSE
313 prncst = dlnorm (drb, .true.)
314 & + 2.0d+00 * (sum + tfn (drb, a))
315 END IF
316
317 RETURN
318 END
319
320 SUBROUTINE student_noncentral_cdf_values (n_data, df, lambda,
321 & x, fx)
322
323 !
324 ! STUDENT_NONCENTRAL_CDF_VALUES returns values of the noncentral Student CDF.
325 !
326 ! Discussion:
327 !
328 ! In Mathematica, the function can be evaluated by:
329 !
330 ! Needs["Statistics`ContinuousDistributions`"]
331 ! dist = NoncentralStudentTDistribution [ df, lambda ]
332 ! CDF [ dist, x ]
333 !
334 ! Mathematica seems to have some difficulty computing this function
335 ! to the desired number of digits.
336 !
337 ! Modified:
338 !
339 ! 25 March 2007
340 !
341 ! Author:
342 !
343 ! John Burkardt
344 !
345 ! Reference:
346 !
347 ! Milton Abramowitz, Irene Stegun,
348 ! Handbook of Mathematical Functions,
349 ! National Bureau of Standards, 1964,
350 ! ISBN: 0-486-61272-4,
351 ! LC: QA47.A34.
352 !
353 ! Stephen Wolfram,
354 ! The Mathematica Book,
355 ! Fourth Edition,
356 ! Cambridge University Press, 1999,
357 ! ISBN: 0-521-64314-7,
358 ! LC: QA76.95.W65.
359 !
360 ! Parameters:
361 !
362 ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the
363 ! first call. On each call, the routine increments N_DATA by 1, and
364 ! returns the corresponding data; when there is no more data, the
365 ! output value of N_DATA will be 0 again.
366 !
367 ! Output, integer DF, real*8 LAMBDA, the parameters of the
368 ! function.
369 !
370 ! Output, real*8 X, the argument of the function.
371 !
372 ! Output, real*8 FX, the value of the function.
373 !
374 IMPLICIT none
375
376 INTEGER n_max
377 PARAMETER (n_max = 30)
378
379 INTEGER df
380 INTEGER df_vec(n_max)
381 REAL*8 fx
382 REAL*8 fx_vec(n_max)
383 REAL*8 lambda
384 REAL*8 lambda_vec(n_max)
385 INTEGER n_data
386 REAL*8 x
387 REAL*8 x_vec(n_max)
388
389 SAVE df_vec
390 SAVE fx_vec
391 SAVE lambda_vec
392 SAVE x_vec
393
394 DATA df_vec /
395 & 1, 2, 3,
396 & 1, 2, 3,
397 & 1, 2, 3,
398 & 1, 2, 3,
399 & 1, 2, 3,
400 & 15, 20, 25,
401 & 1, 2, 3,
402 & 10, 10, 10,
403 & 10, 10, 10,
404 & 10, 10, 10 /
405 DATA fx_vec /
406 & 0.8975836176504333d+00,
407 & 0.9522670169d+00,
408 & 0.9711655571887813d+00,
409 & 0.8231218864d+00,
410 & 0.9049021510d+00,
411 & 0.9363471834d+00,
412 & 0.7301025986d+00,
413 & 0.8335594263d+00,
414 & 0.8774010255d+00,
415 & 0.5248571617d+00,
416 & 0.6293856597d+00,
417 & 0.6800271741d+00,
418 & 0.20590131975d+00,
419 & 0.2112148916d+00,
420 & 0.2074730718d+00,
421 & 0.9981130072d+00,
422 & 0.9994873850d+00,
423 & 0.9998391562d+00,
424 & 0.168610566972d+00,
425 & 0.16967950985d+00,
426 & 0.1701041003d+00,
427 & 0.9247683363d+00,
428 & 0.7483139269d+00,
429 & 0.4659802096d+00,
430 & 0.9761872541d+00,
431 & 0.8979689357d+00,
432 & 0.7181904627d+00,
433 & 0.9923658945d+00,
434 & 0.9610341649d+00,
435 & 0.8688007350d+00 /
436 DATA lambda_vec /
437 & 0.0d+00,
438 & 0.0d+00,
439 & 0.0d+00,
440 & 0.5d+00,
441 & 0.5d+00,
442 & 0.5d+00,
443 & 1.0d+00,
444 & 1.0d+00,
445 & 1.0d+00,
446 & 2.0d+00,
447 & 2.0d+00,
448 & 2.0d+00,
449 & 4.0d+00,
450 & 4.0d+00,
451 & 4.0d+00,
452 & 7.0d+00,
453 & 7.0d+00,
454 & 7.0d+00,
455 & 1.0d+00,
456 & 1.0d+00,
457 & 1.0d+00,
458 & 2.0d+00,
459 & 3.0d+00,
460 & 4.0d+00,
461 & 2.0d+00,
462 & 3.0d+00,
463 & 4.0d+00,
464 & 2.0d+00,
465 & 3.0d+00,
466 & 4.0d+00 /
467 DATA x_vec /
468 & 3.00d+00,
469 & 3.00d+00,
470 & 3.00d+00,
471 & 3.00d+00,
472 & 3.00d+00,
473 & 3.00d+00,
474 & 3.00d+00,
475 & 3.00d+00,
476 & 3.00d+00,
477 & 3.00d+00,
478 & 3.00d+00,
479 & 3.00d+00,
480 & 3.00d+00,
481 & 3.00d+00,
482 & 3.00d+00,
483 & 15.00d+00,
484 & 15.00d+00,
485 & 15.00d+00,
486 & 0.05d+00,
487 & 0.05d+00,
488 & 0.05d+00,
489 & 4.00d+00,
490 & 4.00d+00,
491 & 4.00d+00,
492 & 5.00d+00,
493 & 5.00d+00,
494 & 5.00d+00,
495 & 6.00d+00,
496 & 6.00d+00,
497 & 6.00d+00 /
498
499 IF (n_data .lt. 0) THEN
500 n_data = 0
501 END IF
502
503 n_data = n_data + 1
504
505 IF (n_max .lt. n_data) THEN
506 n_data = 0
507 df = 0
508 lambda = 0.0d+00
509 x = 0.0d+00
510 fx = 0.0d+00
511 ELSE
512 df = df_vec(n_data)
513 lambda = lambda_vec(n_data)
514 x = x_vec(n_data)
515 fx = fx_vec(n_data)
516 END IF
517
518 RETURN
519 END
520
521 REAL*8 FUNCTION tfn (x, fx)
522 !
523 ! TFN calculates the T-function of Owen.
524 !
525 ! Modified:
526 !
527 ! 06 January 2008
528 !
529 ! Author:
530 !
531 ! JC Young, Christoph Minder
532 ! Modifications by John Burkardt
533 !
534 ! Reference:
535 !
536 ! MA Porter, DJ Winstanley,
537 ! Remark AS R30:
538 ! A Remark on Algorithm AS76:
539 ! An Integral Useful in Calculating Noncentral T and Bivariate
540 ! Normal Probabilities,
541 ! Applied Statistics,
542 ! Volume 28, Number 1, 1979, page 113.
543 !
544 ! JC Young, Christoph Minder,
545 ! Algorithm AS 76:
546 ! An Algorithm Useful in Calculating Non-Central T and
547 ! Bivariate Normal Distributions,
548 ! Applied Statistics,
549 ! Volume 23, Number 3, 1974, pages 455-457.
550 !
551 ! Parameters:
552 !
553 ! Input, real*8 X, FX, the parameters of the function.
554 !
555 ! Output, real*8 TFN, the value of the T-function.
556 !
557 IMPLICIT none
558
559 INTEGER ng
560 PARAMETER (ng = 5)
561
562 REAL*8 fx
563 REAL*8 fxs
564 INTEGER i
565 REAL*8 r(ng)
566 REAL*8 r1
567 REAL*8 r2
568 REAL*8 rt
569 REAL*8 tp
570 PARAMETER (tp = 0.159155d+00)
571 REAL*8 tv1
572 PARAMETER (tv1 = 1.0d-35)
573 REAL*8 tv2
574 PARAMETER (tv2 = 15.0d+00)
575 REAL*8 tv3
576 PARAMETER (tv3 = 15.0d+00)
577 REAL*8 tv4
578 PARAMETER (tv4 = 1.0d-05)
579 REAL*8 u(ng)
580 REAL*8 x
581 REAL*8 x1
582 REAL*8 x2
583 REAL*8 xs
584
585 DATA u /
586 & 0.0744372d+00,
587 & 0.2166977d+00,
588 & 0.3397048d+00,
589 & 0.4325317d+00,
590 & 0.4869533d+00 /
591
592 DATA r /
593 & 0.1477621d+00,
594 & 0.1346334d+00,
595 & 0.1095432d+00,
596 & 0.0747257d+00,
597 & 0.0333357d+00 /
598 !
599 ! Test for X near zero.
600 !
601 IF (dabs (x) .lt. tv1) THEN
602 tfn = tp * atan (fx)
603 RETURN
604 END IF
605 !
606 ! Test for large values of abs(X).
607 !
608 IF (tv2 .lt. dabs (x)) THEN
609 tfn = 0.0d+00
610 RETURN
611 END IF
612 !
613 ! Test for FX near zero.
614 !
615 IF (dabs (fx) .lt. tv1) THEN
616 tfn = 0.0d+00
617 RETURN
618 END IF
619 !
620 ! Test whether abs (FX) is so large that it must be truncated.
621 !
622 xs = - 0.5d+00 * x * x
623 x2 = fx
624 fxs = fx * fx
625
626 IF (dlog (1.0d+00 + fxs) - xs * fxs .lt. tv3) THEN
627 GO TO 2
628 END IF
629 !
630 ! Computation of truncation point by Newton iteration.
631 !
632 x1 = 0.5d+00 * fx
633 fxs = 0.25d+00 * fxs
634
635 1 CONTINUE
636
637 rt = fxs + 1.0d+00
638
639 x2 = x1 + (xs * fxs + tv3 - dlog (rt))
640 & / (2.0d+00 * x1 * (1.0d+00 / rt - xs))
641
642 fxs = x2 * x2
643
644 IF (dabs (x2 - x1) .lt. tv4) THEN
645 GO TO 2
646 END IF
647
648 x1 = x2
649
650 GO TO 1
651 !
652 ! Gaussian quadrature.
653 !
654 2 CONTINUE
655
656 rt = 0.0d+00
657
658 DO i = 1, ng
659
660 r1 = 1.0d+00 + fxs * (0.5d+00 + u(i))**2
661 r2 = 1.0d+00 + fxs * (0.5d+00 - u(i))**2
662
663 rt = rt + r(i) * (dexp (xs * r1) / r1
664 & + dexp (xs * r2) / r2)
665
666 END DO
667
668 tfn = rt * x2 * tp
669
670 RETURN
671 END
|
© 2002-2025 J.M. van der Veer (jmvdveer@xs4all.nl)
|
