Importance sampling
system:sage

{{{id=3|
def f(x):
    if x < 0: return 0
    if x > 2: return 0
    return 0.5*x

def f_exact_integral():
    return 1.0

def g(x):
    return x**2 * exp(-x**2)

def g_exact_integral():
    return float(sqrt(pi)/4)
    
def uniform_estimate(f, x_min, x_max, nr_samples = 1000):
    assert x_min < x_max

    s = 0

    for _ in range(nr_samples):
        x = x_max*random()    
        s += f(x)
    
    return 1.0*(x_max - x_min)*s/nr_samples
   
def importance_sampling_estimate(f, x_min, x_max, sample_width, nr_samples = 1000):
    assert x_min < x_max

    sample = [f(x + 0.5*sample_width) for x in srange(x_min, x_max, step = sample_width)]

    sample_sum = sum(sample)
    sample_normalised = [x/sample_sum for x in sample]

    csum = [0]
    csum[0] = sample[0]
    for i in range(1, len(sample)):
        csum.append(csum[i-1] + sample[i])
    
    csum_norm = [x/csum[-1] for x in csum]
   
    s = 0
    for _ in range(nr_samples):
        x = random()
        i = 0
        while csum_norm[i] < x:
            i += 1
        
        actual_point = sample_width*(i + random())

        assert actual_point >= x_min
        assert actual_point <= x_max

        s += sample_width*f(actual_point)/sample_normalised[i]
    
    return s/nr_samples

def stddev(X):
    mu = sum(X)/(1.0*len(X))
    return sqrt(sum([(x - mu)^2 for x in X])/(1.0*len(X)))
///
}}}

{{{id=10|
plot(f, (-1, 3))
///

<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{id=12|
plot(g, (0, 10))
///

<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{id=11|
N = 10
nr_samples = 1000

f_uniforms = [uniform_estimate(f, 0, 2, nr_samples) for _ in range(N)]
f_isamples = [importance_sampling_estimate(f, 0, 2, 0.1, nr_samples) for _ in range(N)]

f_uniforms_err = [abs(x - f_exact_integral()) for x in f_uniforms]
f_isamples_err = [abs(x - f_exact_integral()) for x in f_isamples]

g_uniforms = [uniform_estimate(g, 0, 8, nr_samples) for _ in range(N)]
g_isamples = [importance_sampling_estimate(g, 0, 8, 0.01, nr_samples) for _ in range(N)]

g_uniforms_err = [abs(x - g_exact_integral()) for x in g_uniforms]
g_isamples_err = [abs(x - g_exact_integral()) for x in g_isamples]
///
}}}

{{{id=13|
print stddev(f_uniforms)
print stddev(f_isamples)
///

0.00810877768691104
0.00150839328124713
}}}

{{{id=14|
P_f_uniform_errs = list_plot(f_uniforms_err, rgbcolor=(1,0,0))
P_f_isamples_errs = list_plot(f_isamples_err, rgbcolor=(0,1,0))
plot(P_f_uniform_errs + P_f_isamples_errs)
///

<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{id=16|
P_g_uniform_errs = list_plot(g_uniforms_err, rgbcolor=(1,0,0))
P_g_isamples_errs = list_plot(g_isamples_err, rgbcolor=(0,1,0))
plot(P_g_uniform_errs + P_g_isamples_errs)
///

<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{id=17|
print stddev(g_uniforms)
print stddev(g_isamples)
///

0.0193302936692150
0.000122180625633688
}}}

{{{id=19|

///
}}}

{{{id=20|

///
}}}