#!/usr/bin/env python
# coding: utf-8

# <img src="https://risk-engineering.org/static/img/logo-RE.png" width="100" alt="" style="float:right;margin:15px;">
# 
# This notebook is an element of the [risk-engineering.org courseware](https://risk-engineering.org/). It can be distributed under the terms of the [Creative Commons Attribution-ShareAlike licence](https://creativecommons.org/licenses/by-sa/4.0/).
# 
# ---
# 
# This notebook contains an introduction to the use of [NumPy](https://numpy.org/) for numerical integration using the Monte Carlo method.

# # Monte Carlo methods for numerical integration

# This notebook contains an illustration of the use of Monte Carlo methods for numerical integration. See the [associated course materials](https://risk-engineering.org/monte-carlo-methods/) for an introduction to the use of stochastic simulation methods and to download this content as a Jupyter/Python notebook. Numerical integration is often used to evaluate risk measures in the finance industry. 
# 
# ## First integral 
# 
# **Task**: resolve the integral 
# 
# $\int_1^5 x^2 dx$

# ### Analytical method

# We start by resolving this integral using the standard analytical method, assisted by the [SymPy symbolic mathematics library](https://sympy.org/). 

# In[1]:


import sympy
from sympy.interactive import printing
printing.init_printing(use_latex="mathjax")

x = sympy.Symbol("x")
sympy.integrate(x**2, (x, 1, 5))


# In[2]:


# note that _ references the value of the previous cell
float(_)


# ### Numerical method

# The expression we wish to integrate is very simple here, so calculating its integral analytically is easy (even without resorting to Python’s symbolic mathematics package!). In many cases, however, analytical approaches to integration are not feasible:
# 
# - the expression we wish to integrate is very complicated, possibly without a closed analytical integral
# - it is only known in “black box” form (for instance a software module): we can evaluate it at various points but don’t know the corresponding equation
# 
# In such situations, stochastic simulation (“Monte Carlo”) methods allow us to generate an approximation of the integral, simply by evaluating the expression a large number of times at randomly selected points in the input space and counting the proportion that are less than the integrand at that point. The larger the number of simulations we run, the better the approximation (and note that computers are very very fast today, so on simple problems the simulation will run in a blink of an eye!). 

# In[3]:


import numpy


# In[4]:


N = 100_000
accum = 0
for i in range(N):
    x = numpy.random.uniform(1, 5)
    accum += x**2
measure = 5 - 1
measure * accum/float(N)


# ## Second integral
# 
# **Task**: resolve the integral $\int_3^2 x^2 + 4 x sin(x)$. 

# ### Analytical solution

# In[5]:


x = sympy.Symbol("x")
float(sympy.integrate(x**2 + 4 * x * sympy.sin(x), (x, 2, 3)))


# ### Stochastic solution

# In[6]:


N = 100_000
accum = 0 
for i in range(N):
    x = numpy.random.uniform(2, 3)
    accum += x**2 + 4*x*numpy.sin(x)
measure = 3 - 2
measure * accum/float(N)


# **Exercise**: now undertake the same comparison of analytical and stochastic simulation methods to evaluate the integral
# 
# $$\int_1^3 e^{x^2}$$

# (Hint: the result should be approximately 1464.2)

# ## A 2D integral

# Now we move to a 2D integral, $\int\int cos(x^4) + 3y^2 dx dy$ over $x ∈ [4,6]$ and $y ∈ [0, 1]$. 
# 
# Monte Carlo integration methods become increasingly attractive when compared to other numerical integration methods (such as quadrature) as the number of dimensions increases. 

# Let’s examine the **analytical solution**, again thanks to sympy. 

# In[7]:


x = sympy.Symbol("x")
y = sympy.Symbol("y")
float(sympy.integrate(sympy.cos(x**4) + 3 * y**2, (x, 4, 6), (y, 0, 1)))


# Now compare with **Monte Carlo estimation**:

# In[8]:


N = 100_000
accum = 0
for i in range(N):
    x = numpy.random.uniform(4, 6)
    y = numpy.random.uniform(0, 1)
    accum += numpy.cos(x**4) + 3 * y * y 
measure = 2 * 1
measure * accum/float(N)


# In[ ]: