ParaMonte Fortran 2.0.0
Parallel Monte Carlo and Machine Learning Library
See the latest version documentation.
pm_mathBeta::getBetaInv Interface Reference

Generate and return the regularized Inverse Incomplete Beta Function \(I_x(\alpha, \beta)\) as defined in the details section of pm_mathBeta. More...

Detailed Description

Generate and return the regularized Inverse Incomplete Beta Function \(I_x(\alpha, \beta)\) as defined in the details section of pm_mathBeta.

Parameters
[in]betaInc: The input scalar (or array of the same shape as other array-like arguments) of type real of kind any supported by the processor (e.g., RK, RK32, RK64, or RK128).
[in]alpha: The input scalar (or array of the same shape as other array-like arguments) of the same type and kind as betaInc.
[in]beta: The input scalar (or array of the same shape as other array-like arguments) of the same type and kind as betaInc.
[in]signed: The input scalar (or array of the same shape as other array-like arguments) of type logical of default kind LK.
  1. If signed = .false., the input betaInc must be in range 0 <= betaInc <= 1 and the output betaInv will be the expected incomplete Beta function in range 0 <= betaInv <= 1.
  2. If signed = .true., then the following rules hold:
    1. If the condition betaInc < 0 holds, then the value betaInc = 1 - betaInc < 0 will be used instead of betaInc.
    2. If the output betaInv is near 1, the output will be returned as betaInv = betaInv - 1 < 0 instead of betaInv.
      Therefore, the user is expected to be aware of the convention and apply the necessary correction (addition by 1) before using the output value.
Specifying signed = .true. can lead to considerably more accurate calculations (by orders of magnitudes) for values of betaInc and betaInv that are near 1.
The loss of precision near 1 occurs because of inadequate resolution of real number representations in digital computers near 1 which is orders of magnitude worse than the precision near 0.
(optional, default = .false., following the principle of least surprise.)
Returns
betaInv : The output object of the same type, kind, and rank as highest-rank input argument containing the regularized Inverse Incomplete Beta Function.


Possible calling interfaces ⛓

use pm_mathBeta, only: getBetaInv
betaInv = getBetaInv(betaInc, alpha, beta, signed = signed)
pm_mathBeta::getBetaInv
Generate and return the regularized Inverse Incomplete Beta Function as defined in the details secti...
Definition: pm_mathBeta.F90:678
pm_mathBeta
This module contains classes and procedures for computing the mathematical Beta Function and its inve...
Definition: pm_mathBeta.F90:84
Warning
The conditions 0 <= betaInc .and. 0 <= betaInc must hold for the corresponding input arguments.
The conditions 0 < alpha .and. 0 < beta must hold for the corresponding input arguments.
These conditions are verified only if the library is built with the preprocessor macro CHECK_ENABLED=1.
The procedures under this generic interface will abort the program if the computation of the continued fraction representation of the regularized beta function fails to converge.
The pure procedure(s) documented herein become impure when the ParaMonte library is compiled with preprocessor macro CHECK_ENABLED=1.
By default, these procedures are pure in release build and impure in debug and testing builds.
Remarks
The procedures under discussion are elemental.
See also
getLogBeta
getBetaInv
getBetaLogPDF


Example usage ⛓

1program example
2
3 use pm_kind, only: SK, IK, LK
4 use pm_kind, only: RKG => RK ! all processor kinds are supported.
5 use pm_io, only: display_type
6 use pm_mathBeta, only: getBetaInv
7 use pm_distUnif, only: getUnifRand
8 use pm_mathBeta, only: getBetaInc
9
10 implicit none
11
12 real(RKG), allocatable :: alpha, beta, betaInc(:), betaInv(:)
13
14 type(display_type) :: disp
15 disp = display_type(file = "main.out.F90")
16
17 call disp%skip()
18 call disp%show("betaInc = [0.5_RKG]")
19 betaInc = [0.5_RKG]
20 call disp%show("alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)")
21 alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
22 call disp%show("[alpha, beta]")
23 call disp%show( [alpha, beta] )
24 call disp%show("betaInv = getBetaInv(betaInc, alpha, beta)")
25 betaInv = getBetaInv(betaInc, alpha, beta)
26 call disp%show("betaInv")
27 call disp%show( betaInv )
28 call disp%show("getBetaInc(betaInv, alpha, beta)")
29 call disp%show( getBetaInc(betaInv, alpha, beta) )
30 call disp%skip()
31
32 call disp%skip()
33 call disp%show("betaInc = [0.5_RKG]")
34 betaInc = [0.5_RKG]
35 call disp%show("alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)")
36 alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
37 call disp%show("[alpha, beta]")
38 call disp%show( [alpha, beta] )
39 call disp%show("betaInv = getBetaInv(betaInc, alpha, beta)")
40 betaInv = getBetaInv(betaInc, alpha, beta)
41 call disp%show("betaInv")
42 call disp%show( betaInv )
43 call disp%show("getBetaInc(betaInv, alpha, beta)")
44 call disp%show( getBetaInc(betaInv, alpha, beta) )
45 call disp%skip()
46
47 call disp%skip()
48 call disp%show("betaInc = [1._RKG - epsilon(1._RKG)]")
49 betaInc = [1._RKG - epsilon(1._RKG)]
50 call disp%show("alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)")
51 alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
52 call disp%show("[alpha, beta]")
53 call disp%show( [alpha, beta] )
54 call disp%show("betaInv = getBetaInv(betaInc, alpha, beta)")
55 betaInv = getBetaInv(betaInc, alpha, beta)
56 call disp%show("betaInv")
57 call disp%show( betaInv )
58 call disp%show("getBetaInc(betaInv, alpha, beta)")
59 call disp%show( getBetaInc(betaInv, alpha, beta) )
60 call disp%skip()
61
62 call disp%skip()
63 call disp%show("betaInc = [0._RKG, .5_RKG, 1._RKG]")
64 betaInc = [0._RKG, .5_RKG, 1._RKG]
65 call disp%show("alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)")
66 alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
67 call disp%show("[alpha, beta]")
68 call disp%show( [alpha, beta] )
69 call disp%show("betaInv = getBetaInv(betaInc, alpha, beta)")
70 betaInv = getBetaInv(betaInc, alpha, beta)
71 call disp%show("betaInv")
72 call disp%show( betaInv )
73 call disp%show("getBetaInc(betaInv, alpha, beta)")
74 call disp%show( getBetaInc(betaInv, alpha, beta) )
75 call disp%skip()
76
77 call disp%skip()
78 call disp%show("betaInc = [0._RKG, .5_RKG, 1._RKG]")
79 betaInc = [0._RKG, .5_RKG, 1._RKG]
80 call disp%show("betaInv = getBetaInv(betaInc, alpha = [.1_RKG, 1._RKG, 10._RKG], beta = 3._RKG)")
81 betaInv = getBetaInv(betaInc, alpha = [.1_RKG, 1._RKG, 10._RKG], beta = 3._RKG)
82 call disp%show("betaInv")
83 call disp%show( betaInv )
84 call disp%show("getBetaInc(betaInv, alpha = [.1_RKG, 1._RKG, 10._RKG], beta = 3._RKG)")
85 call disp%show( getBetaInc(betaInv, alpha = [.1_RKG, 1._RKG, 10._RKG], beta = 3._RKG) )
86 call disp%skip()
87
88 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
89 ! Output an example array for visualization.
90 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
91
92 block
93
94 use pm_kind, only: RKG => RKH, RKH
95 use pm_arraySpace, only: setLinSpace
96 integer(IK) , parameter :: NP = 1000
97 real(RKG) , parameter :: alpha(*) = [0.1_RKG, 10._RKG, 1._RKG, 0.1_RKG, 10._RKG], beta(*) = [0.1_RKG, 0.1_RKG, 1._RKG, 10._RKG, 10._RKG]
98 real(RKG) :: betaInv(max(size(alpha), size(beta)))
99 real(RKG) :: betaInvRef(size(betaInv))
100 real(RKG) :: betaInc(NP)
101 integer :: fileUnit, i, j
102
103 call setLinSpace(betaInc, 0._RKG, 1._RKG, fopen = .true._LK, lopen = .true._LK)
104 open(newunit = fileUnit, file = "getBetaInv.RK.txt")
105 do i = 1, NP
106 betaInv = getBetaInv(betaInc(i), alpha, beta)
107 write(fileUnit, "(*(g0,:,','))") betaInc(i), betaInv
108 end do
109 close(fileUnit)
110
111 block
112 use pm_kind, only: RKG => RKS
113 character(*), parameter :: RKSTR = "RKS"
114 real(RKG) :: betaInv(max(size(alpha), size(beta)))
115 open(newunit = fileUnit, file = "getBetaInv."//RKSTR//".abserr.txt")
116 do i = 1, NP
117 betaInv = getBetaInv(real(betaInc(i), RKG), real(alpha, RKG), real(beta, RKG))
118 betaInvRef = getBetaInv(betaInc(i), alpha, beta, signed = .true._LK)
119 write(fileUnit, "(*(g0,:,','))") betaInc(i), abs(betaInv - merge(1 + betaInvRef, betaInvRef, betaInvRef < 0))
120 end do
121 close(fileUnit)
122 end block
123
124 block
125 use pm_kind, only: RKG => RKD
126 character(*), parameter :: RKSTR = "RKD"
127 real(RKG) :: betaInv(max(size(alpha), size(beta)))
128 open(newunit = fileUnit, file = "getBetaInv."//RKSTR//".abserr.txt")
129 do i = 1, NP
130 betaInv = getBetaInv(real(betaInc(i), RKG), real(alpha, RKG), real(beta, RKG))
131 betaInvRef = getBetaInv(betaInc(i), alpha, beta, signed = .true._LK)
132 write(fileUnit, "(*(g0,:,','))") betaInc(i), abs(betaInv - merge(1 + betaInvRef, betaInvRef, betaInvRef < 0))
133 end do
134 close(fileUnit)
135 end block
136
137 block
138 use pm_kind, only: RKG => RKH
139 character(*), parameter :: RKSTR = "RKH"
140 real(RKG) :: betaInv(max(size(alpha), size(beta)))
141 open(newunit = fileUnit, file = "getBetaInv."//RKSTR//".abserr.txt")
142 do i = 1, NP
143 betaInv = getBetaInv(real(betaInc(i), RKG), real(alpha, RKG), real(beta, RKG))
144 betaInvRef = getBetaInv(betaInc(i), alpha, beta, signed = .true._LK)
145 write(fileUnit, "(*(g0,:,','))") betaInc(i), abs(betaInv - merge(1 + betaInvRef, betaInvRef, betaInvRef < 0))
146 end do
147 close(fileUnit)
148 end block
149
150 end block
151
152end program example
pm_arraySpace::setLinSpace
Return the linSpace output argument with size(linSpace) elements of evenly-spaced values over the int...
Definition: pm_arraySpace.F90:324
pm_distUnif::getUnifRand
Generate and return a scalar or a contiguous array of rank 1 of length s1 of randomly uniformly distr...
Definition: pm_distUnif.F90:4150
pm_io::show
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11726
pm_io::skip
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11508
pm_mathBeta::getBetaInc
Generate and return the regularized Incomplete Beta Function as defined in the details section of pm...
Definition: pm_mathBeta.F90:280
pm_arraySpace
This module contains procedures and generic interfaces for generating arrays with linear or logarithm...
Definition: pm_arraySpace.F90:33
pm_distUnif
This module contains classes and procedures for computing various statistical quantities related to t...
Definition: pm_distUnif.F90:274
pm_io
This module contains classes and procedures for input/output (IO) or generic display operations on st...
Definition: pm_io.F90:252
pm_io::disp
type(display_type) disp
This is a scalar module variable an object of type display_type for general display.
Definition: pm_io.F90:11393
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::RKD
integer, parameter RKD
The double precision real kind in Fortran mode. On most platforms, this is an 64-bit real kind.
Definition: pm_kind.F90:568
pm_kind::SK
integer, parameter SK
The default character kind in the ParaMonte library: kind("a") in Fortran, c_char in C-Fortran Intero...
Definition: pm_kind.F90:539
pm_kind::RKH
integer, parameter RKH
The scalar integer constant of intrinsic default kind, representing the highest-precision real kind t...
Definition: pm_kind.F90:858
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_io::display_type
Generate and return an object of type display_type.
Definition: pm_io.F90:10282

Example Unix compile command via Intel ifort compiler ⛓
1#!/usr/bin/env sh
2rm main.exe
3ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example Windows Batch compile command via Intel ifort compiler ⛓
1del main.exe
2set PATH=..\..\..\lib;%PATH%
3ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
4main.exe

Example Unix / MinGW compile command via GNU gfortran compiler ⛓
1#!/usr/bin/env sh
2rm main.exe
3gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example output ⛓
1
2betaInc = [0.5_RKG]
3alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
4[alpha, beta]
5+2.6418922610648354, +4.8484087678096186
6betaInv = getBetaInv(betaInc, alpha, beta)
7betaInv
8+0.33899258858751102
9getBetaInc(betaInv, alpha, beta)
10+0.50000000000000056
11
12
13betaInc = [0.5_RKG]
14alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
15[alpha, beta]
16+4.9099244464307379, +5.5678273004655772
17betaInv = getBetaInv(betaInc, alpha, beta)
18betaInv
19+0.46653360996758036
20getBetaInc(betaInv, alpha, beta)
21+0.49999999999999989
22
23
24betaInc = [1._RKG - epsilon(1._RKG)]
25alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
26[alpha, beta]
27+3.9737188053062686, +7.1403683909559597
28betaInv = getBetaInv(betaInc, alpha, beta)
29betaInv
30+0.99671903397343276
31getBetaInc(betaInv, alpha, beta)
32+0.99999999999999978
33
34
35betaInc = [0._RKG, .5_RKG, 1._RKG]
36alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
37[alpha, beta]
38+6.0549088137146487, +8.8019024201242235
39betaInv = getBetaInv(betaInc, alpha, beta)
40betaInv
41+0.0000000000000000, +0.40329626975593230, +1.0000000000000000
42getBetaInc(betaInv, alpha, beta)
43+0.0000000000000000, +0.49999999999999861, +1.0000000000000000
44
45
46betaInc = [0._RKG, .5_RKG, 1._RKG]
47betaInv = getBetaInv(betaInc, alpha = [.1_RKG, 1._RKG, 10._RKG], beta = 3._RKG)
48betaInv
49+0.0000000000000000, +0.20629947401590026, +1.0000000000000000
50getBetaInc(betaInv, alpha = [.1_RKG, 1._RKG, 10._RKG], beta = 3._RKG)
51+0.0000000000000000, +0.50000000000000011, +1.0000000000000000
52
53

Postprocessing of the example output ⛓
1#!/usr/bin/env python
2
3import matplotlib.pyplot as plt
4import pandas as pd
5import numpy as np
6import glob
7import sys
8import os
9
10fontsize = 17
11
12kind = "RK"
13label = [ r"$\alpha, \beta = .1, .1$"
14 , r"$\alpha, \beta = 10, .1$"
15 , r"$\alpha, \beta = 1., 1.$"
16 , r"$\alpha, \beta = .1, 10$"
17 , r"$\alpha, \beta = 10, 10$"
18 ]
19
20parent = os.path.basename(os.path.dirname(__file__))
21pattern = parent + "*.txt"
22
23fileList = glob.glob(pattern)
24for file in fileList:
25
26 df = pd.read_csv(file, delimiter = ",")
27
28 fig = plt.figure(figsize = 1.25 * np.array([6.4, 4.8]), dpi = 200)
29 ax = plt.subplot()
30
31 for i in range(1,len(df.values[0,:]+1)):
32
33 plt.plot( df.values[:, 0]
34 , df.values[:,i]
35 , linewidth = 2
36 )
37
38 plt.xticks(fontsize = fontsize - 2)
39 plt.yticks(fontsize = fontsize - 2)
40 if "abserr" in file:
41 nbit = file.split(".")[1][2:]
42 ax.set_xlabel("X : Regularized Incomplete Beta Function", fontsize = fontsize)
43 ax.set_ylabel("{}-bits Absolute Error: X - BetaInc(BetaInv(X))".format(nbit), fontsize = fontsize)
44 else:
45 ax.set_xlabel("x", fontsize = fontsize)
46 ax.set_ylabel("Regularized Inverse Beta Function", fontsize = fontsize)
47
48 ax.legend ( label
49 , fontsize = fontsize
50 #, loc = "center left"
51 #, bbox_to_anchor = (1, 0.5)
52 )
53
54 plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
55 ax.tick_params(axis = "y", which = "minor")
56 ax.tick_params(axis = "x", which = "minor")
57 plt.tight_layout()
58
59 plt.savefig(file.replace(".txt",".png"))

Visualization of the example output ⛓
Test:
test_pm_mathBeta


Final Remarks ⛓

If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

  1. If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
  2. If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Copyright
Computational Data Science Lab
Author:
Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan

Definition at line 678 of file pm_mathBeta.F90.

The documentation for this interface was generated from the following file: