! bayes_beta.f90 20 June 2000
!
program bayes_beta
!
!*******************************************************************************
!
!! BAYES_BETA does a simple demonstration of Bayesian statistics.
!
!
! Modified:
!
! 06 December 1999
!
! Author:
!
! John Burkardt
!
integer ntoss_obs
integer obs_num
real theta1_true
!
write ( *, * ) ' '
write ( *, * ) 'BAYES_BETA'
write ( *, * ) ' Simple Bayesian Statistics demonstrations.'
write ( *, * ) ' '
theta1_true = 0.75
ntoss_obs = 5
obs_num = 4
call test01 ( ntoss_obs, obs_num, theta1_true )
write ( *, * ) ' '
write ( *, * ) 'BAYES_BETA:'
write ( *, * ) ' Normal end of execution.'
stop
end
subroutine test01 ( ntoss_obs, obs_num, theta1_true )
!
!*******************************************************************************
!
!! TEST01 does a simple demonstration of Bayesian statistics.
!
real alpha1_init
real alpha1_post
real alpha1_prior
real alpha2_init
real alpha2_post
real alpha2_prior
integer i
integer iseed
integer nhead_obs
integer nhead_post
integer nhead_prior
integer ntail_obs
integer ntail_post
integer ntail_prior
integer ntoss_obs
integer ntoss_post
integer ntoss_prior
integer obs_i
integer obs_num
integer step_num
real theta1_init
real theta1_post
real theta1_prior
real theta1_true
real theta2_init
real theta2_post
real theta2_prior
real theta2_true
real uniform_01_sample
real x
!
write ( *, * ) ' '
write ( *, * ) 'BAYES_BETA'
write ( *, * ) ' Simple Bayesian Statistics demonstrations.'
write ( *, * ) ' '
write ( *, * ) ' Suppose we''re watching a "system" and trying'
write ( *, * ) ' to analyze its behavior. Each time we observe'
write ( *, * ) ' the system, it flips a coin a certain number of'
write ( *, * ) ' times, and reports the number of heads and'
write ( *, * ) ' tails. We want to estimate THETA1 and THETA2,'
write ( *, * ) ' the probabilities of heads and of tails.'
write ( *, * ) ' '
write ( *, * ) ' We treat the values of THETA1 and THETA2 as'
write ( *, * ) ' random variables themselves, controlled by a'
write ( *, * ) ' Beta probability density function, which has'
write ( *, * ) ' parameters ALPHA1 and ALPHA2. We make an '
write ( *, * ) ' arbitrary or informed guess for initial values'
write ( *, * ) ' of ALPHA1 and ALPHA2. We observe the system,'
write ( *, * ) ' and adjust ALPHA1 and ALPHA2 using Bayes''s'
write ( *, * ) ' Law. We continue until we are satisfied'
write ( *, * ) ' that our estimates seem to have converged.'
!
! Get the initial estimates of the parameters ALPHA1 and ALPHA2.
!
alpha1_init = 0.5
alpha2_init = 0.5
call beta_max ( alpha1_init, alpha2_init, theta1_init, theta2_init )
alpha1_post = alpha1_init
alpha2_post = alpha2_init
nhead_post = 0
ntail_post = 0
ntoss_post = 0
theta1_post = theta1_init
theta2_post = theta2_init
theta1_true = max ( theta1_true, 0.0 )
theta1_true = min ( theta1_true, 1.0 )
theta2_true = 1.0 - theta1_true
call get_seed ( iseed )
!
! Report on run parameters:
!
write ( *, * ) ' '
write ( *, * ) 'Run parameters:'
write ( *, * ) ' '
write ( *, '(a,2g14.6)' ) ' THETA1/THETA2: ', theta1_true, theta2_true
write ( *, * ) ' Number of observations to make = ', obs_num
write ( *, * ) ' Number of tosses per observation = ', ntoss_obs
write ( *, * ) ' The random seed value is ', iseed
!
! Report on initial data.
!
write ( *, * ) ' '
write ( *, * ) ' Initial data:'
write ( *, * ) ' '
write ( *, '(a,2g14.6)' ) ' ALPHA1/ALPHA2 ', alpha1_init, alpha2_init
write ( *, '(a,2g14.6)' ) ' THETA1/THETA2: ', theta1_init, theta2_init
step_num = 20
call beta_plot ( step_num, alpha1_init, alpha2_init )
!
! Make observations of the system.
!
do obs_i = 1, obs_num
alpha1_prior = alpha1_post
alpha2_prior = alpha2_post
nhead_prior = nhead_post
ntail_prior = ntail_post
ntoss_prior = ntoss_post
theta1_prior = theta1_post
theta2_prior = theta2_post
nhead_obs = 0
ntail_obs = 0
!
! Flip the coin NTOSS_OBS times.
!
do i = 1, ntoss_obs
x = uniform_01_sample ( iseed )
if ( x <= theta1_true ) then
nhead_obs = nhead_obs + 1
else
ntail_obs = ntail_obs + 1
end if
end do
!
! Use the observations to adjust our estimates of the system.
!
alpha1_post = alpha1_prior + nhead_obs
alpha2_post = alpha2_prior + ntail_obs
nhead_post = nhead_prior + nhead_obs
ntail_post = ntail_prior + ntail_obs
ntoss_post = ntoss_prior + ntoss_obs
call beta_max ( alpha1_post, alpha2_post, theta1_post, theta2_post )
!
! Print out the data.
!
write ( *, * ) ' '
write ( *, * ) 'BAYES_BETA - After observation ', obs_i
write ( *, * ) ' '
write ( *, * ) ' Prior observation data:'
write ( *, * ) ' '
write ( *, * ) ' Tosses: ', ntoss_prior
write ( *, * ) ' Heads/Tails: ', nhead_prior, ntail_prior
write ( *, '(a,2g14.6)' ) ' ALPHA1/ALPHA2: ', alpha1_prior, &
alpha2_prior
write ( *, '(a,2g14.6)' ) ' THETA1/THETA2: ', theta1_prior, &
theta2_prior
write ( *, * ) ' '
write ( *, * ) ' Observation data:'
write ( *, * ) ' '
write ( *, * ) ' Tosses: ', ntoss_obs
write ( *, * ) ' Heads/Tails: ', nhead_obs, ntail_obs
write ( *, * ) ' '
write ( *, * ) ' Post observation data:'
write ( *, * ) ' '
write ( *, * ) ' Tosses: ', ntoss_post
write ( *, * ) ' Heads/Tails: ', nhead_post, ntail_post
write ( *, '(a,2g14.6)' ) ' ALPHA1/ALPHA2: ', alpha1_post, alpha2_post
write ( *, '(a,2g14.6)' ) ' THETA1/THETA2: ', theta1_post, theta2_post
end do
step_num = 20
call beta_plot ( step_num, alpha1_post, alpha2_post )
return
end
function beta ( x, y )
!
!*******************************************************************************
!
!! BETA returns the value of the Beta function.
!
!
! Formula:
!
! BETA(X,Y) = ( GAMMA(X) * GAMMA(Y) ) / GAMMA(X+Y)
!
! Properties:
!
! BETA(X,Y) = BETA(Y,X).
! BETA(X,Y) = Integral ( 0 <= T <= 1 ) T**(X-1) (1-T)**(Y-1) dT.
!
! Restrictions:
!
! Both X and Y must be greater than 0.
!
! Modified:
!
! 16 June 1999
!
! Author:
!
! John Burkardt
! Parameters:
!
! Input, real X, Y, the two parameters that define the Beta function.
!
! Output, real BETA, the value of the Beta function.
!
real beta
real gamma_log
real x
real y
!
if ( x <= 0.0 .or. y <= 0.0 ) then
write ( *, * ) ' '
write ( *, * ) 'BETA - Fatal error!'
write ( *, * ) ' Both X and Y must be greater than 0.'
stop
end if
beta = exp ( gamma_log ( x ) + gamma_log ( y ) - gamma_log ( x + y ) )
return
end
subroutine beta_max ( alpha1, alpha2, t1, t2 )
!
!*******************************************************************************
!
!! BETA_MAX returns the most likely values of T1 and T2 for a Beta PDF.
!
!
! Modified:
!
! 31 March 1999
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, real ALPHA1, ALPHA2, the parameter values for the Beta
! probability density function.
!
! Output, real T1, T2, the most likely values of T1 and T2, which
! are controlled by the Beta probability density function.
!
real alpha1
real alpha2
real t1
real t2
!
t1 = alpha1 / ( alpha1 + alpha2 )
t2 = alpha2 / ( alpha1 + alpha2 )
return
end
function beta_pdf ( t, x, y )
!
!*******************************************************************************
!
!! BETA_PDF returns the value of the Beta probability density function.
!
!
! Formula:
!
! BETA_PDF(T,X,Y) = T**(X-1) * (1-T)**(Y-1) / BETA(X,Y).
!
! Properties:
!
! BETA_PDF(T,X,Y) = BETA_PDF(T,Y,X).
! BETA_PDF(0,X,Y) = 0.
! BETA_PDF(1,X,Y) = 0.
! Integral ( 0 <= T <= 1 ) BETA_PDF(T,X,Y) dT = 1.
! Expected value of T is X / ( X + Y ).
!
! Restrictions:
!
! T must be in [0,1], X and Y must be greater than 0.
!
! Modified:
!
! 16 June 1999
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, real T, the value at which the Beta probability density
! function is to be evaluated.
!
! Input, real X, Y, the two parameters that define the Beta function.
! X and Y must be greater than 0.
!
! Output, real BETA_PDF, the value of the Beta function.
!
real beta_log
real beta_pdf
real gamma_log
real t
real x
real y
!
if ( x <= 0.0 .or. y <= 0.0 ) then
write ( *, * ) ' '
write ( *, * ) 'BETA_PDF - Fatal error!'
write ( *, * ) ' Both X and Y must be greater than 0.'
stop
end if
if ( t < 0.0 .or. t > 1.0 ) then
write ( *, * ) ' '
write ( *, * ) 'BETA_PDF - Fatal error!'
write ( *, * ) ' 0 <= T <= 1 is required.'
write ( *, * ) ' Your value of T is ', t
stop
end if
!
! If you use the straightforward formula for the Beta PDF, you
! will get intermediate overflow with values of X and Y as small
! as 200, although the values of BETA_PDF are reasonable.
!
if ( t == 0.0 .or. t == 1.0 ) then
beta_pdf = 0.0
else
beta_log = ( x - 1.0 ) * log ( t ) + ( y - 1.0 ) * log ( 1.0 - t ) &
+ gamma_log ( x + y ) - gamma_log ( x ) - gamma_log ( y )
beta_pdf = exp ( beta_log )
end if
return
end
subroutine beta_plot ( step_num, alpha1, alpha2 )
!
!*******************************************************************************
!
!! BETA_PLOT "plots" the Beta distribution for given parameters.
!
!
! Modified:
!
! 06 December 1999
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer STEP_NUM, the number of steps to take from 0 to 1
! in the plotting.
!
! Input, real ALPHA1, ALPHA2, the parameters of the distribution.
!
real alpha1
real alpha2
real beta_pdf
real pdf
integer step_i
integer step_num
real t
!
write ( *, * ) ' '
write ( *, * ) 'Plot of Beta distribution,'
write ( *, * ) ' '
write ( *, * ) ' ALPHA1 = ', alpha1
write ( *, * ) ' ALPHA2 = ', alpha2
write ( *, * ) ' '
write ( *, * ) ' T BETA(T,ALPHA1,ALPHA2)'
write ( *, * ) ' '
do step_i = 0, step_num
t = real ( step_i ) / real ( step_num )
pdf = beta_pdf ( t, alpha1, alpha2 )
write ( *, '(2g14.6)' ) t, pdf
end do
return
end
function gamma_log ( x )
!
!*******************************************************************************
!
!! GAMMA_LOG calculates the natural logarithm of GAMMA ( X ) for positive X.
!
!
! Discussion:
!
! Computation is based on an algorithm outlined in references 1 and 2.
! The program uses rational functions that theoretically approximate
! LOG(GAMMA(X)) to at least 18 significant decimal digits. The
! approximation for X > 12 is from reference 3, while approximations
! for X < 12.0 are similar to those in reference 1, but are unpublished.
! The accuracy achieved depends on the arithmetic system, the compiler,
! intrinsic functions, and proper selection of the machine-dependent
! constants.
!
! Modified:
!
! 16 June 1999
!
! Authors:
!
! W. J. Cody and L. Stoltz
! Argonne National Laboratory
!
! References:
!
! # 1)
! W. J. Cody and K. E. Hillstrom,
! Chebyshev Approximations for the Natural Logarithm of the Gamma Function,
! Math. Comp.
! Volume 21, 1967, pages 198-203.
!
! # 2)
! K. E. Hillstrom,
! ANL/AMD Program ANLC366S, DGAMMA/DLGAMA,
! May 1969.
!
! # 3)
! Hart, Et. Al.,
! Computer Approximations,
! Wiley and sons, New York, 1968.
!
! Parameters:
!
! Input, real X, the argument of the Gamma function. X must be positive.
!
! Output, real GAMMA_LOG, the logarithm of the Gamma function of X.
! If X <= 0.0, or if overflow would occur, the program returns the
! value XINF, the largest representable floating point number.
!
!*******************************************************************************
!
! Explanation of machine-dependent constants
!
! BETA - radix for the floating-point representation.
!
! MAXEXP - the smallest positive power of BETA that overflows.
!
! XBIG - largest argument for which LN(GAMMA(X)) is representable
! in the machine, i.e., the solution to the equation
! LN(GAMMA(XBIG)) = BETA**MAXEXP.
!
! XINF - largest machine representable floating-point number;
! approximately BETA**MAXEXP.
!
! EPS - The smallest positive floating-point number such that
! 1.0+EPS .GT. 1.0
!
! FRTBIG - Rough estimate of the fourth root of XBIG
!
!
! Approximate values for some important machines are:
!
! BETA MAXEXP XBIG
!
! CRAY-1 (S.P.) 2 8191 9.62E+2461
! Cyber 180/855
! under NOS (S.P.) 2 1070 1.72E+319
! IEEE (IBM/XT,
! SUN, etc.) (S.P.) 2 128 4.08E+36
! IEEE (IBM/XT,
! SUN, etc.) (D.P.) 2 1024 2.55D+305
! IBM 3033 (D.P.) 16 63 4.29D+73
! VAX D-Format (D.P.) 2 127 2.05D+36
! VAX G-Format (D.P.) 2 1023 1.28D+305
!
!
! XINF EPS FRTBIG
!
! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 3.13E+615
! Cyber 180/855
! under NOS (S.P.) 1.26E+322 3.55E-15 6.44E+79
! IEEE (IBM/XT,
! SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.42E+9
! IEEE (IBM/XT,
! SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.25D+76
! IBM 3033 (D.P.) 7.23D+75 2.22D-16 2.56D+18
! VAX D-Format (D.P.) 1.70D+38 1.39D-17 1.20D+9
! VAX G-Format (D.P.) 8.98D+307 1.11D-16 1.89D+76
!
real, parameter :: d1 = - 5.772156649015328605195174e-1
real, parameter :: d2 = 4.227843350984671393993777e-1
real, parameter :: d4 = 1.791759469228055000094023e0
real, parameter :: EPS = 1.19e-7
real, parameter :: FRTBIG = 1.42e9
real, parameter :: PNT68 = 0.6796875
real, parameter :: SQRTPI = 0.9189385332046727417803297
real, parameter :: XBIG = 4.08e36
real, parameter :: XINF = 3.401e38
!
real, parameter, dimension ( 7 ) :: c = (/ &
-1.910444077728e-03, &
8.4171387781295e-04, &
-5.952379913043012e-04, &
7.93650793500350248e-04, &
-2.777777777777681622553e-03, &
8.333333333333333331554247e-02, &
5.7083835261e-03 /)
real corr
integer i
real gamma_log
real, parameter, dimension ( 8 ) :: p1 = (/ &
4.945235359296727046734888e0, &
2.018112620856775083915565e2, &
2.290838373831346393026739e3, &
1.131967205903380828685045e4, &
2.855724635671635335736389e4, &
3.848496228443793359990269e4, &
2.637748787624195437963534e4, &
7.225813979700288197698961e3 /)
real, parameter, dimension ( 8 ) :: p2 = (/ &
4.974607845568932035012064e0, &
5.424138599891070494101986e2, &
1.550693864978364947665077e4, &
1.847932904445632425417223e5, &
1.088204769468828767498470e6, &
3.338152967987029735917223e6, &
5.106661678927352456275255e6, &
3.074109054850539556250927e6 /)
real, parameter, dimension ( 8 ) :: p4 = (/ &
1.474502166059939948905062e4, &
2.426813369486704502836312e6, &
1.214755574045093227939592e8, &
2.663432449630976949898078e9, &
2.940378956634553899906876e10, &
1.702665737765398868392998e11, &
4.926125793377430887588120e11, &
5.606251856223951465078242e11 /)
real, parameter, dimension ( 8 ) :: q1 = (/ &
6.748212550303777196073036e1, &
1.113332393857199323513008e3, &
7.738757056935398733233834e3, &
2.763987074403340708898585e4, &
5.499310206226157329794414e4, &
6.161122180066002127833352e4, &
3.635127591501940507276287e4, &
8.785536302431013170870835e3 /)
real, parameter, dimension ( 8 ) :: q2 = (/ &
1.830328399370592604055942e2, &
7.765049321445005871323047e3, &
1.331903827966074194402448e5, &
1.136705821321969608938755e6, &
5.267964117437946917577538e6, &
1.346701454311101692290052e7, &
1.782736530353274213975932e7, &
9.533095591844353613395747e6 /)
real, parameter, dimension ( 8 ) :: q4 = (/ &
2.690530175870899333379843e3, &
6.393885654300092398984238e5, &
4.135599930241388052042842e7, &
1.120872109616147941376570e9, &
1.488613728678813811542398e10, &
1.016803586272438228077304e11, &
3.417476345507377132798597e11, &
4.463158187419713286462081e11 /)
real res
real x
real xden
real xm1
real xm2
real xm4
real xnum
real xsq
!
! Return immediately if the argument is out of range.
!
if ( x <= 0.0 .or. x > XBIG ) then
gamma_log = XINF
return
end if
if ( x <= EPS ) then
res = - log ( x )
else if ( x <= 1.5 ) then
if ( x < PNT68 ) then
corr = - log ( x )
xm1 = x
else
corr = 0.0
xm1 = ( x - 0.5 ) - 0.5
end if
if ( x <= 0.5 .or. x >= PNT68 ) then
xden = 1.0
xnum = 0.0
do i = 1, 8
xnum = xnum * xm1 + p1(i)
xden = xden * xm1 + q1(i)
end do
res = corr + ( xm1 * ( d1 + xm1 * ( xnum / xden ) ) )
else
xm2 = ( x - 0.5 ) - 0.5
xden = 1.0
xnum = 0.0
do i = 1, 8
xnum = xnum * xm2 + p2(i)
xden = xden * xm2 + q2(i)
end do
res = corr + xm2 * ( d2 + xm2 * ( xnum / xden ) )
end if
else if ( x <= 4.0 ) then
xm2 = x - 2.0
xden = 1.0
xnum = 0.0
do i = 1, 8
xnum = xnum * xm2 + p2(i)
xden = xden * xm2 + q2(i)
end do
res = xm2 * ( d2 + xm2 * ( xnum / xden ) )
else if ( x <= 12.0 ) then
xm4 = x - 4.0
xden = - 1.0
xnum = 0.0
do i = 1, 8
xnum = xnum * xm4 + p4(i)
xden = xden * xm4 + q4(i)
end do
res = d4 + xm4 * ( xnum / xden )
else
res = 0.0
if ( x <= FRTBIG ) then
res = c(7)
xsq = x * x
do i = 1, 6
res = res / xsq + c(i)
end do
end if
res = res / x
corr = log ( x )
res = res + SQRTPI - 0.5 * corr
res = res + x * ( corr - 1.0 )
end if
gamma_log = res
return
end
subroutine get_seed ( iseed )
!
!*******************************************************************************
!
!! GET_SEED returns a seed for the random number generator.
!
!
! Discussion:
!
! The seed depends on the current time, and ought to be (slightly)
! different every millisecond. Once the seed is obtained, a random
! number generator should be called a few times to further process
! the seed.
!
! If a FORTRAN77 compiler is used, then a different set of routines
! must be invoked in order to get the date and time.
!
! Modified:
!
! 23 March 2000
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Output, integer ISEED, a pseudorandom seed value.
!
integer, parameter :: I_MAX = 2147483647
!
integer iseed
double precision temp
integer values(8)
!
! FORTRAN90 declarations
!
character ( len = 10 ) time90
character ( len = 8 ) today90
character ( len = 5 ) zone
!
! FORTRAN77 declarations
!
! integer day
! integer hour
! integer minute
! integer month
! integer second
! character*8 time77
! character*9 today77
! integer year
!
! FORTRAN90
!
call date_and_time ( today90, time90, zone, values )
!
! FORTRAN77
!
! call date ( today77 )
! call s_to_ymd ( today77, 'DD-NNN-YY', year, month, day )
! call time ( time77 )
! call s_to_hms ( time77, 'hh:mm:ss', hour, minute, second )
!
! if ( year == 99 ) then
! year = 1900 + year
! else
! year = 2000 + year
! end if
!
! values(1) = year
! values(2) = month
! values(3) = day
! values(4) = 0
! values(5) = hour
! values(6) = minute
! values(7) = second
! values(8) = 0
!
temp = 0.0
temp = temp + dble ( values(2) - 1 ) / 11.0
temp = temp + dble ( values(3) - 1 ) / 30.0
temp = temp + dble ( values(5) ) / 23.0
temp = temp + dble ( values(6) ) / 59.0
temp = temp + dble ( values(7) ) / 59.0
temp = temp + dble ( values(8) ) / 999.0
temp = temp / 6.0
if ( temp <= 0.0 ) then
temp = 1.0 / 3.0
else if ( temp >= 1.0 ) then
temp = 2.0 / 3.0
end if
iseed = int ( dble ( I_MAX ) * temp )
!
! Never use a seed of 0 or I_MAX.
!
if ( iseed == 0 ) then
iseed = 1
end if
if ( iseed == I_MAX ) then
iseed = I_MAX - 1
end if
return
end
function uniform_01_sample ( iseed )
!
!*******************************************************************************
!
!! UNIFORM_01_SAMPLE is a portable random number generator.
!
!
! Formula:
!
! ISEED = ISEED * (7**5) mod (2**31 - 1)
! UNIFORM_01_SAMPLE = ISEED * / ( 2**31 - 1 )
!
! Parameters:
!
! Input/output, integer ISEED, the integer "seed" used to generate
! the output random number, and updated in preparation for the
! next one. ISEED should not be zero.
!
! Output, real UNIFORM_01_SAMPLE, a random value between 0 and 1.
!
integer, parameter :: ia = 16807
integer, parameter :: ib15 = 32768
integer, parameter :: ib16 = 65536
integer, parameter :: ip = 2147483647
!
integer iprhi
integer iseed
integer ixhi
integer k
integer leftlo
integer loxa
real uniform_01_sample
!
! Don't let ISEED be 0.
!
if ( iseed == 0 ) then
iseed = ip
end if
!
! Get the 15 high order bits of ISEED.
!
ixhi = iseed / ib16
!
! Get the 16 low bits of ISEED and form the low product.
!
loxa = ( iseed - ixhi * ib16 ) * ia
!
! Get the 15 high order bits of the low product.
!
leftlo = loxa / ib16
!
! Form the 31 highest bits of the full product.
!
iprhi = ixhi * ia + leftlo
!
! Get overflow past the 31st bit of full product.
!
k = iprhi / ib15
!
! Assemble all the parts and presubtract IP. The parentheses are
! essential.
!
iseed = ( ( ( loxa - leftlo * ib16 ) - ip ) + ( iprhi - k * ib15 ) * ib16 ) &
+ k
!
! Add IP back in if necessary.
!
if ( iseed < 0 ) then
iseed = iseed + ip
end if
!
! Multiply by 1 / (2**31-1).
!
uniform_01_sample = real ( iseed ) * 4.656612875e-10
return
end