Visual and analytic solution comparison
=======================================

The idea here is to compare the various analytic functions with
approximate solutions, to verify that the analytic implementation is
correct. Ideally, we want to see two things:

1. As the number of approximation points increases, the error should go
   down and tend towards zero.

2. With a sufficient number of approximation points, the analytic and
   approximate solutions should appear visually identical.

.. code:: python

    # imports 
    %matplotlib inline
    from bayesian_quadrature import BQ, bq_c, gauss_c, util
    from gp import GaussianKernel
    
    # set loglevel to debug
    import logging
    logger = logging.getLogger("bayesian_quadrature")
    logger.setLevel("DEBUG")
.. code:: python

    # seed the numpy random generator, so we always get the same randomness
    np.random.seed(8728)
.. code:: python

    # parameters
    options = {
        'n_candidate': 10,
        'x_mean': 0.0,
        'x_var': 10.0,
        'candidate_thresh': 0.5,
        'kernel': GaussianKernel,
        'optim_method': 'L-BFGS-B'
    }
    
    # different number of approximation points to try, to see if the
    # approximation converges, and the range over which we will compute
    # the approximation points
    N = np.logspace(np.log10(30), np.log10(500), 100)
    xmin, xmax = -10, 10
Here we define some helper functions to plot the error/visual
comparisons, and to print information about the worst error and
approximate/analytic solution values.

.. code:: python

    def plot_0d(calc, approx):
        fig, ax = plt.subplots()
        ax.plot(N, np.abs(calc - approx), lw=2, color='r')
        ax.set_title("Error")
        ax.set_xlabel("# approximation points")
        ax.set_xlim(N.min(), N.max())
        util.set_scientific(ax, -5, 4)
.. code:: python

    def print_0d(calc, approx):
        err = np.abs(calc - approx[-1])
        print "error : %s" % err
        print "approx: %s" % approx[-1]
        print "calc  : %s" % calc
.. code:: python

    def plot_1d(calc, approx):
        fig, (ax1, ax2) = plt.subplots(1, 2)
    
        ax1.plot(N, np.abs(calc[:, None] - approx).sum(axis=0), lw=2, color='r')
        ax1.set_title("Error")
        ax1.set_xlabel("# approximation points")
        ax1.set_xlim(N.min(), N.max())
        util.set_scientific(ax1, -5, 4)
    
        ax2.plot(approx[:, -1], lw=2, label="approx")
        ax2.plot(calc, '--', lw=2, label="calc")
        ax2.set_title("Comparison")
        ax2.legend(loc=0)
        util.set_scientific(ax2, -5, 4)
    
        fig.set_figwidth(8)
        plt.tight_layout()
.. code:: python

    def print_1d(calc, approx):
        maxerr = np.abs(calc - approx[..., -1]).max()
        approx_sum, calc_sum = approx[..., -1].sum(), calc.sum()
        print "max error  : %s" % maxerr
        print "sum(approx): %s" % approx_sum
        print "sum(calc)  : %s" % calc_sum
.. code:: python

    def plot_2d(calc, approx):
        fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4)
    
        ax1.plot(N, np.abs(calc[:, :, None] - approx).sum(axis=0).sum(axis=0), lw=2, color='r')
        ax1.set_title("Error")
        ax1.set_xlabel("# approximation points")
        ax1.set_xlim(N.min(), N.max())
        util.set_scientific(ax1, -5, 4)
    
        vmax = max(approx[:, :, -1].max(), calc.max())
    
        ax2.imshow(approx[:, :, -1], cmap='gray', interpolation='nearest', vmin=0, vmax=vmax)
        ax2.set_title("Approximation")
    
        ax3.imshow(calc, cmap='gray', interpolation='nearest', vmin=0, vmax=vmax)
        ax3.set_title("Calculated")
    
        i = 0
        ax4.plot(approx[i, :, -1], lw=2, label='approx')
        ax4.plot(calc[i], '--', lw=2, label='calc')
        ax4.set_title("Comparison at index %d" % i)
        ax4.legend(loc=0)
        util.set_scientific(ax4, -5, 4)
    
        fig.set_figwidth(14)
        plt.tight_layout()
Now we create the Bayesian Quadrature object that we'll be testing
against.

.. code:: python

    # choose points and evaluate them under a standard normal distribution
    x = np.linspace(-5, 5, 9)
    f_y = lambda x: scipy.stats.norm.pdf(x, 0, 1)
    y = f_y(x)
    
    # create the bayesian quadrature object
    bq = BQ(x, y, **options)
    bq.init(
        params_tl=(15, 2, 0.), 
        params_l=(0.2, 1.3, 0.))
    
    # plot the BQ approximation
    fig, axes = bq.plot(f_l=f_y, xmin=xmin, xmax=xmax)

.. parsed-literal::

    DEBUG:bayesian_quadrature:Choosing candidate points



.. image:: visual-tests_files/visual-tests_12_1.png


.. code:: python

    fig, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True)
    xmin, xmax = -10, 10
    bq.plot_l(ax1, f_l=f_y, xmin=xmin, xmax=xmax)
    bq.plot_expected_squared_mean(ax2, xmin=xmin, xmax=xmax)
    bq.plot_expected_variance(ax3, xmin=xmin, xmax=xmax)
    fig.set_figheight(8)
    plt.tight_layout()


.. image:: visual-tests_files/visual-tests_13_0.png


.. code:: python

    x_l = np.array(bq.gp_l.x[None], order='F')
    h_l = bq.gp_l.K.h
    w_l = np.array([bq.gp_l.K.w], order='F')
    
    x_tl = np.array(bq.gp_log_l.x[None], order='F')
    h_tl = bq.gp_log_l.K.h
    w_tl = np.array([bq.gp_log_l.K.w], order='F')
    
    xm = bq.options['x_mean']
    xC = bq.options['x_cov']
int\_K
------

Test
:math:`\int K(x^\prime, x) \mathcal{N}(x^\prime | \mu, \Sigma)\ \mathrm{d}x^\prime`

.. code:: python

    calc = np.empty(bq.nsc, order='F')
    gauss_c.int_K(calc, x_l, h_l, w_l, xm, xC)
    
    approx = np.empty((bq.nsc, N.size), order='F')
    for i, n in enumerate(N):
        xo = np.linspace(xmin, xmax, n)
        Kxxo = np.array(bq.gp_l.Kxxo(xo), order='F')
        gauss_c.approx_int_K(
            approx[:, i], 
            np.array(xo[None], order='F'),
            Kxxo, xm, xC)
.. code:: python

    plot_1d(calc, approx)
    print_1d(calc, approx)

.. parsed-literal::

    max error  : 1.55678221442e-09
    sum(approx): 0.0339891854186
    sum(calc)  : 0.033989188741



.. image:: visual-tests_files/visual-tests_17_1.png


int\_K1\_K2
-----------

Test
:math:`\int K_1(x_1, x^\prime) K_2(x^\prime, x_2) \mathcal{N}(x^\prime | \mu, \Sigma)\ \mathrm{d}x^\prime`

.. code:: python

    calc = np.empty((bq.nsc, bq.ns), order='F')
    gauss_c.int_K1_K2(calc, x_l, x_tl, h_l, w_l, h_tl, w_tl, xm, xC)
    
    approx = np.empty((bq.nsc, bq.ns, N.size), order='F')
    for i, n in enumerate(N):
        xo = np.linspace(xmin, xmax, n)
        K1xxo = np.array(bq.gp_l.Kxxo(xo), order='F')
        K2xxo = np.array(bq.gp_log_l.Kxxo(xo), order='F')
        gauss_c.approx_int_K1_K2(
            approx[:, :, i], 
            np.array(xo[None], order='F'), 
            K1xxo, K2xxo, xm, xC)
.. code:: python

    plot_2d(calc, approx)
    print_1d(calc, approx)

.. parsed-literal::

    max error  : 0.000933347326451
    sum(approx): 5.55089408134
    sum(calc)  : 5.55090592396



.. image:: visual-tests_files/visual-tests_20_1.png


int\_int\_K1\_K2\_K1
--------------------

Test
:math:`\int \int K_1(x, x_1^\prime) K_2(x_1^\prime, x_2^\prime) K_1(x_2^\prime, x) \mathcal{N}(x_1^\prime | \mu, \Sigma) \mathcal{N}(x_2^\prime | \mu, \Sigma)\ \mathrm{d}x_1^\prime\ \mathrm{d}x_2^\prime`

.. code:: python

    calc = np.empty((bq.nsc, bq.nsc), order='F')
    gauss_c.int_int_K1_K2_K1(calc, x_l, h_l, w_l, h_tl, w_tl, xm, xC)
        
    approx = np.empty((bq.nsc, bq.nsc, N.size), order='F')
    for i, n in enumerate(N):
        xo = np.linspace(xmin, xmax, n)
        K1xxo = np.array(bq.gp_l.Kxxo(xo), order='F')
        K2xoxo = np.array(bq.gp_log_l.Kxoxo(xo), order='F')
        gauss_c.approx_int_int_K1_K2_K1(
            approx[:, :, i], 
            np.array(xo[None], order='F'),
            K1xxo, K2xoxo, xm, xC)
.. code:: python

    plot_2d(calc, approx)
    print_1d(calc, approx)

.. parsed-literal::

    max error  : 7.56401320431e-12
    sum(approx): 0.0230979040561
    sum(calc)  : 0.0230979040904



.. image:: visual-tests_files/visual-tests_23_1.png


int\_int\_K1\_K2
----------------

Test
:math:`\int \int K_1(x_2^\prime, x_1^\prime) K_2(x_1^\prime, x) \mathcal{N}(x_1^\prime | \mu, \Sigma) \mathcal{N}(x_2^\prime | \mu, \Sigma)\ \mathrm{d}x_1^\prime\ \mathrm{d}x_2^\prime`

.. code:: python

    calc = np.empty(bq.ns, order='F')
    gauss_c.int_int_K1_K2(calc, x_tl, h_l, w_l, h_tl, w_tl, xm, xC)
    
    approx = np.empty((bq.ns, N.size), order='F')
    for i, n in enumerate(N):
        xo = np.linspace(xmin, xmax, n)
        K1xoxo = np.array(bq.gp_l.Kxoxo(xo), order='F')
        K2xxo = np.array(bq.gp_log_l.Kxxo(xo), order='F')
        gauss_c.approx_int_int_K1_K2(
            approx[:, i],
            np.array(xo[None], order='F'),
            K1xoxo, K2xxo, xm, xC)
.. code:: python

    plot_1d(calc, approx)
    print_1d(calc, approx)

.. parsed-literal::

    max error  : 2.90292116595e-07
    sum(approx): 0.576959029431
    sum(calc)  : 0.576959869272



.. image:: visual-tests_files/visual-tests_26_1.png


int\_int\_K
-----------

Test
:math:`\int \int K(x_1^\prime, x_2^\prime) \mathcal{N}(x_1^\prime | \mu, \Sigma) \mathcal{N}(x_2^\prime | \mu, \Sigma)\ \mathrm{d}x_1^\prime\ \mathrm{d}x_2^\prime`

.. code:: python

    calc = gauss_c.int_int_K(1, h_l, w_l, xm, xC)
    
    approx = np.empty(N.size, order='F')
    for i, n in enumerate(N):
        xo = np.linspace(xmin, xmax, n)
        Kxoxo = np.array(bq.gp_l.Kxoxo(xo), order='F')
        approx[i] = gauss_c.approx_int_int_K(
            np.array(xo[None], order='F'), Kxoxo, xm, xC)
.. code:: python

    plot_0d(calc, approx)
    print_0d(calc, approx)

.. parsed-literal::

    error : 1.01454767464e-07
    approx: 0.00342631605819
    calc  : 0.00342641751296



.. image:: visual-tests_files/visual-tests_29_1.png


Mean of Z
---------

Test :math:`E[Z]`.

.. code:: python

    calc = bq._exact_Z_mean()
    approx = np.empty(N.size)
    for i, n in enumerate(N):
        xo = np.linspace(xmin, xmax, n)
        approx[i] = bq._approx_Z_mean(xo)
.. code:: python

    plot_0d(calc, approx)
    print_0d(calc, approx)

.. parsed-literal::

    error : 1.8803056015e-08
    approx: 0.119771024599
    calc  : 0.119771005796



.. image:: visual-tests_files/visual-tests_32_1.png


Variance of Z
-------------

Test :math:`\mathrm{Var}(Z)`.

.. code:: python

    calc = bq._exact_Z_var()
    approx = np.empty(N.size)
    for i, n in enumerate(N):
        xo = np.linspace(xmin, xmax, n)
        approx[i] = bq._approx_Z_var(xo)
.. code:: python

    plot_0d(calc, approx)
    print_0d(calc, approx)

.. parsed-literal::

    error : 6.18632061484e-11
    approx: 5.97977452085e-07
    calc  : 5.98039315292e-07



.. image:: visual-tests_files/visual-tests_35_1.png