Monte Carlo Method Using Octave GNU Tool
1. Abbreviation
2. Introduction
The Monte Carlo method is a numerical method for solving mathematical problems by modeling random variables using pseudo random numbers. I used the Octave Tool for my calculations. See the MonteCarlo.m
Here I will give some information from the theory of probability, which I will use in the future. The law of large numbers: let there be a random variable ξ and N measurements/samples of this random variable are carried out:
then the mathematical expectation is calculated:
in this case the calculation error (4) is determined by the standard deviation:
For N => +∞, then (4) gives the exact value of the mathematical expectation.
Note:
The variance (1) is calculated using the mathematical statistics formula:
The central limit theorem of probability theory: If a random variable depends on a large number of random factors, then it is described by a distribution close to normal (Gaussian) distribution. For example, let there be N random variables with the same distribution law, then the sum of these random variables is described by a normal distribution for N => +∞:
Then the probability is within:
The Probability Density Function of the normal (Gaussian) distribution is:
Figure 2-1: Probability Density Function of the normal (Gaussian) Distribution with the σ=1, E(ξ)=0
The Octave/Matlab Tool function is used to create an array of pseudo random numbers with a normal distribution σ=1, E(ξ)=0:
Gaussian_Sigma1_Expected0 = randn(1, Size, “double”); (12)
where
Size – random array size
“double” – double float (64 bits) or single float (32 bits)
For arbitrary values σ, E(ξ) the formula is used:
Random_SigmaValue_ExpectedValue = SigmaValue*randn(1, Size, “double”) + ExpectedValue; (13)
To calculate the probabilities of the normal distribution, the Gauss error function is used, which is supported by Octave/Matlab Tools:
If there is a normal distribution with E(ξ) = a and standard deviation σ, then the probability that the random variable takes values is -x…x:
Let’s calculate the function (15) using the Octave Tool for -3σ … 3σ. See the MonteCarlo.m:
Conclusion:
For a normal distribution, all values of a random variable with almost a single probability lie within: -3σ … 3σ. This fact is actively used for modeling and is called the three sigma rule.
The Monte Carlo method often uses an uniform distribution. Let’s first consider an uniform distribution for a random variable given in an open interval (0,1) with a probability density function:
f(x) = 1 if 0 < x < 1
f(x) = 0 if x ≤ 0 or x ≥ 1 (17)
Expected value:
Variance:
Figure 2-2: Probability Density Function of the Uniform Distribution
To create an array of pseudo random numbers of the uniform distribution (0,1), the Octave/Matlab Tool function is used:
UniformDistribution_0_1 = rand(1, Size, “double”); (20)
where
Size – random array size
“double” – double float with 64 bits or single float with 32 bits
Pseudo random numbers of the uniform distribution over the interval (a,b) is calculated by the Octave/Matlab Tool:
UniformDistribution_a_b = (b – a)*rand(1, Size, “double”) + a; (21)
Probability Density Function of the uniform distribution over the interval (a,b):
f(x) = 1/(b-a) if a < x < b
f(x) = 0 if x ≤ a or x ≥ b (22)
Expected value:
Variance:
.
3. Simulation of a Voltage Divider with the Monte Carlo Method
Let there be a voltage divider. See the Figure 3-1:
Figure 3-1: Voltage Divider
Scheme Parameters:
Let’s model this circuit and estimate the impedance Rall = R1 + R2 and the output voltage Uout using arrays of pseudo random numbers of size N=100. (25) contains three random variables having a normal distribution.
R1:
R2:
Uin:
Next using the Octave Tool the pseudo random arrays (26), (27), (28) are generated. Then the Rall and Uout arrays are calculated using the formulas:
Then according to the formulas (4), (6), (2) expected values, variance values and standard deviation values of the Rall and Uout are calculated. Using the three sigma rule, the relative error in percent is calculated. See the MonteCarlo.m
One of the simulation results at N=100:
Using the formula (5), it is possible to estimate the accuracy of the E(Rall) and E(Uout) simulation:
Rall:
Uout:
Each new simulation always gives a different result, but the standard deviations must correspond to (32), (33). To increase the accuracy of the simulation, it is necessary to increase N.
Note
Appendix 1 provides another estimation method for (31) by calculating the differentials of the function, but this method gives a pessimistic less accurate result. I recommend using the Monte Carlo method for this.
.
4. Conclusions
-
The evaluation of processes/models influenced by a variety of random factors is solved efficiently using the Monte Carlo method
-
The Monte Carlo method can be used not only for simulating models in which random factors of influence are present, but also for solving differential equations, calculating certain integrals, etc. Appendix 2 provides a simple example of the integral calculating. Unfortunately, the accuracy of these calculations is often low
.
5. Octave GNU file MonteCarlo.m
The script generates plots of the probability density functions of the normal and uniform random distributions, as well as simulates a probabilistic Voltage Divider model and calculates a definite integral using the Monte Carlo method.
.
6. Download the MonteCarlo.m
You can download the files:
MonteCarlo.m
with the button:
.
7. Literature / References
[1] I. M. Sobol “Numerical Monte Carlo Methods“, Publishing House: Nauka, Moscow, 1973
[2] B. R. Levin “Theoretical foundations of statistical radio engineering“, First Book, Publishing House “Sovietskoe radio“, Moscow, 1969
[3] B. V. Gnedenko, A. Y. Hinchin “An elementary introduction to probability theory”, Publishing House: Nauka, Moscow, 1970
.
Appendix 1. Voltage Divider using the Estimate of the Function Differential
Let there be a function:
The differential from this function:
Moving from differential to delta:
Note:
Such a transition is possible if the delta value is small and the nonlinearity of the function f() can be neglected.
Now let’s go back to the Voltage Divider example. See the section 3.
Rall = R1 + R2 (A1-4)
Then from (A1-3):
ΔRall = ΔR1 + ΔR2 (A1-5)
From (26), (27):
ΔRall = 20Ω + 5Ω = 25Ω
(ΔRall/Rall) 100% = (25Ω/150Ω) 100% = 16.667% (A1-6)
Let’s make an estimate of the output voltage:
Then from (A1-3):
From (26), (27), (28):
Note:
The obtained results of deviation (A1-6), (A1-9) from the mean are significantly higher compared to the result of the Monte Carlo method (31). But it is the result obtained by the Monte Carlo method that is more realistic!
.
Appendix 2. Calculation of the definite integral by the Monte Carlo Method
Let there be a problem of calculating a definite integral:
Let there be a random variable ξ given by (a,b) with a probability density function: f(x). Now let’s define a new random variable:
The expected value of the η is calculated by the formula:
On the other hand the expected value of the η can be calculated by N times simulation of the random variable ξ:
Then:
Thus the integral (A2-1) can be calculated by the formula (A2-5) using an arbitrary probability distribution f(x). The easiest way is to take a uniform distribution on (a,b). See the (22).
Let’s estimate the accuracy of the calculation by the formula (A2-5). This sum can be interpreted as the sum of N independent random variables. It follows from the Central limit theorem of probability theory that the sum for large N gives a normal process for which the three sigma rule applies, where sigma is calculated by the formula (5):
Note
In this case V(η) can be calculated by the formula (6)
Using the three sigma rule we get a relative error:
Let’s consider an example of calculating the integral:
Note:
The exact solution of the integral:
Let there be a random variable ξ given by (0,π) with a probability density function of an uniform distribution:
Then from (A2-2):
And the integral (A2-8) will be written as:
This integral can be calculated using the formula (A2-5):
And using the Octave Tool:
Let’s estimate the accuracy of the integral calculation (A2-14) by the formula (A2-6) and (A2-7) using the Octave Tool:
