Search code examples
pythonstatisticsprobabilitysamplingmontecarlo

Integral using importance sampling formula


I am trying to do an integral via importance sampling formula. The algorithm is as follows: Code

The value for the integral should give: 0.838932960013382

Also I must use the next probability distribution to generate 1,000,000 random numbers between 0 and 1. I also use the next weight function. Finally with these numbers I must calculate this formula. But I am getting the numerical value wrong, and I am not sure about the calculation for the 1.000.000 random numbers.


Solution

  • You don't compute weights normalization and inverse properly, check your math

    Code below, Python 3.9, Win 10 x64

    import numpy as np
    from scipy import integrate
    
    def f(x):
        return 1.0/((np.exp(x)+1.0)*np.sqrt(x))
    
    def w(x):
        return 0.5/np.sqrt(x)
    
    def inv(x):
        return x*x
    
    rng = np.random.default_rng()
    
    N = 100000
    
    x = rng.random(N)
    
    p = inv(x)
    
    q = f(p)/w(p)
    
    print(np.mean(q))
    
    print(integrate.quad(f, 0, 1))
    

    prints

    0.8389948486429488
    (0.8389329600133858, 2.0727863869751673e-13)
    

    looks good?