HP Prime: Monte Carlo Integration
Integration by Random Points
There are many ways to calculate the area under the curve f(x) in the interval [a, b]. Two of the common are calculating the antiderivative and using the Simpson's Rule. Another method is the Monte Carlo method. Unlike Simpson's Rule, where the intervals are fixed, the Monte Carlo method picks random points in interval. The integral is calculated using the simple formula:
∫ f(x) dx ≈ (b - a) ÷ (n - 1) * Σ( f(x_i), from i = 0 to n)
x_i is a random point in the interval [a, b]
n = number of points, n > 2
How are the random numbers are picked are determined by various methods: pseudo-random generators, using a uniform distribution, using a normal distribution, etc.
The program MONTE for the HP Prime uses the RANDOM function. The HP Prime allows to pick a random (real) number in the interval [a , b] by the syntax RANDOM(a,b).
Note: I use the Function app function AREA to calculate the integral using any of the system variables F0 through F9. In trade, the program uses global variables instead of local variables.
HP Prime Program: MONTE
Code:
EXPORT MONTE()
BEGIN
// 2023-06-03 EWS
STARTAPP("Function");
// use input box for global variables
INPUT({A,B,E},"Monte Carlo of F1(X)",
{"lower:","upper:","# places:"});
// set radians
HAngle:=0;
// input check
IF A≥B THEN
MSGBOX("Lower limit must be
less than upper limit.");
KILL;
END;
// AREA must be stored to a global var.
I:=Function.AREA(F1(X),A,B);
// initialize variables
N:=1;
S:=0;
// repeat loop
REPEAT
X:=RANDOM(A,B);
S:=S+F1(X);
N:=N+1;
IF N>1 THEN
J:=(B-A)*S/(N-1);
END;
UNTIL (N>1) AND (ABS(I-J)<ALOG(−E));
// results
PRINT();
PRINT("Results");
PRINT("Actual: "+STRING(I));
PRINT("Approx: "+STRING(J));
PRINT("# terms used: "+STRING(N));
// return home
STARTVIEW(−1);
END;
Examples
Estimate the integral of F1(X) to three decimal places. Your mileage may vary.
∫ F1(X) dX from X = A to X = B
Example 1:
F1(X) = X^2 + 1
A = 0, B = 3
Exact Integral = 12
Trial 1: 12.0004750664, 64 points
Trial 2: 11.9996447323, 272 points
Trial 3: 12.0001030297, 9 points
Example 2:
F1(X) = COS(X - 1)
A = 1, B = π/2 + 1
Exact Integral = 1
Trial 1: 1.00079319497, 411 points
Trial 2: .999705900656, 187 points
Trial 3: 1.00071067474, 737 points
Example 3:
F1(X) = SIN X/X
A = 0, B = 4
Exact Integral = Si(4) ≈ 1.75820313895
Trial 1: 1.75789723089, 205 points
Trial 2: 1.75729880586, 133 points
Trial 3: 1.75832553216, 223 points
Notes
While the Monte Carlo method is easy to calculate, it is difficult to get an accurate answer. The method requires a lot of calculation points, and how many really depends on what random numbers are picked. It's really the luck of a draw. This is good for a short approximation but I recommend the Simpson's Rule, Trapezoid Rule, or when possible and feasible, finding the antiderivative instead.
Source
Cumer, Victor. "The basics of Monte Carlo integration" Towards Data Science. Medium. October 26, 2020. Last Retrieved June 4, 2023. https://towardsdatascience.com/the-basics-of-monte-carlo-integration-5fe16b40482d
Eddie
All original content copyright, © 2011-2023. Edward Shore. Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited. This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.