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

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

Detailed Description

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

Parameters
[in]betaInv: The output scalar (or array of the same shape as other array-like arguments) of the same type and kind as the input argument betaInc containing the regularized Inverse Incomplete Beta Function.
[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),
containing the value at which the function must be computed.
[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]logFuncBeta: The input scalar (or array of the same shape as other array-like arguments) of the same type and kind as betaInc, representing the natural logarithm of the Beta Function as returned by getLogBeta(alpha, beta).
Providing this argument can lead to significant runtime performance gains when the interface is called repeatedly for different betaInc values but with identical alpha and beta.
[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.
[out]info: The output scalar (or array of the same shape as other array-like arguments) of type integer of default kind IK.
  1. It is set to 0 if and only if the algorithm returns successfully.
  2. It is set to a non-zero if and only if the computation of the inverse Incomplete Beta Function at the specified value fails to converge.
    See the documentation of the corresponding argument info of setBetaInc for possible output values and their meanings.
Returns


Possible calling interfaces

call setBetaInv(betaInv, betaInc, alpha, beta, logFuncBeta, signed, info)
Return the regularized Inverse Incomplete Beta Function as defined in the details section of pm_math...
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. betaInc <= 1 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 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
setBetaInv
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_mathBeta, only: setBetaInv
8 use pm_distUnif, only: getUnifRand
9 use pm_mathBeta, only: getBetaInc
10 use pm_mathBeta, only: getLogBeta
11 use pm_arrayResize, only: setResized
12
13 implicit none
14
15 logical(LK) :: signed
16 integer(IK) :: info, infos(3)
17 real(RKG) :: alpha, beta, betaInc, betaInv
18 real(RKG), allocatable :: alphas(:), betas(:), betaIncs(:), betaInvs(:)
19
20 type(display_type) :: disp
21 disp = display_type(file = "main.out.F90")
22
23 signed = .false.
24
25 call disp%skip()
26 call disp%show("signed")
27 call disp%show( signed )
28 call disp%show("betaInc = 0.5_RKG")
29 betaInc = 0.5_RKG
30 call disp%show("alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)")
31 alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
32 call disp%show("[alpha, beta]")
33 call disp%show( [alpha, beta] )
34 call disp%show("call setBetaInv(betaInv, betaInc, alpha, beta, getLogBeta(alpha, beta), signed, info)")
35 call setBetaInv(betaInv, betaInc, alpha, beta, getLogBeta(alpha, beta), signed, info)
36 call disp%show("if (info /= 0) error stop 'Beta inversion failed.'")
37 if (info /= 0) error stop 'Beta inversion failed.'
38 call disp%show("betaInv")
39 call disp%show( betaInv )
40 call disp%show("getBetaInc(betaInv, alpha, beta, signed)")
41 call disp%show( getBetaInc(betaInv, alpha, beta, signed) )
42 call disp%skip()
43
44 call disp%skip()
45 call disp%show("signed")
46 call disp%show( signed )
47 call disp%show("betaInc = 0.5_RKG")
48 betaInc = 0.5_RKG
49 call disp%show("alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)")
50 alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
51 call disp%show("[alpha, beta]")
52 call disp%show( [alpha, beta] )
53 call disp%show("call setBetaInv(betaInv, betaInc, alpha, beta, getLogBeta(alpha, beta), signed, info)")
54 call setBetaInv(betaInv, betaInc, alpha, beta, getLogBeta(alpha, beta), signed, info)
55 call disp%show("if (info /= 0) error stop 'Beta inversion failed.'")
56 if (info /= 0) error stop 'Beta inversion failed.'
57 call disp%show("betaInv")
58 call disp%show( betaInv )
59 call disp%show("getBetaInc(betaInv, alpha, beta, signed)")
60 call disp%show( getBetaInc(betaInv, alpha, beta, signed) )
61 call disp%skip()
62
63 call disp%skip()
64 call disp%show("signed")
65 call disp%show( signed )
66 call disp%show("betaInc = 1._RKG - epsilon(1._RKG)")
67 betaInc = 1._RKG - epsilon(1._RKG)
68 call disp%show("alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)")
69 alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
70 call disp%show("[alpha, beta]")
71 call disp%show( [alpha, beta] )
72 call disp%show("call setBetaInv(betaInv, betaInc, alpha, beta, getLogBeta(alpha, beta), signed, info)")
73 call setBetaInv(betaInv, betaInc, alpha, beta, getLogBeta(alpha, beta), signed, info)
74 call disp%show("if (info /= 0) error stop 'Beta inversion failed.'")
75 if (info /= 0) error stop 'Beta inversion failed.'
76 call disp%show("betaInv")
77 call disp%show( betaInv )
78 call disp%show("getBetaInc(betaInv, alpha, beta, signed)")
79 call disp%show( getBetaInc(betaInv, alpha, beta, signed) )
80 call disp%skip()
81
82 call disp%skip()
83 call disp%show("signed")
84 call disp%show( signed )
85 call disp%show("betaIncs = [0._RKG, .5_RKG, 1._RKG]")
86 betaIncs = [0._RKG, .5_RKG, 1._RKG]
87 call disp%show("alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)")
88 alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
89 call disp%show("[alpha, beta]")
90 call disp%show( [alpha, beta] )
91 call disp%show("call setResized(betaInvs, size(betaIncs, 1, IK))")
92 call setResized(betaInvs, size(betaIncs, 1, IK))
93 call disp%show("call setBetaInv(betaInvs, betaIncs, alpha, beta, getLogBeta(alpha, beta), signed, infos)")
94 call setBetaInv(betaInvs, betaIncs, alpha, beta, getLogBeta(alpha, beta), signed, infos)
95 call disp%show("if (any(infos /= 0)) error stop 'Beta inversion failed.'")
96 if (any(infos /= 0)) error stop 'Beta inversion failed.'
97 call disp%show("betaInvs")
98 call disp%show( betaInvs )
99 call disp%show("getBetaInc(betaInvs, alpha, beta, signed)")
100 call disp%show( getBetaInc(betaInvs, alpha, beta, signed) )
101 call disp%skip()
102
103 call disp%skip()
104 call disp%show("signed")
105 call disp%show( signed )
106 call disp%show("betaIncs = [0._RKG, .5_RKG, 1._RKG]")
107 betaIncs = [0._RKG, .5_RKG, 1._RKG]
108 call disp%show("alphas = [.1_RKG, 1._RKG, 10._RKG]")
109 alphas = [.1_RKG, 1._RKG, 10._RKG]
110 call disp%show("beta = 3._RKG")
111 beta = 3._RKG
112 call disp%show("call setResized(betaInvs, size(betaIncs, 1, IK))")
113 call setResized(betaInvs, size(betaIncs, 1, IK))
114 call disp%show("call setBetaInv(betaInvs, betaIncs, alphas, beta, getLogBeta(alpha, beta), signed, infos)")
115 call setBetaInv(betaInvs, betaIncs, alphas, beta, getLogBeta(alpha, beta), signed, infos)
116 call disp%show("if (any(infos /= 0)) error stop 'Beta inversion failed.'")
117 if (any(infos /= 0)) error stop 'Beta inversion failed.'
118 call disp%show("betaInvs")
119 call disp%show( betaInvs )
120 call disp%show("getBetaInc(betaInvs, alpha = [.1_RKG, 1._RKG, 10._RKG], beta = 3._RKG, signed = signed)")
121 call disp%show( getBetaInc(betaInvs, alpha = [.1_RKG, 1._RKG, 10._RKG], beta = 3._RKG, signed = signed) )
122 call disp%skip()
123
124 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
125 ! Output an example array for visualization.
126 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
127
128 block
129
130 use pm_kind, only: RKG => RKH, RKH
131 use pm_arraySpace, only: setLinSpace
132 integer(IK) , parameter :: NP = 1000
133 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]
134 real(RKG) :: betaInv(max(size(alpha), size(beta)))
135 real(RKG) :: betaInvRef(size(betaInv))
136 real(RKG) :: betaInc(NP)
137 integer(IK) :: infos(size(betaInv))
138 integer :: fileUnit, i, j
139
140 call setLinSpace(betaInc, 0._RKG, 1._RKG, fopen = .true._LK, lopen = .true._LK)
141 open(newunit = fileUnit, file = "setBetaInv.RK.txt")
142 do i = 1, NP
143 call setBetaInv(betaInv, betaInc(i), alpha, beta, getLogBeta(alpha, beta), signed, infos)
144 if (any(infos /= 0)) error stop "Beta inversion failed."
145 write(fileUnit, "(*(g0,:,','))") betaInc(i), betaInv
146 end do
147 close(fileUnit)
148
149 block
150 use pm_kind, only: RKG => RKS
151 character(*), parameter :: RKSTR = "RKS"
152 real(RKG) :: betaInv(max(size(alpha), size(beta)))
153 open(newunit = fileUnit, file = "setBetaInv."//RKSTR//".abserr.txt")
154 do i = 1, NP
155 call setBetaInv(betaInv, real(betaInc(i), RKG), real(alpha, RKG), real(beta, RKG), getLogBeta(real(alpha, RKG), real(beta, RKG)), signed, infos)
156 if (any(infos /= 0)) error stop "Beta inversion failed."
157 betaInvRef = getBetaInv(betaInc(i), alpha, beta, signed = .true._LK)
158 write(fileUnit, "(*(g0,:,','))") betaInc(i), abs(betaInv - merge(1 + betaInvRef, betaInvRef, betaInvRef < 0))
159 end do
160 close(fileUnit)
161 end block
162
163 block
164 use pm_kind, only: RKG => RKD
165 character(*), parameter :: RKSTR = "RKD"
166 real(RKG) :: betaInv(max(size(alpha), size(beta)))
167 open(newunit = fileUnit, file = "setBetaInv."//RKSTR//".abserr.txt")
168 do i = 1, NP
169 call setBetaInv(betaInv, real(betaInc(i), RKG), real(alpha, RKG), real(beta, RKG), getLogBeta(real(alpha, RKG), real(beta, RKG)), signed, infos)
170 if (any(infos /= 0)) error stop "Beta inversion failed."
171 betaInvRef = getBetaInv(betaInc(i), alpha, beta, signed = .true._LK)
172 write(fileUnit, "(*(g0,:,','))") betaInc(i), abs(betaInv - merge(1 + betaInvRef, betaInvRef, betaInvRef < 0))
173 end do
174 close(fileUnit)
175 end block
176
177 block
178 use pm_kind, only: RKG => RKH
179 character(*), parameter :: RKSTR = "RKH"
180 real(RKG) :: betaInv(max(size(alpha), size(beta)))
181 open(newunit = fileUnit, file = "setBetaInv."//RKSTR//".abserr.txt")
182 do i = 1, NP
183 call setBetaInv(betaInv, real(betaInc(i), RKG), real(alpha, RKG), real(beta, RKG), getLogBeta(real(alpha, RKG), real(beta, RKG)), signed, infos)
184 if (any(infos /= 0)) error stop "Beta inversion failed."
185 betaInvRef = getBetaInv(betaInc(i), alpha, beta, signed = .true._LK)
186 write(fileUnit, "(*(g0,:,','))") betaInc(i), abs(betaInv - merge(1 + betaInvRef, betaInvRef, betaInvRef < 0))
187 end do
188 close(fileUnit)
189 end block
190
191 end block
192
193end program example
Allocate or resize (shrink or expand) an input allocatable scalar string or array of rank 1....
Return the linSpace output argument with size(linSpace) elements of evenly-spaced values over the int...
Generate and return a scalar or a contiguous array of rank 1 of length s1 of randomly uniformly distr...
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11726
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11508
Generate and return the regularized Incomplete Beta Function as defined in the details section of pm...
Generate and return the regularized Inverse Incomplete Beta Function as defined in the details secti...
Generate and return the natural logarithm of the Beta Function as defined in the details section of ...
This module contains procedures and generic interfaces for resizing allocatable arrays of various typ...
This module contains procedures and generic interfaces for generating arrays with linear or logarithm...
This module contains classes and procedures for computing various statistical quantities related to t...
This module contains classes and procedures for input/output (IO) or generic display operations on st...
Definition: pm_io.F90:252
type(display_type) disp
This is a scalar module variable an object of type display_type for general display.
Definition: pm_io.F90:11393
This module defines the relevant Fortran kind type-parameters frequently used in the ParaMonte librar...
Definition: pm_kind.F90:268
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
integer, parameter LK
The default logical kind in the ParaMonte library: kind(.true.) in Fortran, kind(....
Definition: pm_kind.F90:541
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
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
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
integer, parameter RKH
The scalar integer constant of intrinsic default kind, representing the highest-precision real kind t...
Definition: pm_kind.F90:858
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
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
2signed
3F
4betaInc = 0.5_RKG
5alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
6[alpha, beta]
7+3.2976681223369355, +6.2550693246846247
8call setBetaInv(betaInv, betaInc, alpha, beta, getLogBeta(alpha, beta), signed, info)
9if (info /= 0) error stop 'Beta inversion failed.'
10betaInv
11+0.33400809038700580
12getBetaInc(betaInv, alpha, beta, signed)
13+0.50000000000000111
14
15
16signed
17F
18betaInc = 0.5_RKG
19alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
20[alpha, beta]
21+2.3533391252636435, +4.3429897156975619
22call setBetaInv(betaInv, betaInc, alpha, beta, getLogBeta(alpha, beta), signed, info)
23if (info /= 0) error stop 'Beta inversion failed.'
24betaInv
25+0.33589224526626843
26getBetaInc(betaInv, alpha, beta, signed)
27+0.49999999999999961
28
29
30signed
31F
32betaInc = 1._RKG - epsilon(1._RKG)
33alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
34[alpha, beta]
35+4.9678587351010206, +9.3837541022825572
36call setBetaInv(betaInv, betaInc, alpha, beta, getLogBeta(alpha, beta), signed, info)
37if (info /= 0) error stop 'Beta inversion failed.'
38betaInv
39+0.98940797702415051
40getBetaInc(betaInv, alpha, beta, signed)
41+0.99999999999999978
42
43
44signed
45F
46betaIncs = [0._RKG, .5_RKG, 1._RKG]
47alpha = getUnifRand(0.1_RKG, 10._RKG); beta = getUnifRand(0.1_RKG, 10._RKG)
48[alpha, beta]
49+2.3932013480958121, +2.1735288873594900
50call setResized(betaInvs, size(betaIncs, 1, IK))
51call setBetaInv(betaInvs, betaIncs, alpha, beta, getLogBeta(alpha, beta), signed, infos)
52if (any(infos /= 0)) error stop 'Beta inversion failed.'
53betaInvs
54+0.0000000000000000, +0.52784376386833431, +1.0000000000000000
55getBetaInc(betaInvs, alpha, beta, signed)
56+0.0000000000000000, +0.49999999999999967, +1.0000000000000000
57
58
59signed
60F
61betaIncs = [0._RKG, .5_RKG, 1._RKG]
62alphas = [.1_RKG, 1._RKG, 10._RKG]
63beta = 3._RKG
64call setResized(betaInvs, size(betaIncs, 1, IK))
65call setBetaInv(betaInvs, betaIncs, alphas, beta, getLogBeta(alpha, beta), signed, infos)
66if (any(infos /= 0)) error stop 'Beta inversion failed.'
67betaInvs
68+0.0000000000000000, +0.28855052598281295E-1, +1.0000000000000000
69getBetaInc(betaInvs, alpha = [.1_RKG, 1._RKG, 10._RKG], beta = 3._RKG, signed = signed)
70+0.0000000000000000, +0.84091340736003467E-1, +1.0000000000000000
71
72

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"
22fileList = glob.glob(pattern)
23for file in fileList:
24
25 df = pd.read_csv(file, delimiter = ",")
26
27 fig = plt.figure(figsize = 1.25 * np.array([6.4, 4.8]), dpi = 200)
28 ax = plt.subplot()
29
30 for i in range(1,len(df.values[0,:]+1)):
31
32 plt.plot( df.values[:, 0]
33 , df.values[:,i]
34 , linewidth = 2
35 )
36
37 plt.xticks(fontsize = fontsize - 2)
38 plt.yticks(fontsize = fontsize - 2)
39 if "abserr" in file:
40 nbit = file.split(".")[1][2:]
41 ax.set_xlabel("X : Regularized Incomplete Beta Function", fontsize = fontsize)
42 ax.set_ylabel("{}-bits Absolute Error: X - BetaInc(BetaInv(X))".format(nbit), fontsize = fontsize)
43 else:
44 ax.set_xlabel("x", fontsize = fontsize)
45 ax.set_ylabel("Regularized Inverse Beta Function", fontsize = fontsize)
46
47 ax.legend ( label
48 , fontsize = fontsize
49 #, loc = "center left"
50 #, bbox_to_anchor = (1, 0.5)
51 )
52
53 plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
54 ax.tick_params(axis = "y", which = "minor")
55 ax.tick_params(axis = "x", which = "minor")
56 plt.tight_layout()
57
58 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.

Author:
Amir Shahmoradi, Oct 16, 2009, 11:14 AM, Michigan

Definition at line 827 of file pm_mathBeta.F90.


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