The Kumaraswamy Distribution

Reference the Kumaraswamy Wikipedia Page.

Imports

from qmcpy import *
from scipy.special import gamma
from numpy import *
import matplotlib
from matplotlib import pyplot
%matplotlib inline
pyplot.rc('font', size=16)          # controls default text sizes
pyplot.rc('axes', titlesize=16)     # fontsize of the axes title
pyplot.rc('axes', labelsize=16)     # fontsize of the x and y labels
pyplot.rc('xtick', labelsize=16)    # fontsize of the tick labels
pyplot.rc('ytick', labelsize=16)    # fontsize of the tick labels
pyplot.rc('legend', fontsize=16)    # legend fontsize
pyplot.rc('figure', titlesize=16)   # fontsize of the figure title

2D Density Plot

def plt_2d_kumaraswamy(kamaraswamy,n):
    print(kamaraswamy)
    n_mesh = 502
    # Points
    t = kamaraswamy.gen_samples(n)
    # PDF
    mesh = zeros(((n_mesh)**2,3),dtype=float)
    grid_tics = linspace(0,1,n_mesh)
    x_mesh,y_mesh = meshgrid(grid_tics,grid_tics)
    mesh[:,0] = x_mesh.flatten()
    mesh[:,1] = y_mesh.flatten()
    mesh[:,2] = kamaraswamy._weight(mesh[:,:2])
    z_mesh = mesh[:,2].reshape((n_mesh,n_mesh))
    # plots
    fig,ax = pyplot.subplots(figsize=(7,7),nrows=1,ncols=1)
    #   colors
    clevel = arange(mesh[:,2].min(),mesh[:,2].max(),.025)
    cmap = matplotlib.colors.LinearSegmentedColormap.from_list("", [(.95,.95,.95),(0,0,1)]) #cmap = pyplot.get_cmap('Blues')
    #   contour + scatter plot
    ax.contourf(x_mesh,y_mesh,z_mesh,clevel,cmap=cmap,extend='both')
    ax.scatter(t[:,0],t[:,1],s=5,color='w')
    #   axis
    for nsew in ['top','bottom','left','right']: ax.spines[nsew].set_visible(False)
    ax.xaxis.set_ticks_position('none')
    ax.yaxis.set_ticks_position('none')
    ax.set_aspect(1)
    ax.set_xlim([0,1])
    ax.set_xticks([0,1])
    ax.set_ylim([0,1])
    ax.set_yticks([0,1])
    #   labels
    ax.set_xlabel('$T_1$')
    ax.set_ylabel('$T_2$')
    ax.set_title('Kamaraswamy PDF and Random Samples')
    #   metas
    fig.tight_layout()
kumaraswamy = Kumaraswamy(Sobol(2),a=[2,3],b=[3,5]) # won't plot correctly with multiple composed transforms
plt_2d_kumaraswamy(kumaraswamy,n=2**8)
Kumaraswamy (TrueMeasure Object)
    a               [2 3]
    b               [3 5]
../_images/kumaraswamy_6_1.png

1D Kumaraswamy Expected Value

a,b = 2,6
abs_tol = 1e-6
cf = CustomFun(true_measure=Kumaraswamy(Sobol(1),a,b), g=lambda x:x)
sol,data = CubQMCSobolG(cf,abs_tol=abs_tol).integrate()
print(data)
true_val = b*gamma(1+1/a)*gamma(b)/gamma(1+1/a+b)
error = abs(true_val-sol)
if error>abs_tol: raise Exception("QMC error %.3f not within tolerance."%error)
Solution: 0.3410
CustomFun (Integrand Object)
Sobol (DiscreteDistribution Object)
    d               1
    randomize       1
    graycode        0
    seed            96725
    mimics          StdUniform
    dim0            0
Kumaraswamy (TrueMeasure Object)
    a               2^(1)
    b               6
CubQMCSobolG (StoppingCriterion Object)
    abs_tol         1.00e-06
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
LDTransformData (AccumulateData Object)
    n_total         2^(17)
    solution        0.341
    error_bound     4.69e-07
    time_integrate  0.033

Importance Samling with a Single Kumaraswamy

# compose with a Kumaraswamy transformation
abs_tol = 1e-4
cf = CustomFun(Uniform(Kumaraswamy(Sobol(1,seed=7),a=.5,b=.5)), g=lambda x:x)
sol,data = CubQMCSobolG(cf,abs_tol=abs_tol).integrate()
print(data)
true_val = .5 # expected value of standard uniform
error = abs(true_val-sol)
if error>abs_tol: raise Exception("QMC error %.3f not within tolerance."%error)
Solution: 0.5000
CustomFun (Integrand Object)
Sobol (DiscreteDistribution Object)
    d               1
    randomize       1
    graycode        0
    seed            61616
    mimics          StdUniform
    dim0            0
Uniform (TrueMeasure Object)
    lower_bound     0
    upper_bound     1
    transform       Kumaraswamy (TrueMeasure Object)
                       a               2^(-1)
                       b               2^(-1)
CubQMCSobolG (StoppingCriterion Object)
    abs_tol         1.00e-04
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
LDTransformData (AccumulateData Object)
    n_total         2^(11)
    solution        0.500
    error_bound     1.40e-05
    time_integrate  0.035
# plot the above functions
x = linspace(0,1,1000).reshape(-1,1)
x = x[1:-1] # plotting locations
# plot
fig,ax = pyplot.subplots(ncols=2,figsize=(15,5))
#    density
rho = cf.true_measure.transform._weight(x)
tfvals = cf.true_measure.transform._transform_r(x)
ax[0].hist(tfvals,density=True,bins=50,color='c')
ax[0].plot(x,rho,color='g',linewidth=2)
ax[0].set_title('Sampling density')
#    functions
gs = cf.g(x)
fs = cf.f(x)
ax[1].plot(x,gs,color='r',label='g')
ax[1].plot(x,fs,color='b',label='f')
ax[1].legend()
ax[1].set_title('functions');
# metas
for i in range(2):
    ax[i].set_xlim([0,1])
    ax[i].set_xticks([0,1])
../_images/kumaraswamy_11_0.png

Importance Sampling with 2 (Composed) Kumaraswamys

# compose multiple Kumaraswamy distributions
abs_tol = 1e-3
dd = Sobol(1,seed=7)
t1 = Kumaraswamy(dd,a=.5,b=.5) # first transformation
t2 = Kumaraswamy(t1,a=5,b=2) # second transformation
cf = CustomFun(Uniform(t2), g=lambda x:x)
sol,data = CubQMCSobolG(cf,abs_tol=abs_tol).integrate() # expected value of standard uniform
print(data)
true_val = .5
error = abs(true_val-sol)
if error>abs_tol: raise Exception("QMC error %.3f not within tolerance."%error)
Solution: 0.4998
CustomFun (Integrand Object)
Sobol (DiscreteDistribution Object)
    d               1
    randomize       1
    graycode        0
    seed            61616
    mimics          StdUniform
    dim0            0
Uniform (TrueMeasure Object)
    lower_bound     0
    upper_bound     1
    transform       Kumaraswamy (TrueMeasure Object)
                       a               5
                       b               2^(1)
                       transform       Kumaraswamy (TrueMeasure Object)
                                          a               2^(-1)
                                          b               2^(-1)
CubQMCSobolG (StoppingCriterion Object)
    abs_tol         0.001
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
LDTransformData (AccumulateData Object)
    n_total         2^(11)
    solution        0.500
    error_bound     5.19e-04
    time_integrate  0.002
x = linspace(0,1,1000).reshape(-1,1)
x = x[1:-1] # plotting locations
# plot
fig,ax = pyplot.subplots(ncols=2,figsize=(15,5))
#    density
tfvals = cf.true_measure.transform._transform_r(x)
ax[0].hist(tfvals,density=True,bins=50,color='c')
ax[0].set_title('Sampling density')
#    functions
gs = cf.g(x)
fs = cf.f(x)
ax[1].plot(x,gs,color='r',label='g')
ax[1].plot(x,fs,color='b',label='f')
ax[1].legend()
for i in range(2):
    ax[i].set_xlim([0,1])
    ax[i].set_xticks([0,1])
ax[1].set_title('functions');
../_images/kumaraswamy_14_0.png

Can we Improve the Keister function?

abs_tol = 1e-4
d = 1
# standard method
keister_std = Keister(Sobol(d))
sol_std,data_std = CubQMCSobolG(keister_std,abs_tol=abs_tol).integrate()
data_std
Solution: 1.3804
Keister (Integrand Object)
Sobol (DiscreteDistribution Object)
    d               1
    randomize       1
    graycode        0
    seed            57958
    mimics          StdUniform
    dim0            0
Lebesgue (TrueMeasure Object)
    transform       Gaussian (TrueMeasure Object)
                       mean            0
                       covariance      2^(-1)
                       decomp_type     pca
CubQMCSobolG (StoppingCriterion Object)
    abs_tol         1.00e-04
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
LDTransformData (AccumulateData Object)
    n_total         2^(14)
    solution        1.380
    error_bound     1.58e-05
    time_integrate  0.011
# Kumaraswamy importance sampling
dd = Sobol(d,seed=7)
t1 = Kumaraswamy(dd,a=.8,b=.8) # first transformation
keister_kuma = Keister(Gaussian(t1,covariance=1/2))
sol_kuma,data_kuma = CubQMCSobolG(keister_kuma,abs_tol=abs_tol).integrate()
print(data_kuma)
Solution: 1.3804
Keister (Integrand Object)
Sobol (DiscreteDistribution Object)
    d               1
    randomize       1
    graycode        0
    seed            61616
    mimics          StdUniform
    dim0            0
Lebesgue (TrueMeasure Object)
    transform       Gaussian (TrueMeasure Object)
                       mean            0
                       covariance      2^(-1)
                       decomp_type     pca
                       transform       Kumaraswamy (TrueMeasure Object)
                                          a               0.800
                                          b               0.800
CubQMCSobolG (StoppingCriterion Object)
    abs_tol         1.00e-04
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
LDTransformData (AccumulateData Object)
    n_total         2^(11)
    solution        1.380
    error_bound     9.59e-05
    time_integrate  0.004
t_frac = data_kuma.time_integrate/data_std.time_integrate
n_frac = data_kuma.n_total/data_std.n_total
print('Kuma importance sampling took %.1f%% of the time and %.1f%% of the samples compared to default keister.'%(t_frac*100,n_frac*100))
Kuma importance sampling took 39.5% of the time and 12.5% of the samples compared to default keister.
x = linspace(0,1,1000).reshape(-1,1)
x = x[1:-1] # plotting locations
# plot
fig,ax = pyplot.subplots(figsize=(10,8))
#    functions
fs = [keister_std.f,keister_kuma.f]
labels = ['Default Keister','Kuma Keister']
colors = ['m','c']
for f,label,color in zip(fs,labels,colors): ax.plot(x,f(x),color=color,label=label)
ax.legend()
ax.set_xlim([0,1])
ax.set_xticks([0,1])
ax.set_title('functions');
../_images/kumaraswamy_20_0.png