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.
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.
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).
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 
Non ci capisco niente Paolo. Sono proprio analfabeta in questo campo. 😦
LikeLike