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

Gamma distribution.png

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:

hazard rate.png

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:

gamma function.png

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.

Gamma function table.png
Table 1. Values of the gamma function in [0.1, 1.99], with a step of 0.01.
Gamma.JPG
Figure 1. The gamma function between 0 and 4, as calculated by my code in Fortran. Plotted by Octave.

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 α.

Gamma distribution plots.png
Figure 2. Densities (left) and distribution functions (right) for α = 1.5, 2.0, 3.5 and λ = 1.5.

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.

Cattura.PNG
Figure 3. Density and distribution function for α = 20 and λ=2.

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).

Hazard rate plot.png
Figure 4. Hazard rate for α = 0.8, λ = 1.5 (left) and α = 3.0, λ=1.5 (right).

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

 

 

 

Annunci

2 pensieri su “The Gamma distribution

Rispondi

Inserisci i tuoi dati qui sotto o clicca su un'icona per effettuare l'accesso:

Logo WordPress.com

Stai commentando usando il tuo account WordPress.com. Chiudi sessione /  Modifica )

Google+ photo

Stai commentando usando il tuo account Google+. Chiudi sessione /  Modifica )

Foto Twitter

Stai commentando usando il tuo account Twitter. Chiudi sessione /  Modifica )

Foto di Facebook

Stai commentando usando il tuo account Facebook. Chiudi sessione /  Modifica )

Connessione a %s...