Saturday, August 12, 2023

HP Prime: Monte Carlo Integration

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.