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:
- for α>1 the probability of failure of a mechanical component (or death of an organism) increases as the time passes by (increasing hazard rate);
- for α=1 the probability of failure is constant (constant hazard rate);
- 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_gammaUSE DISLIN !libreria grafica
USE mod_f_gamma !modulo con procedureIMPLICIT 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 cicloCALL 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(*,*) chiusuraEND 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 significativeCONTAINS
!————————————————————————————-
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-SimpsonIMPLICIT 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 gammaciclo_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,1funzione (i) = ( (d_x*REAL(i))**(alpha(j)-1.) )*( e**((-1)*d_x*REAL(i)) )
index = i
IF (funzione(i)<=0.0001) EXIT ciclo_funzioneEND 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_integrazionef_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_gammaRETURN
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 procedureIMPLICIT 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(*,*) chiusuraEND 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 precisioneCONTAINS
!————————————————————————————-
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-SimpsonIMPLICIT 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 gammaciclo_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,1funzione (i) = ( (d_x*REAL(i))**(alpha(j)-1.) )*( e**((-1)*d_x*REAL(i)) )
index = i
IF (funzione(i)<=0.0001) EXIT ciclo_funzioneEND 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_integrazionef_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 ordinataCALL 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 defaultCALL PAGERA !chiedo una linea di contorno per il piano xy
CALL DUPLX !chiedo un font a doppio spessoreCALL AXSPOS (300,1400) !coordinate dell’angolo in basso a sinistra
CALL AXSLEN (1500,1000) !la lunghezza, in pixel, dei due assiCALL LABDIG (2,’X’) !numero di cifre decimali negli assi
CALL LABDIG (3,’Y’) !numero di cifre decimali negli assiCALL NAME (‘X’,’x’) !assegno un nome all’asse delle x
CALL NAME (label_y,’y’) !assegno un nome all’asse delle yCALL 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 procedureIMPLICIT 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(*,*) chiusuraEND 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 significativeCONTAINS
!————————————————————————————-
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-SimpsonIMPLICIT 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 gammaciclo_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,1funzione (i) = ( (d_x*REAL(i))**(alpha(j)-1.) )*( e**((-1)*d_x*REAL(i)) )
index = i
IF (funzione(i)<=0.0001) EXIT ciclo_funzioneEND 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_integrazionef_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 ordinataCALL METAFL (‘PDF’) !scelgo il formato dell’output
CALL SETPAG (‘DA4P’) !formato della pagina
CALL SCRMOD (‘revers’) !scritta nera su fondo biancoCALL 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 spessoreCALL AXSPOS (300,1400) !coordinate dell’angolo in basso a sinistra
CALL AXSLEN (1500,1000) !la lunghezza, in pixel, dei due assiCALL LABDIG (2,’X’) !numero di cifre decimali negli assi
CALL LABDIG (3,’Y’) !numero di cifre decimali negli assiCALL NAME (‘X’,’x’) !assegno un nome all’asse delle x
CALL NAME (label_y,’y’) !assegno un nome all’asse delle yCALL 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 |
paolo maccallini 