#  Code for Monte Carlo class, Jonathan Goodman
#  http://www.math.nyu.edu/faculty/goodman/teaching/MonteCarlo20/
#  ImportanceSamplingDemo.py
#  The author gives permission for anyone to use this publicly posted 
#  code for any purpose.  The code was written for teaching, not research
#  or commercial use.  It has not been tested thoroughly and probably has
#  serious bugs.  Results may be inaccurate, incorrect, or just wrong.

# Illustrate importance sampling for estimating moments of a Gaussian.
# See week 2 notes for an explanation in mathematical notation.
# The notation in the code is the same as the notation in the notes.

from numpy.random import Generator, PCG64    #  numpy randon number generator
import numpy as np
  

seed = 12344             # change this to get different answers
bg   = PCG64(seed)       # instantiate a bit generator
rg   = Generator(bg)     # instantiate a random number generator with ...
                         # ... this bit generator
n   = 4                  # don't hard wire constants
tn  = 2*n                # tn for "two times n"
N = 10000                # number of samples per MC estimate

#         Calclate the exact moment, formula from the wee2 notes
A = 1.
for k in range(1,(n+1)):
    A = A * (2*k-1)       #  (2n-1)!! = 1*3*5*...*(2n-1)

sigs = [1, 1.5, 2, 2.5, 3, 3.5, 4, 5, 7, 10]   # experiment with these sigma values

print("\n   importance sampling for the " + str(tn) + "-th moment with " + str(N)
    +     " samples \n")

print("   sigma    estimate      error       error bar   relative error \n")

for sig in sigs:        #  do a calculation for each element of L

    Vsum   = 0.          #  sum of V for all paths with this l
    VsqSum = 0.          #  sum of the squares
    Zsamps = rg.normal( 0, sig, N)    # make N samples all at once
    
    a = -.5*(1.-(1/(sig*sig)))      # avoid recalculating this constant
    for Z in Zsamps:     #  loop over the samples
        
        V      = np.power( Z, tn)*np.exp(a*Z**2)  # Z^(2n) with importance weight
        Vsum   = Vsum   + V
        VsqSum = VsqSum + V**2
    
    Ahat  = Vsum/N                           # sample mean
    sVhat = np.sqrt( VsqSum/N - Ahat**2)     # estimate stardard deviation of V
    sAhat = sig*sVhat / np.sqrt(N)           # sigma factors (see week 2)
    Ahat  = sig*Ahat
    out = "{sig:7.2f}   {Ahat:10.3e}    {diff:10.3e}   {sAhat:10.3e}    {re:10.3e}"
    diff = Ahat - A
    re   = diff/A
    out = out.format(sig = sig, Ahat = Ahat, diff = diff, sAhat = sAhat, re=re)
    print(out)