betran(3f) - [M_datapac:RANDOM] generate beta random numbers
subroutine BETRAN (N,Alpha,Beta,Iseed,X)
INTEGER,intent(in) :: N
REAL(kind=wp),intent(in) :: Alpha
REAL(kind=wp),intent(in) :: Beta
INTEGER,intent(inout) :: Iseed
REAL(kind=wp),intent(out) :: X(:)
BETRAN(3f) generates a random sample of size N from the beta
distribution with shape parameters ALPHA and BETA.
The prototype beta distribution used herein has
mean = ALPHA/(ALPHA+BETA)
and
standard_deviation=sqrt((ALPHA*BETA)/((ALPHA+BETA)**2)*(ALPHA+BETA+1))
This distribution is defined for all X between 0.0 (inclusively)
and 1.0 (inclusively) and has the probability density function
f(x) = (1/constant) * x**(alpha-1) * (1.0-x)**(beta-1)
where the constant = the beta function evaluated at the values ALPHA
and BETA.
N The desired integer number of random numbers to be generated.
ALPHA The value of the first shape parameter.
ALPHA should be greater than or equal to 1.0.
BETA The value of the second shape parameter.
BETA should be greater than or equal to 1.0.
ISEED An integer iseed value. Should be set to a non-negative value to start a new sequence of values. Will be set to -1 on return to indicate the next call should continue the current random sequence walk.
X A random sample of size N from the beta distribution
with shape parameter values ALPHA and BETA.
A vector (of dimension at least N) into which
the generated random sample will be placed.
Sample program:
program demo_betran
use M_datapac, only : betran
implicit none
! call betran(x,y)
end program demo_betran
Results:
The original DATAPAC library was written by James Filliben of the
Statistical Engineering Division, National Institute of Standards
and Technology.
John Urban, 2022.05.31
CC0-1.0
Greenwood, ‘A Fast Generator for Beta-distributed Random Variables’, Compstat 1974, Proceedings in Computational Statistics, Vienna, September, 1974, pages 19-27.
Tocher, The Art of Simulation, 1963, pages 24-27.
Hammersley and Handscomb, Monte Carlo Methods, 1964, pages 36-37.
Johnson and Kotz, Continuous Univariate Distributions –2, 1970, pages 37-56.
Hastings and Peacock, Statistical Distributions–A Handbook For Students and Practitioners, 1975, pages 30-35.
National Bureau of Standards Applied Mathematics Series 55, 1964, page 952.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | N | |||
| real(kind=wp), | intent(in) | :: | Alpha | |||
| real(kind=wp), | intent(in) | :: | Beta | |||
| integer, | intent(inout) | :: | Iseed | |||
| real(kind=wp), | intent(out) | :: | X(:) |
SUBROUTINE BETRAN(N,Alpha,Beta,Iseed,X) INTEGER,intent(in) :: N REAL(kind=wp),intent(in) :: Alpha REAL(kind=wp),intent(in) :: Beta INTEGER,intent(inout) :: Iseed REAL(kind=wp),intent(out) :: X(:) REAL(kind=wp) :: a1, a2, arg, b1, b2, funct, term, u(10), xg, xg01, xg02, xg1, xg2 REAL(kind=wp) :: xn01, xn02, xn(1) INTEGER :: i real(kind=wp),parameter :: athird = 0.33333333_wp real(kind=wp),parameter :: sqrt3= 1.73205081_wp ! ! ***** STILL NEEDS ALGORITHM WORK ****** ! !-----START POINT----------------------------------------------------- ! ! CHECK THE INPUT ARGUMENTS FOR ERRORS ! IF ( N<1 ) THEN WRITE (G_IO,99001) 99001 FORMAT (' ','***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO THE BETRAN SUBROUTINE IS NON-POSITIVE *****') WRITE (G_IO,99002) N 99002 FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',I0,' *****') RETURN ELSEIF ( Alpha<1.0_wp ) THEN WRITE (G_IO,99003) 99003 FORMAT (' ','***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE BETRAN SUBROUTINE IS SMALLER THAN 1.0 *****') WRITE (G_IO,99005) Alpha RETURN ELSEIF ( Beta<1.0_wp ) THEN WRITE (G_IO,99004) 99004 FORMAT (' ','***** FATAL ERROR--THE THIRD INPUT ARGUMENT TO THE BETRAN SUBROUTINE IS SMALLER THAN 1.0 *****') WRITE (G_IO,99005) Beta RETURN ELSE ! ! GENERATE N BETA RANDOM NUMBERS ! BY USING THE FACT THAT IF X1 IS A GAMMA VARIATE WITH PARAMETER ALPHA ! AND IF X2 IS A GAMMA VARIATE WITH PARAMETER BETA, ! THEN THE RATIO X1/(X1+X2) IS A BETA VARIATE WITH PARAMETERS ALPHA AND BETA. ! ! TO GENERATE N GAMMA DISTRIBUTION RANDOM NUMBERS, ! USE GREENWOOD'S REJECTION ALGORITHM-- ! 1) GENERATE A NORMAL RANDOM NUMBER; ! 2) TRANSFORM THE NORMAL VARIATE TO AN APPROXIMATE ! GAMMA VARIATE USING THE WILSON-HILFERTY ! APPROXIMATION (SEE THE JOHNSON AND KOTZ ! REFERENCE, page 176); ! 3) FORM THE REJECTION FUNCTION VALUE, BASED ! ON THE PROBABILITY DENSITY FUNCTION VALUE ! OF THE ACTUAL DISTRIBUTION OF THE PSEUDO-GAMMA ! VARIATE, AND THE PROBABILITY DENSITY FUNCTION VALUE ! OF A TRUE GAMMA VARIATE. ! 4) GENERATE A UNIFORM RANDOM NUMBER; ! 5) IF THE UNIFORM RANDOM NUMBER IS LESS THAN ! THE REJECTION FUNCTION VALUE, THEN ACCEPT ! THE PSEUDO-RANDOM NUMBER AS A GAMMA VARIATE; ! IF THE UNIFORM RANDOM NUMBER IS LARGER THAN ! THE REJECTION FUNCTION VALUE, THEN REJECT ! THE PSEUDO-RANDOM NUMBER AS A GAMMA VARIATE. ! a1 = 1.0_wp/(9.0_wp*Alpha) b1 = SQRT(a1) xn01 = -sqrt3 + b1 xg01 = Alpha*(1.0_wp-a1+b1*xn01)**3 a2 = 1.0_wp/(9.0_wp*Beta) b2 = SQRT(a2) xn02 = -sqrt3 + b2 xg02 = Beta*(1.0_wp-a2+b2*xn02)**3 ! DO i = 1 , N DO ! CALL NORRAN(1,Iseed,xn(1:1)) xg = Alpha*(1.0_wp-a1+b1*xn(1))**3 IF ( xg>=0.0_wp ) THEN term = (xg/xg01)**(Alpha-athird) arg = 0.5_wp*xn(1)*xn(1) - xg - 0.5_wp*xn01*xn01 + xg01 funct = term*EXP(arg) CALL UNIRAN(1,Iseed,u) IF ( u(1)<=funct ) THEN xg1 = xg DO ! CALL NORRAN(1,Iseed,xn(1:1)) xg = Beta*(1.0_wp-a2+b2*xn(1))**3 IF ( xg>=0.0_wp ) THEN term = (xg/xg02)**(Beta-athird) arg = 0.5_wp*xn(1)*xn(1) - xg - 0.5_wp*xn02*xn02 + xg02 funct = term*EXP(arg) CALL UNIRAN(1,Iseed,u) IF ( u(1)<=funct ) THEN xg2 = xg ! X(i) = xg1/(xg1+xg2) GOTO 50 ENDIF ENDIF ENDDO ENDIF ENDIF ENDDO ! 50 ENDDO ENDIF 99005 FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',E15.8,' *****') ! END SUBROUTINE BETRAN