Monte Carlo Simulation to estimate e

[Quasar Chunawala]

Hi all,

Some of you might already know, how to write a Monte-Carlo simulation to estimate [imath]\pi[/imath]. We use the uniform random number generator to produce a pair of pseudo-random numbers [imath](x_i,y_i)[/imath] in [imath][0,1][/imath]. After [imath]N[/imath] draws, we find the proportion of points lying in the unit circle, that is, those with [imath]x_i^2 + y_i^2 \leq 1[/imath] and divide this by the total number of draws. If [imath]N[/imath] is large, this ratio approaches the geometrical probability of the event that a random point drawn from the unit square lies in a unit quarter-circle: [imath]\pi/4[/imath].

While browsing Math SE, I recently came across a somewhat more challenging puzzle/tasks. How do you estimate the Euler's constant [imath]e \approx 2.71828\ldots [/imath] using the same generator?

My initial idea : If [imath]X[/imath] is [imath]Expo(1)[/imath] random variable, its CDF is [imath]F_X(x) = 1 - e^{-x}[/imath]. By the inverse transform method, [imath]F_X^{-1}(U)[/imath] is an exponential random variable. I generated a sample of [imath]U[0,1][/imath] random numbers, and applied the inverse transform. Since, [imath]F_X(1)=1-e^{-1}[/imath], we can find a count of the number of sample points less than or equal to unity. Their proportion should approach [imath]1-e^{-1}[/imath].

import numpy as np
import math

def generate_e(N):
    sum = 0

    # Draw N samples from a uniform distribution
    U = np.random.uniform(low=1e-9,high=1.0,size=N)

    # Using the inverse transform technique, log(1/(\lambda(1-u)))~Expo(\lambda) distribution
    lambda_ = 1.0
    X = -np.log(lambda_ *(1.0 - U))

    X2 = X[X<=1.0]
    sum = len(X2)

    e = 1 / (1 - (sum / N))
    return e

estimate = generate_e(10**6)
print(estimate)

However, this approach is really bad (very slow convergence). Also, you aren't allowed to use log and circular functions. Here's the Stats SE thread, that discusses a much simpler and cool idea. (Also, I still need to prove that the random variable [imath]\xi[/imath] defined in the first answer to the OP's question has expectation [imath]e[/imath]. It'll make for an interesting exercise.)

 

[Daniel Duffy]

V2 u(x) = exp(x) satisfies ODE

du/dx = u
u(0) = 1

Use FDM and compute up to x = 1. QED
(Padé approximant, essentially)

https://www.researchgate.net/publication/228590312_Nineteen_Dubious_Ways_to_Compute_the_Exponential_of_a_Matrix_Twenty-Five_Years_Later

 

[Daniel Duffy]

And I have another 20 methods here

exp(5) = [$]e^5[$] - wilmott.com

forum.wilmott.com

And use compound interest (that's where e came from way back).

 

[Quasar Chunawala]

V2 u(x) = exp(x) satisfies ODE

du/dx = u
u(0) = 1

Use FDM and compute up to x = 1. QED
(Padé approximant, essentially)

https://www.researchgate.net/publication/228590312_Nineteen_Dubious_Ways_to_Compute_the_Exponential_of_a_Matrix_Twenty-Five_Years_Later

I looked up Pade's approximant on Wiki; so [imath]P_N^M(x)[/imath] is essentially a quotient of two polynomials of degree [imath]M,N[/imath] whose power series ties up upto the first [imath]M+N+1[/imath] terms with the power series of the function.

For the function [imath]f(x)=e^x[/imath], I computed

[math]
P_3^3(x) = \frac{1+\frac{1}{2}x + \frac{1}{10}x^2 + \frac{1}{120}x^3}{1-\frac{1}{2}x + \frac{1}{10}x^2 - \frac{1}{120}x^3}
[/math]

I get [imath]P_3^3(1) = 2.71830986[/imath], having a relative error [imath]1.03 \times 10^{-5}[/imath], still better than Taylor's series truncated upto terms of degree [imath]6[/imath]; a relative error [imath]8.32 \times 10^{-5}[/imath].

It's like the Pade's approximant of order [imath]3[/imath] is trying to kill [imath]6[/imath] terms of the Taylor series. Pretty impressive speedup!

Btw Daniel, the original puzzle was to use a non-deterministic way to approximate [imath]e[/imath].

 

[Quasar Chunawala]

V2 u(x) = exp(x) satisfies ODE

du/dx = u
u(0) = 1

Use FDM and compute up to x = 1. QED
(Padé approximant, essentially)

https://www.researchgate.net/publication/228590312_Nineteen_Dubious_Ways_to_Compute_the_Exponential_of_a_Matrix_Twenty-Five_Years_Later

Great nuggets, I will explore this. Btw, one of the authors appears to be Van Loan from the infamous matrix methods book.

 

[Daniel Duffy]

Great nuggets, I will explore this. Btw, one of the authors appears to be Van Loan from the infamous matrix methods book.

And Clive Moler is Mr. Matlab.

 

[Daniel Duffy]

I looked up Pade's approximant on Wiki; so [imath]P_N^M(x)[/imath] is essentially a quotient of two polynomials of degree [imath]M,N[/imath] whose power series ties up upto the first [imath]M+N+1[/imath] terms with the power series of the function.

For the function [imath]f(x)=e^x[/imath], I computed

[math]
P_3^3(x) = \frac{1+\frac{1}{2}x + \frac{1}{10}x^2 + \frac{1}{120}x^3}{1-\frac{1}{2}x + \frac{1}{10}x^2 - \frac{1}{120}x^3}
[/math]

I get [imath]P_3^3(1) = 2.71830986[/imath], having a relative error [imath]1.03 \times 10^{-5}[/imath], still better than Taylor's series truncated upto terms of degree [imath]6[/imath]; a relative error [imath]8.32 \times 10^{-5}[/imath].

It's like the Pade's approximant of order [imath]3[/imath] is trying to kill [imath]6[/imath] terms of the Taylor series. Pretty impressive speedup!

Btw Daniel, the original puzzle was to use a non-deterministic way to approximate [imath]e[/imath].

That is Pade(q, p), more generally.
Very important for PDE and ODE

Pade(0,1) implicit Euler
Pade(1,0) explicit Euler
Pade(1,1) Crank Nicolson

Gives a great foothold into FDM. You could run them as well to investigate accuracy O(dt^(p+q+1)) as dt -> 0.

 

[Daniel Duffy]

https://www.ma.imperial.ac.uk/~ajacquie/IC_Num_Methods/IC_Num_Methods_Docs/Literature/DuffyCN.pdf

Can compute e by FD schemes (20), (21), (22). How does accuracy compare?
(it's about 4 lines of code).

 

[Daniel Duffy]

Updated version

 

[Quasar Chunawala]

That is Pade(q, p), more generally.
Very important for PDE and ODE

Pade(0,1) implicit Euler
Pade(1,0) explicit Euler
Pade(1,1) Crank Nicolson

Gives a great foothold into FDM. You could run them as well to investigate accuracy O(dt^(p+q+1)) as dt -> 0.

Simple foward Euler and backward Euler

import numpy as np
import math
import matplotlib.pyplot as plt

plt.style.use('seaborn-poster')

f = lambda t, u: u          # ODE
h = 0.10                    # Step-size
A = 0.0
B = 1.0                     # Region
u_0 = 1                     # Initial condition

def forwardEuler(f,A,B,h,u_0):
    t = np.arange(A, B + h, h)  # Define the mesh points
    u = np.zeros(len(t))
    u[0] = u_0

    for k in range(0,len(t)-1):
        u[k+1] = u[k] + h*f(t[k],u[k])

    return u,t

def backwardEuler(f,A,B,h,u_0):
    t = np.arange(A, B + h, h)  # Define the mesh points
    u = np.zeros(len(t))
    u[0] = u_0

    for k in range(0,len(t)-1):
        u[k+1] = u[k]/(1-h)

    return u,t

The error is clearly [imath]O(h)[/imath].

Absolute error in fwd difference for h = 0.1, err = 0.124539368359045
Absolute error in fwd difference for h = 0.05, err = 0.06498412331462378
Absolute error in bkwd difference for h = 0.1, err = 0.14969016233339527
Absolute error in bkwd difference for h = 0.05, err = 0.07122798905721561

 

[Daniel Duffy]

Very nice.
The Crank Nicolcon Pade(1,1) is 2nd order and is used a _lot_ in finance.

Now, we can discretise Black Scholes pde(S,t) in S to get an ODE(t) and use the above ODE solvers and others (Boost C++ odeint, Python odeitnt).
aka Method of Lines (MOL).