quasimontecarlo.py

import ghalton
from pylab import *
from scipy.stats import gamma


def qmc_gamma(alpha, N):
    M = len(alpha)
    sequencer = ghalton.Halton(2*M)
    valid = ones(shape=(N,), dtype=bool)
    hs = array(sequencer.get(N)).transpose()
    X = zeros(shape=(M,N), dtype=float)
    for i in range(M):
        a = pow(2*alpha[i]-1, -0.5)
        b = alpha[i] - log(4)
        c = alpha[i] + pow(a, -1.0)

        Y = a*log(hs[2*i]/(1-hs[2*i]))
        X[i] = alpha[i]*exp(Y)
        Z = hs[2*i]*hs[2*i]*hs[2*i+1]
        R = b + c*Y - X[i]
        cond1 = R + 2.5040774 - 4.5*Z > 0
        cond2 = R > log(Z)
        valid = logical_and(valid, logical_or(cond1, cond2))
    return X[:, valid]

def mc_gamma(alpha, N):
    return array([ gamma(a,1).rvs(N) for a in alpha])

def mc(x):
    return mean( maximum(maximum(x[0]-x[1], x[0]-x[2]), 0) )

def f(x):
    return mean( maximum(maximum(x[0]-x[1], x[0]-x[2]), 0) > 0.1 )

result = []
Nvals = arange(10000, 500000, 5000)
for N in Nvals:
    #result.append( [ mc(qmc_gamma([3,4,5],N)), mc(mc_gamma([3,4,5], N)) ] )
    result.append( [ f(qmc_gamma([3,4,5],N)), f(mc_gamma([3,4,5], N)) ] )

result = array(result)

plot(Nvals, result[:,0], 'bo', label="Quasi-Monte Carlo")
plot(Nvals, result[:,1], 'ro', label="Monte carlo")
xlabel("Number of samples")
ylabel("Estimated value")
legend()
show()