The Gamma distribution

1_Introduction

In what follows I will present a set of codes that I have written in Octave and in Fortran for the study of random variables that follow a gamma distribution. In another post, we will discuss the beta distribution.

2_The gamma distribution

A random variable is said to have a gamma distribution with parameters (α, λ) if its density is given by

with α, λ>0 and where Γ(α) is the gamma function, which is quite a pain in the lower back, as you will see. This random variable can be used to describe the probability that the duration of life of an entity (whether it is a mechanical component or a bacterium) depends on time. The gamma distribution allows to describe the following situations:

1. for α>1 the probability of failure of a mechanical component (or death of an organism) increases as the time passes by (increasing hazard rate);
2. for α=1 the probability of failure is constant (constant hazard rate);
3. for α<1 the failure is less likely as time passes by (decreasing hazard rate).

The hazard rate is the limit below:

where F_X(t) is the distribution function of the gamma distribution. I have not been able to find an analytic expression for the hazard rate of the gamma distribution, but I have been able to obtain it with the computer-aided calculation (see below). But in order to calculate F_X(t), we need to calculate the gamma function…

3_The gamma function

Γ(α) (the gamma function) is an integral function that cannot be calculated analytically:

It has to be calculated with numerical integration, but a problem is that Γ(α) is an indefinite integral. One way to face this problem is to stop the numerical integration when the addend is below a certain cut off (let’s say 0.0001). It is worth noting that an important feature of this function is that if you know its value in an interval such as [x, x+1], then you can derive Γ(α) for any other possible value of α. This is due to the recursive formula: Γ(α+1) = αΓ(α). That said, I have calculated Γ(α) in [1, 2] and then I have used these values to define Γ(α) in [0, 4] with the .exe file (written in Fortran) available here. Download it and put it in a folder. When you run it, it gives a .txt file with the values of Γ(α) in [0, 4] and two .mat files that can be used to build the plot of Γ(α) (see figure 1). Values of Γ(α) in [0, 2] are tabulated below, while those of you who like the very long list of numbers can find a huge table with values of the gamma function in [0, 4] at the end of this article.

You can find below the script of main_f_gamma.exe and its module.

!Data: 1/03/2018
PROGRAM main_f_gamma

USE DISLIN !libreria grafica
USE mod_f_gamma !modulo con procedure

IMPLICIT NONE

REAL(KIND=num), DIMENSION (1:400):: f_gamma !valori della funzione gamma fra 1 e 4
REAL(KIND=num), DIMENSION (1:400):: alpha !valori della ascissa fra 1 e 4
CHARACTER(len=10):: chiusura !serve a chiudere il programma
INTEGER:: i,j !indici di ciclo

CALL sub_f_gamma_1_2 (f_gamma,alpha) !calcola la funzione gamma fra 1 e 2

!calcolo della funzione gamma e della ascissa fra 2 e 4

ciclo: DO i=201,400,1
alpha(i) = 1. + alpha(i-100)
f_gamma(i) = alpha(i-100)*f_gamma(i-100)
END DO ciclo

!calcolo della funzione gamma e della ascissa fra 0 e 1

ciclo: DO i=2,100,1
alpha(i) = alpha(100+i)-1.
f_gamma(i) = f_gamma(100+i)/alpha(i)
END DO ciclo

!scrivo la funzione gamma tra 0 e 4 su un file di testo

OPEN(UNIT=20, FILE=’f_gamma.txt’, STATUS=’REPLACE’, ACTION=’WRITE’)
WRITE(20,*) “Gamma function”
WRITE(20,*) ” “
WRITE(20,*) ” alpha”,”|”,”Gamma(alpha)”
WRITE(20,*) “_________________________________”

ciclo_write: DO j=2,400,1
WRITE(20,100) alpha(j),”|”,f_gamma(j)
END DO ciclo_write
100 FORMAT (1X,1F8.4,A,1F8.4)

!salva alpha e f_gamma su un file

OPEN(UNIT=20, FILE=’f_gamma.mat’, STATUS=’REPLACE’, ACTION=’WRITE’)
WRITE(20,200) f_gamma
OPEN(UNIT=20, FILE=’alpha.mat’, STATUS=’REPLACE’, ACTION=’WRITE’)
WRITE(20,200) alpha
200 FORMAT (1X,400(F8.4,/))

WRITE(*,*) “To close this program press any key, then press RETURN.”
READ(*,*) chiusura

END PROGRAM main_f_gamma

!Data: 1/03/2018

MODULE mod_f_gamma

!sezione dichiarativa

IMPLICIT NONE

!costanti

REAL,PARAMETER:: e = 2.71828182 !numero di Nepero
REAL,PARAMETER:: d_alpha = 0.01 !incremento delle ascisse della funzione gamma
INTEGER,PARAMETER:: num = SELECTED_REAL_KIND (p=13) !aumento il numero di cifre significative

CONTAINS

!————————————————————————————-

SUBROUTINE sub_f_gamma_1_2 (f_gamma, alpha)

!Esegue l’integrale necessario al calcolo della funzione gamma
!per valori compresi fra 1 e 2 e registra questi valori in un array.
!L’integrazione viene eseguita con il metodo di Cavalieri-Simpson

IMPLICIT NONE

!dichiaro gli argomenti fittizi

REAL(KIND=num),INTENT(OUT),DIMENSION(1:400):: f_gamma !valori della funz. gamma in [1,2]
REAL(KIND=num),INTENT(OUT),DIMENSION(1:400):: alpha !contiene l’ascissa della funzione gamma

!dichiaro le variabili locali

REAL(KIND=num),DIMENSION(1:1000000):: funzione !è la funzione integranda
REAL(KIND=num),DIMENSION(1:1000000):: primitiva !è l’integrale illimitato
REAL(KIND=num),PARAMETER:: d_x = 0.001 !intervallo di integrazione
INTEGER:: i,j !indice del ciclo
INTEGER:: index = 1 !contiene l’estremo superiore della integrazione numerica

!sezione esecutiva

!assegno i valori di funzione e ascissa, eseguo l’integrazione numerica, assegno i valori di f_gamma_1_2

alpha(101) = 1. !assegna il primo valore di alpha
alpha(102) = alpha(101)+d_alpha !assegna il primo valore di alpha
f_gamma(101) = 1. !assegna il primo valore di gamma

ciclo_esterno: DO j=102,200,1 !assegna i valori a f_gamma_1_2

!assegna i valori a funzione

ciclo_funzione: DO i=1,1000000,1

funzione (i) = ( (d_x*REAL(i))**(alpha(j)-1.) )*( e**((-1)*d_x*REAL(i)) )
index = i
IF (funzione(i)<=0.0001) EXIT ciclo_funzione

END DO ciclo_funzione

!assegna i valori a primitiva

primitiva (1) = 0.
primitiva (1+2) = primitiva (1) + ( ( funzione(1+2) + funzione(1+1)*4. + funzione(1) )*(d_x/3.) )
primitiva (2) = primitiva (3)/2.

ciclo_integrazione: DO i=2,index-2,1 !integrazione numerica

primitiva (i+2) = primitiva (i) + ( ( funzione(i+2) + funzione(i+1)*4. + funzione(i) )*(d_x/3.) )

END DO ciclo_integrazione

f_gamma(j) = primitiva(index) !assegna il valore j-mo a f_gamma_1_2

alpha (j+1) = alpha (j) + d_alpha

END DO ciclo_esterno

WRITE(*,*) alpha
WRITE(*,*) f_gamma

RETURN

END SUBROUTINE sub_f_gamma_1_2

!————————————————————————————-

END MODULE mod_f_gamma

4_The distribution function

The distribution function of the gamma density cannot be calculated analytically so, again, I had to make good use of my poor coding skills and I wrote a code in Octave (here, it uses the two .mat files generated by the code in Fortran mentioned above!) and one in Fortran (available here). In figure 2 (code in Octave) you can find the density and the distribution function for different values of α.

The code in Fortran, on the other hand, asks the user for the values of α and λ and gives a plot for the density and one for the distribution function. Moreover, it gives a .txt file with the values of the distribution function. Note that you can only place values for α that are below 25, otherwise you generate an overflow condition. I am not sure why you have overflow condition with a value of α that is below 34, but this is another story. Nevertheless, for those who are interested, Γ(n) = (n-1)! and the biggest factorial that can be handled by a computer in single precision mode (where 8 bits are used for the exponent) is 33!. So, we have a little mystery here… The reason why I have stuck with single precision is that DISLIN – the graphics library – seems to be unable to handle double precision numbers.  If you are interested in the gamma function, the density and the distribution function for α = 20 and λ=2, then the code in Fortran gives for Γ(20) a value of 1.216451E+17 and the two plots in figure 3.

In what follows you find the scripts of the Fortran .exe mentioned in this paragraph.

PROGRAM main_legge_gamma
USE DISLIN !libreria grafica
USE mod_legge_gamma !modulo con procedure

IMPLICIT NONE

REAL(KIND=num):: alph, lamb !valore di alpha e di lambda
REAL(KIND=num):: f_g !il valore della funzione gamma
REAL(KIND=num):: estremo !estremo superiore di x
REAL(KIND=num):: delta !intervallo di integrazione numerica
REAL(KIND=num), DIMENSION (1:500000):: f_gamma !valori della funzione gamma fra 0 e 5000
REAL(KIND=num), DIMENSION (1:500000):: alpha !valori della ascissa fra 0 e 5000
REAL(KIND=num), DIMENSION (1:1000):: d_g !contiene la densita’
REAL(KIND=num), DIMENSION (1:1000):: f_r !funzione di ripartizione
REAL(KIND=num), DIMENSION (1:1000):: x !ascissa
CHARACTER(len=10):: chiusura !serve a chiudere il programma
INTEGER:: i,j,k !indici di ciclo
INTEGER:: ind !il valore di alpha in cui calcolo la funz. gamma

!calcolo della funzione gamma fra 1 e 2

CALL sub_f_gamma_1_2 (f_gamma,alpha) !calcola la funzione gamma fra 1 e 2

WRITE(*,*) “This program calculates the gamma distribution, for alpha and lambda specified”
WRITE(*,*) “by the user. It plots the diagram of the density and distribution function, and”
WRITE(*,*) “writes a file .txt with some values of the distribution function.”
WRITE(*,*) ” “
WRITE(*,*) “Please, type the value of alpha (a positive number between 0.01 and 25.0).”
WRITE(*,*) “Use dot as decimal separator.”
READ(*,*) alph
WRITE(*,*) ” “
WRITE(*,*) “Please, type the value of lambda (a positive number).”
WRITE(*,*) “Use dot as decimal separator.”
READ(*,*) lamb

!calcolo della funzione gamma e della ascissa fra 2 e 25

ciclo_2_5000: DO i=201,2500,1
alpha(i) = 1. + alpha(i-100)
f_gamma(i) = alpha(i-100)*f_gamma(i-100)
END DO ciclo_2_5000

!calcolo della funzione gamma e della ascissa fra 0 e 1

ciclo_0_1: DO i=2,100,1
alpha(i) = alpha(100+i)-1.
f_gamma(i) = f_gamma(100+i)/alpha(i)
END DO ciclo_0_1

!cerco il valore di alpha specificato dall’utente

ciclo_alph: DO i=1,2500,1
ind = i
IF (ABS(alpha(i)-alph)<=0.01) EXIT ciclo_alph
END DO ciclo_alph

!assegno il valore alla funzione gamma

f_g = f_gamma(ind)

WRITE(*,*) “The value of the gamma function in”, alph, “is”, f_g

!calcolo la densita’ gamma

IF (alph <= 1) THEN
estremo = 2.
delta = estremo/1000.
ciclo_densita_1: DO k=1,1000,1
x(k)= delta*k
d_g(k) = ( ( lamb**alph )/( f_g ) )*( x(k)**( alph – 1. ) )*( e**( (-1)*lamb*x(k) ) )
END DO ciclo_densita_1
ELSE
estremo = ( (alph-1.)/lamb )*3.
delta = estremo/1000.
ciclo_densita_1: DO k=1,1000,1
x(k)= delta*(k-1.)
d_g(k) = ( ( lamb**alph )/( f_g ) )*( x(k)**( alph – 1. ) )*( e**( (-1)*lamb*x(k) ) )
END DO ciclo_densita_1
END IF

!calcolo la funzione di ripartizione

f_r(1) = 0.
f_r(3) = f_r(1) + delta*( d_g(1) + ( 4*d_g(2) ) + d_g(3) )/3.
f_r(2) = f_r(3)*0.5
ciclo_f_r: DO k=2,998,1
f_r(k+2) = f_r(k) + delta*( d_g(k) + ( 4*d_g(k+1) ) + d_g(k+2) )/3.;
END DO ciclo_f_r

!scrivo la funzione di ripartizione su un file di testo

OPEN(UNIT=20, FILE=’f_ripartizione.txt’, STATUS=’REPLACE’, ACTION=’WRITE’)

WRITE(20,*) “The gamma distribution for alpha =”, alph, “and lambda =”, lamb
WRITE(20,*) ” “
WRITE(20,*) “The value of the gamma function in”, alph, “is”, f_g
WRITE(20,*) ” “
WRITE(20,*) “distribution function”
WRITE(20,*) ” “
WRITE(20,*) ” x “,”|”,”distribution function”
WRITE(20,*) “_________________________________”

ciclo_write: DO j=1,1000,1
WRITE(20,100) x(j),”|”,f_r(j)
END DO ciclo_write
100 FORMAT (1X,1F8.4,A,1F8.4)

!disegno il grafico della densita’

CALL diagramma (x, d_g, ‘Density ‘)

!disegna il grafico della funzione di ripartizione

CALL diagramma (x, f_r, ‘distribution function’)

WRITE(*,*) “To close this program press any key, then press RETURN.”
READ(*,*) chiusura

END PROGRAM main_legge_gamma

!Data: 1/03/2018

MODULE mod_legge_gamma

!sezione dichiarativa

IMPLICIT NONE

!costanti

REAL,PARAMETER:: e = 2.71828182 !numero di Nepero
REAL,PARAMETER:: d_alpha = 0.01 !incremento delle ascisse della funzione gamma
INTEGER,PARAMETER:: num = 1 !singola precisione

CONTAINS

!————————————————————————————-

SUBROUTINE sub_f_gamma_1_2 (f_gamma, alpha)

!Esegue l’integrale necessario al calcolo della funzione gamma
!per valori compresi fra 1 e 2 e registra questi valori in un array.
!L’integrazione viene eseguita con il metodo di Cavalieri-Simpson

IMPLICIT NONE

!dichiaro gli argomenti fittizi

REAL(KIND=num),INTENT(OUT),DIMENSION(:):: f_gamma !valori della funz. gamma in [1,2]
REAL(KIND=num),INTENT(OUT),DIMENSION(:):: alpha !contiene l’ascissa della funzione gamma

!dichiaro le variabili locali

REAL(KIND=num),DIMENSION(1:1000000):: funzione !è la funzione integranda
REAL(KIND=num),DIMENSION(1:1000000):: primitiva !è l’integrale illimitato
REAL(KIND=num),PARAMETER:: d_x = 0.001 !intervallo di integrazione
INTEGER:: i,j !indice del ciclo
INTEGER:: index = 1 !contiene l’estremo superiore della integrazione numerica

!sezione esecutiva

!assegno i valori di funzione e ascissa, eseguo l’integrazione numerica, assegno i valori di f_gamma_1_2

alpha(101) = 1. !assegna il primo valore di alpha
alpha(102) = alpha(101)+d_alpha !assegna il primo valore di alpha
f_gamma(101) = 1. !assegna il primo valore di gamma

ciclo_esterno: DO j=102,200,1 !assegna i valori a f_gamma_1_2

!assegna i valori a funzione

ciclo_funzione: DO i=1,1000000,1

funzione (i) = ( (d_x*REAL(i))**(alpha(j)-1.) )*( e**((-1)*d_x*REAL(i)) )
index = i
IF (funzione(i)<=0.0001) EXIT ciclo_funzione

END DO ciclo_funzione

!assegna i valori a primitiva

primitiva (1) = 0.
primitiva (1+2) = primitiva (1) + ((funzione(1+2)+funzione(1+1)*4.+funzione(1))*(d_x/3.))
primitiva (2) = primitiva (3)/2.

ciclo_integrazione: DO i=2,index-2,1 !integrazione numerica

primitiva(i+2)= primitiva (i)+((funzione(i+2)+funzione(i+1)*4.+funzione(i) )*(d_x/3.))

END DO ciclo_integrazione

f_gamma(j) = primitiva(index) !assegna il valore j-mo a f_gamma_1_2

alpha (j+1) = alpha (j) + d_alpha

END DO ciclo_esterno

RETURN

END SUBROUTINE sub_f_gamma_1_2

!————————————————————————————-

SUBROUTINE diagramma (x, funzione,label_y)

!disegno il grafico di funzione

!dichiaro gli argomenti fittizi

REAL(KIND=num),INTENT(IN),DIMENSION(1:1000):: x !ascissa
REAL(KIND=num),INTENT(IN),DIMENSION(1:1000):: funzione !funzione
CHARACTER (len=20), INTENT(IN):: label_y !nome delle ordinate

!dichiaro le variabili locali

REAL(KIND=num),DIMENSION(1:1000):: max_x !estremo ascissa
REAL(KIND=num),DIMENSION(1:1000):: max_funz !estremo funzione

!sezione esecutiva

max_x = MAXVAL(x) !valore massimo dlla ascissa
max_funz = MAXVAL(funzione) !valore massimo della ordinata

CALL METAFL (‘PDF’) !scelgo il formato dell’output
CALL SETPAG (‘DA4P’) !formato della pagina
CALL SCRMOD (‘revers’) !scritta nera su fondo bianco
CALL DISINI !inizializza DISLIN impostando parametri di default

CALL PAGERA !chiedo una linea di contorno per il piano xy
CALL DUPLX !chiedo un font a doppio spessore

CALL AXSPOS (300,1400) !coordinate dell’angolo in basso a sinistra
CALL AXSLEN (1500,1000) !la lunghezza, in pixel, dei due assi

CALL LABDIG (2,’X’) !numero di cifre decimali negli assi
CALL LABDIG (3,’Y’) !numero di cifre decimali negli assi

CALL NAME (‘X’,’x’) !assegno un nome all’asse delle x
CALL NAME (label_y,’y’) !assegno un nome all’asse delle y

CALL GRAF (0.,max_x,0.,max_x/5.,0.,max_funz+max_funz/5.,0.,max_funz/5.) !ambiti degli assi
CALL TEXMOD (‘ON’) !modalita’ TeX
CALL GRID (2,2) !impone una griglia
CALL CURVE (x,funzione,1000)

CALL DISFIN !termina DISLIN

END SUBROUTINE diagramma

!————————————————————————————-

END MODULE mod_legge_gamma

5_The hazard rate

As mentioned in paragraph 2, the hazard rate of the gamma distribution cannot be easily studied. I forgot to mention that it is relatively easy to prove that it is constant for α=1. But what happens for α>1 and for α<1? In order to answer, I have written a code (download) that plots it and it seems to be possible to say that the hazard rate increases for α>1 and decreases for α<1. This means that the gamma distribution is suitable for modeling each of the three kinds of patterns for the probability of failure as the time passes by. In figure 4 you find an example of the hazard rate for α<1 (left) and an example for α>1 (right).

This is the script of the Fortran exe program mentioned in this paragraph:

!Data: 7/03/2018

PROGRAM main_legge_gamma_avaria

USE DISLIN !libreria grafica
USE mod_legge_gamma_avaria !modulo con procedure

IMPLICIT NONE

REAL(KIND=num):: alph, lamb !valore di alpha e di lambda
REAL(KIND=num):: f_g !il valore della funzione gamma
REAL(KIND=num):: estremo !estremo superiore di x
REAL(KIND=num):: delta !intervallo di integrazione numerica
REAL(KIND=num), DIMENSION (1:500000):: f_gamma !valori della funzione gamma fra 0 e 5000
REAL(KIND=num), DIMENSION (1:500000):: alpha !valori della ascissa fra 0 e 5000
REAL(KIND=num), DIMENSION (1:1000):: d_g !contiene la densita’
REAL(KIND=num), DIMENSION (1:1000):: f_r !funzione di ripartizione
REAL(KIND=num), DIMENSION (1:1000):: x !ascissa
REAL(KIND=num), DIMENSION (1:1000):: h !tasso di avaria
CHARACTER(len=10):: chiusura !serve a chiudere il programma
INTEGER:: i,j,k !indici di ciclo
INTEGER:: ind !il valore di alpha in cui calcolo la funz. gamma

!calcolo della funzione gamma fra 1 e 2

CALL sub_f_gamma_1_2 (f_gamma,alpha) !calcola la funzione gamma fra 1 e 2

WRITE(*,*) “This program calculates the gamma distribution, for alpha and lambda specified”
WRITE(*,*) “by the user. It plots the diagram of density, distribution function, failure rate”
WRITE(*,*) “and writes a file .txt with some values of the distribution function.”
WRITE(*,*) ” “
WRITE(*,*) “Please, type the value of alpha (a positive number between 0.01 and 25.0).”
WRITE(*,*) “Use dot as decimal separator.”
READ(*,*) alph
WRITE(*,*) ” “
WRITE(*,*) “Please, type the value of lambda (a positive number).”
WRITE(*,*) “Use dot as decimal separator.”
READ(*,*) lamb

!calcolo della funzione gamma e della ascissa fra 2 e 25

ciclo_2_5000: DO i=201,2500,1
alpha(i) = 1. + alpha(i-100)
f_gamma(i) = alpha(i-100)*f_gamma(i-100)
END DO ciclo_2_5000

!calcolo della funzione gamma e della ascissa fra 0 e 1

ciclo_0_1: DO i=2,100,1
alpha(i) = alpha(100+i)-1.
f_gamma(i) = f_gamma(100+i)/alpha(i)
END DO ciclo_0_1

!cerco il valore di alpha specificato dall’utente

ciclo_alph: DO i=1,2500,1
ind = i
IF (ABS(alpha(i)-alph)<=0.01) EXIT ciclo_alph
END DO ciclo_alph

!assegno il valore alla funzione gamma

f_g = f_gamma(ind)

WRITE(*,*) “The value of the gamma function in”, alph, “is:”, f_g

!calcolo la densita’ gamma

IF (alph <= 1) THEN
estremo = 2.
delta = estremo/1000.
ciclo_densita_1: DO k=1,1000,1
x(k)= delta*k
d_g(k) = ( ( lamb**alph )/( f_g ) )*( x(k)**( alph – 1. ) )*( e**( (-1)*lamb*x(k) ) )
END DO ciclo_densita_1
ELSE
estremo = ( (alph-1.)/lamb )*3.
delta = estremo/1000.
ciclo_densita_1: DO k=1,1000,1
x(k)= delta*(k-1.)
d_g(k) = ( ( lamb**alph )/( f_g ) )*( x(k)**( alph – 1. ) )*( e**( (-1)*lamb*x(k) ) )
END DO ciclo_densita_1
END IF

!calcolo la funzione di ripartizione

f_r(1) = 0.
f_r(3) = f_r(1) + delta*( d_g(1) + ( 4*d_g(2) ) + d_g(3) )/3.
f_r(2) = f_r(3)*0.5
ciclo_f_r: DO k=2,998,1
f_r(k+2) = f_r(k) + delta*( d_g(k) + ( 4*d_g(k+1) ) + d_g(k+2) )/3.;
END DO ciclo_f_r

!scrivo la funzione di ripartizione su un file di testo

OPEN(UNIT=20, FILE=’f_ripartizione.txt’, STATUS=’REPLACE’, ACTION=’WRITE’)

WRITE(20,*) “The gamma distribution for alpha =”, alph, “and lambda =”, lamb
WRITE(20,*) ” “
WRITE(20,*) “The value of the gamma function in”, alph, “is”, f_g
WRITE(20,*) ” ”
WRITE(20,*) “distribution function”
WRITE(20,*) ” “
WRITE(20,*) ” x “,”|”,”distribution function”
WRITE(20,*) “__________________________________”

ciclo_write: DO j=1,1000,1
WRITE(20,100) x(j),”|”,f_r(j)
END DO ciclo_write
100 FORMAT (1X,1F8.4,A,1F8.4)

!calcolo il tasso di avaria

ciclo_avaria: DO k=1,1000,1
h(k) = d_g(k)/( 1. – f_r(k) )
END DO ciclo_avaria

!disegno il grafico della densita’

CALL diagramma (x, d_g, ‘Density ‘)

!disegna il grafico della funzione di ripartizione

CALL diagramma (x, f_r, ‘distribution function’)

!disegno il grafico del tasso di avaria

CALL diagramma (x, h, ‘Failure rate ‘)

WRITE(*,*) “To close this program press any key, then press RETURN.”
READ(*,*) chiusura

END PROGRAM main_legge_gamma_avaria

!Data: 1/03/2018

MODULE mod_legge_gamma_avaria

!sezione dichiarativa

IMPLICIT NONE

!costanti

REAL,PARAMETER:: e = 2.71828182 !numero di Nepero
REAL,PARAMETER:: d_alpha = 0.01 !incremento delle ascisse della funzione gamma
INTEGER,PARAMETER:: num = 1 !aumento il numero di cifre significative

CONTAINS

!————————————————————————————-

SUBROUTINE sub_f_gamma_1_2 (f_gamma, alpha)

!Esegue l’integrale necessario al calcolo della funzione gamma
!per valori compresi fra 1 e 2 e registra questi valori in un array.
!L’integrazione viene eseguita con il metodo di Cavalieri-Simpson

IMPLICIT NONE

!dichiaro gli argomenti fittizi

REAL(KIND=num),INTENT(OUT),DIMENSION(:):: f_gamma !valori della funz. gamma in [1,2]
REAL(KIND=num),INTENT(OUT),DIMENSION(:):: alpha !contiene l’ascissa della funzione gamma

!dichiaro le variabili locali

REAL(KIND=num),DIMENSION(1:1000000):: funzione !è la funzione integranda
REAL(KIND=num),DIMENSION(1:1000000):: primitiva !è l’integrale illimitato
REAL(KIND=num),PARAMETER:: d_x = 0.001 !intervallo di integrazione
INTEGER:: i,j !indice del ciclo
INTEGER:: index = 1 !contiene l’estremo superiore della integrazione numerica

!sezione esecutiva

!assegno i valori di funzione e ascissa, eseguo l’integrazione numerica, assegno i valori di f_gamma_1_2

alpha(101) = 1. !assegna il primo valore di alpha
alpha(102) = alpha(101)+d_alpha !assegna il primo valore di alpha
f_gamma(101) = 1. !assegna il primo valore di gamma

ciclo_esterno: DO j=102,200,1 !assegna i valori a f_gamma_1_2

!assegna i valori a funzione

ciclo_funzione: DO i=1,1000000,1

funzione (i) = ( (d_x*REAL(i))**(alpha(j)-1.) )*( e**((-1)*d_x*REAL(i)) )
index = i
IF (funzione(i)<=0.0001) EXIT ciclo_funzione

END DO ciclo_funzione

!assegna i valori a primitiva

primitiva (1) = 0.
primitiva (1+2) = primitiva (1) + ( ( funzione(1+2) + funzione(1+1)*4. + funzione(1) )*(d_x/3.) )
primitiva (2) = primitiva (3)/2.

ciclo_integrazione: DO i=2,index-2,1 !integrazione numerica

primitiva (i+2) = primitiva (i) + ( ( funzione(i+2) + funzione(i+1)*4. + funzione(i) )*(d_x/3.) )

END DO ciclo_integrazione

f_gamma(j) = primitiva(index) !assegna il valore j-mo a f_gamma_1_2

alpha (j+1) = alpha (j) + d_alpha

END DO ciclo_esterno

RETURN

END SUBROUTINE sub_f_gamma_1_2

!————————————————————————————-

SUBROUTINE diagramma (x, funzione,label_y)

!disegno il grafico di funzione

!dichiaro gli argomenti fittizi

REAL(KIND=num),INTENT(IN),DIMENSION(1:1000):: x !ascissa
REAL(KIND=num),INTENT(IN),DIMENSION(1:1000):: funzione !funzione
CHARACTER (len=20), INTENT(IN):: label_y !nome delle ordinate

!dichiaro le variabili locali

REAL(KIND=num),DIMENSION(1:1000):: max_x !estremo ascissa
REAL(KIND=num),DIMENSION(1:1000):: max_funz !estremo funzione

!sezione esecutiva

max_x = MAXVAL(x) !valore massimo dlla ascissa
max_funz = MAXVAL(funzione) !valore massimo della ordinata

CALL METAFL (‘PDF’) !scelgo il formato dell’output
CALL SETPAG (‘DA4P’) !formato della pagina
CALL SCRMOD (‘revers’) !scritta nera su fondo bianco

CALL DISINI !inizializza DISLIN impostando parametri di default

CALL PAGERA !chiedo una linea di contorno per il piano xy
CALL DUPLX !chiedo un font a doppio spessore

CALL AXSPOS (300,1400) !coordinate dell’angolo in basso a sinistra
CALL AXSLEN (1500,1000) !la lunghezza, in pixel, dei due assi

CALL LABDIG (2,’X’) !numero di cifre decimali negli assi
CALL LABDIG (3,’Y’) !numero di cifre decimali negli assi

CALL NAME (‘X’,’x’) !assegno un nome all’asse delle x
CALL NAME (label_y,’y’) !assegno un nome all’asse delle y

CALL GRAF (0.,max_x,0.,max_x/5.,0.,max_funz+max_funz/5.,0.,max_funz/5.) !ambiti degli assi
CALL TEXMOD (‘ON’) !modalita’ TeX per le formule
CALL GRID (2,2) !impone una griglia
CALL CURVE (x,funzione,1000)

CALL DISFIN !termina DISLIN

END SUBROUTINE diagramma

!————————————————————————————-

END MODULE mod_legge_gamma_avaria

A table for the gamma function in [0.01, 3.99] is reported below. In each column, the value of α is indicated on the left while the corresponding value for Γ(α) is indicated on the right.

0.0100| 99.3302   0.0200| 49.3945 0.0300| 32.7556 0.0400| 24.4404 0.0500| 19.4548 0.0600| 16.1337 0.0700| 13.7639 0.0800| 11.9888 0.0900| 10.6096 0.1000|  9.5079 0.1100|  8.6081 0.1200|  7.8592 0.1300|  7.2268 0.1400|  6.6857 0.1500|  6.2175 0.1600|  5.8089 0.1700|  5.4490 0.1800|  5.1300 0.1900|  4.8450 0.2000|  4.5893 0.2100|  4.3585 0.2200|  4.1492 0.2300|  3.9587 0.2400|  3.7845 0.2500|  3.6246 0.2600|  3.4776 0.2700|  3.3418 0.2800|  3.2161 0.2900|  3.0994 0.3000|  2.9909 0.3100|  2.8898 0.3200|  2.7952 0.3300|  2.7067 0.3400|  2.6237 0.3500|  2.5457 0.3600|  2.4723 0.3700|  2.4031 0.3800|  2.3378 0.3900|  2.2762 0.4000|  2.2178 0.4100|  2.1625 0.4200|  2.1101 0.4300|  2.0602 0.4400|  2.0129 0.4500|  1.9679 0.4600|  1.9249 0.4700|  1.8840 0.4800|  1.8451 0.4900|  1.8078 0.5000|  1.7722 0.5100|  1.7382 0.5200|  1.7056 0.5300|  1.6744 0.5400|  1.6446 0.5500|  1.6159 0.5600|  1.5884 0.5700|  1.5621 0.5800|  1.5367 0.5900|  1.5124 0.6000|  1.4890 0.6100|  1.4665 0.6200|  1.4449

0.6300|  1.4240   0.6400|  1.4040 0.6500|  1.3846 0.6600|  1.3660 0.6700|  1.3480 0.6800|  1.3307 0.6900|  1.3140 0.7000|  1.2979 0.7100|  1.2823 0.7200|  1.2673 0.7300|  1.2528 0.7400|  1.2388 0.7500|  1.2253 0.7600|  1.2122 0.7700|  1.1996 0.7800|  1.1873 0.7900|  1.1755 0.8000|  1.1641 0.8100|  1.1530 0.8200|  1.1424 0.8300|  1.1320 0.8400|  1.1220 0.8500|  1.1124 0.8600|  1.1030 0.8700|  1.0939 0.8800|  1.0852 0.8900|  1.0767 0.9000|  1.0685 0.9100|  1.0606 0.9200|  1.0529 0.9300|  1.0455 0.9400|  1.0383 0.9500|  1.0313 0.9600|  1.0246 0.9700|  1.0181 0.9800|  1.0118 0.9900|  1.0058 1.0000|  1.0000 1.0100|  0.9933 1.0200|  0.9879 1.0300|  0.9827 1.0400|  0.9776 1.0500|  0.9727 1.0600|  0.9680 1.0700|  0.9635 1.0800|  0.9591 1.0900|  0.9549 1.1000|  0.9508 1.1100|  0.9469 1.1200|  0.9431 1.1300|  0.9395 1.1400|  0.9360 1.1500|  0.9326 1.1600|  0.9294 1.1700|  0.9263 1.1800|  0.9234 1.1900|  0.9206 1.2000|  0.9179 1.2100|  0.9153 1.2200|  0.9128 1.2300|  0.9105 1.2400|  0.9083

1.2500|  0.9062   1.2600|  0.9042 1.2700|  0.9023 1.2800|  0.9005 1.2900|  0.8988 1.3000|  0.8973 1.3100|  0.8958 1.3200|  0.8945 1.3300|  0.8932 1.3400|  0.8920 1.3500|  0.8910 1.3600|  0.8900 1.3700|  0.8892 1.3800|  0.8884 1.3900|  0.8877 1.4000|  0.8871 1.4100|  0.8866 1.4200|  0.8862 1.4300|  0.8859 1.4400|  0.8857 1.4500|  0.8855 1.4600|  0.8855 1.4700|  0.8855 1.4800|  0.8856 1.4900|  0.8858 1.5000|  0.8861 1.5100|  0.8865 1.5200|  0.8869 1.5300|  0.8874 1.5400|  0.8881 1.5500|  0.8888 1.5600|  0.8895 1.5700|  0.8904 1.5800|  0.8913 1.5900|  0.8923 1.6000|  0.8934 1.6100|  0.8946 1.6200|  0.8958 1.6300|  0.8971 1.6400|  0.8985 1.6500|  0.9000 1.6600|  0.9016 1.6700|  0.9032 1.6800|  0.9049 1.6900|  0.9067 1.7000|  0.9085 1.7100|  0.9105 1.7200|  0.9125 1.7300|  0.9146 1.7400|  0.9167 1.7500|  0.9190 1.7600|  0.9213 1.7700|  0.9237 1.7800|  0.9261 1.7900|  0.9287 1.8000|  0.9313 1.8100|  0.9340 1.8200|  0.9367 1.8300|  0.9396 1.8400|  0.9425 1.8500|  0.9455 1.8600|  0.9486

1.8700|  0.9517   1.8800|  0.9550 1.8900|  0.9583 1.9000|  0.9617 1.9100|  0.9651 1.9200|  0.9687 1.9300|  0.9723 1.9400|  0.9760 1.9500|  0.9798 1.9600|  0.9836 1.9700|  0.9876 1.9800|  0.9916 1.9900|  0.9957 2.0000|  1.0000 2.0100|  1.0032 2.0200|  1.0076 2.0300|  1.0121 2.0400|  1.0167 2.0500|  1.0214 2.0600|  1.0261 2.0700|  1.0309 2.0800|  1.0358 2.0900|  1.0408 2.1000|  1.0459 2.1100|  1.0510 2.1200|  1.0563 2.1300|  1.0616 2.1400|  1.0670 2.1500|  1.0725 2.1600|  1.0781 2.1700|  1.0838 2.1800|  1.0896 2.1900|  1.0955 2.2000|  1.1014 2.2100|  1.1075 2.2200|  1.1136 2.2300|  1.1199 2.2400|  1.1263 2.2500|  1.1327 2.2600|  1.1393 2.2700|  1.1459 2.2800|  1.1526 2.2900|  1.1595 2.3000|  1.1665 2.3100|  1.1735 2.3200|  1.1807 2.3300|  1.1880 2.3400|  1.1953 2.3500|  1.2028 2.3600|  1.2104 2.3700|  1.2181 2.3800|  1.2260 2.3900|  1.2339 2.4000|  1.2420 2.4100|  1.2501 2.4200|  1.2584 2.4300|  1.2668 2.4400|  1.2754 2.4500|  1.2840 2.4600|  1.2928 2.4700|  1.3017 2.4800|  1.3107

2.4900|  1.3199   2.5000|  1.3292 2.5100|  1.3386 2.5200|  1.3481 2.5300|  1.3578 2.5400|  1.3676 2.5500|  1.3776 2.5600|  1.3877 2.5700|  1.3979 2.5800|  1.4083 2.5900|  1.4188 2.6000|  1.4294 2.6100|  1.4403 2.6200|  1.4512 2.6300|  1.4623 2.6400|  1.4736 2.6500|  1.4850 2.6600|  1.4966 2.6700|  1.5083 2.6800|  1.5202 2.6900|  1.5323 2.7000|  1.5445 2.7100|  1.5569 2.7200|  1.5694 2.7300|  1.5822 2.7400|  1.5951 2.7500|  1.6082 2.7600|  1.6214 2.7700|  1.6349 2.7800|  1.6485 2.7900|  1.6623 2.8000|  1.6763 2.8100|  1.6905 2.8200|  1.7049 2.8300|  1.7194 2.8400|  1.7342 2.8500|  1.7492 2.8600|  1.7644 2.8700|  1.7797 2.8800|  1.7953 2.8900|  1.8111 2.9000|  1.8271 2.9100|  1.8434 2.9200|  1.8598 2.9300|  1.8765

2.9400|  1.8934   2.9500|  1.9106 2.9600|  1.9279 2.9700|  1.9455 2.9800|  1.9634 2.9900|  1.9815 3.0000|  2.0000 3.0100|  2.0165 3.0200|  2.0354 3.0300|  2.0547 3.0400|  2.0741 3.0500|  2.0938 3.0600|  2.1138 3.0700|  2.1340 3.0800|  2.1545 3.0900|  2.1753 3.1000|  2.1963 3.1100|  2.2177 3.1200|  2.2393 3.1300|  2.2612 3.1400|  2.2835 3.1500|  2.3059 3.1600|  2.3288 3.1700|  2.3519 3.1800|  2.3753 3.1900|  2.3991 3.2000|  2.4231 3.2100|  2.4476 3.2200|  2.4723 3.2300|  2.4974 3.2400|  2.5228 3.2500|  2.5486 3.2600|  2.5747 3.2700|  2.6012 3.2800|  2.6280 3.2900|  2.6552 3.3000|  2.6828 3.3100|  2.7109 3.3200|  2.7392 3.3300|  2.7680 3.3400|  2.7971 3.3500|  2.8266 3.3600|  2.8566 3.3700|  2.8870 3.3800|  2.9178

3.3900|  2.9490   3.4000|  2.9807 3.4100|  3.0128 3.4200|  3.0454 3.4300|  3.0784 3.4400|  3.1119 3.4500|  3.1459 3.4600|  3.1803 3.4700|  3.2152 3.4800|  3.2506 3.4900|  3.2865 3.5000|  3.3229 3.5100|  3.3598 3.5200|  3.3973 3.5300|  3.4352 3.5400|  3.4737 3.5500|  3.5128 3.5600|  3.5524 3.5700|  3.5926 3.5800|  3.6333 3.5900|  3.6746 3.6000|  3.7165 3.6100|  3.7591 3.6200|  3.8022 3.6300|  3.8459 3.6400|  3.8903 3.6500|  3.9353 3.6600|  3.9809 3.6700|  4.0272 3.6800|  4.0742 3.6900|  4.1218 3.7000|  4.1702 3.7100|  4.2192 3.7200|  4.2689 3.7300|  4.3194 3.7400|  4.3705 3.7500|  4.4225 3.7600|  4.4751 3.7700|  4.5286 3.7800|  4.5828 3.7900|  4.6378 3.8000|  4.6936 3.8100|  4.7503 3.8200|  4.8077 3.8300|  4.8660

3.8400|  4.9251   3.8500|  4.9852 3.8600|  5.0461 3.8700|  5.1079 3.8800|  5.1706 3.8900|  5.2342 3.9000|  5.2987 3.9100|  5.3642 3.9200|  5.4307 3.9300|  5.4982 3.9400|  5.5667 3.9500|  5.6361 3.9600|  5.7067 3.9700|  5.7782 3.9800|  5.8508 3.9900|  5.9245