source code for Tokyo.SciPy

bv.py

from __future__ import division

import numpy as np

class _callableArray(np.ndarray):
    def __call__(self, *arg):
        return np.array([f(*arg) for f in self])

def carray(funcs):
    buf = np.array(funcs)
    return _callableArray(shape=buf.shape,
                          dtype=buf.dtype,
                          buffer=buf)

N = 25
M = 25
L = 100

SD = 0.3

sample_xs = np.linspace(0, 1, N)

def sin2pix(xs): return np.sin(2 * np.pi * xs)

def calc_lambs(min, max, per_one):
    return np.exp(np.linspace(min, max, (max - min) * per_one + 1))

lambs = calc_lambs(-2, 3, 4)

setting = {'h': sin2pix, 'sd': SD, 'sample_xs': sample_xs,
           'lambs': lambs}

def sample(seed=1, setting=setting, L=L):
    h = setting['h']
    sd = setting['sd']
    sample_xs = setting['sample_xs']
    return np.random.normal(h(sample_xs), sd, (L, len(sample_xs)))

tss = sample()

def gauss_func(mu, s):
    return lambda x: np.exp(-(x - mu) ** 2 / (2 * s ** 2))

def sigma(a): return 1 / (1 + np.exp(-a))
def sigmoid_func(mu, s):
    return lambda x: sigma((x - mu) / s)

def basis_funcs(s, func=gauss_func, M=M, offset=False):
    mus = np.arange(0, 1, 1 / (M - 1))
    if offset: mus += 0.5 / (M - 1)
    return carray([np.vectorize(lambda x: 1)] + [func(mu, s) for mu in mus])

gauss005 = basis_funcs(0.05)

def PhiT_and_PhiT_Phi(phis, xs):
    PhiT = phis(xs)
    PhiT_Phi = np.dot(PhiT, PhiT.T)
    return (PhiT, PhiT_Phi)

def make_weights_func(PhiT, PhiT_Phi, lamb):
    (dim, _) = PhiT_Phi.shape
    lambI = lamb * np.eye(dim)
    A = lambI + PhiT_Phi
    def weights_func(ts):
        b = np.dot(PhiT, ts)
        return np.linalg.solve(A, b)
    return weights_func

def wss(PhiT, PhiT_Phi, lamb, tss):
    weights_func = make_weights_func(PhiT, PhiT_Phi, lamb)
    return np.apply_along_axis(weights_func, 1, tss)

def mean_and_var_ys(wss, phis_xs):
    yss = np.dot(wss, phis_xs)
    return (np.mean(yss, axis=0), np.var(yss, axis=0))

def dataset(basis_funcs, basis_funcs_name, setting=setting, tss=tss):
    xs = setting['sample_xs']
    (PhiT, PhiT_Phi) = PhiT_and_PhiT_Phi(basis_funcs, xs)
    lambs = setting['lambs']
    phis_xs = basis_funcs(xs)

    # not good!
    wsss = [wss(PhiT, PhiT_Phi, lamb, tss) for lamb in lambs]
    mean_yss, var_yss = zip(*[mean_and_var_ys(wss_, phis_xs) for wss_ in wsss])
                                      
    return {
        'basis_funcs'     : basis_funcs,
        'basis_funcs_name': basis_funcs_name,
        'setting'         : setting,
        'tss'             : tss,
        'wsss'            : np.array(wsss),
        'mean_yss'        : np.array(mean_yss),
        'var_yss'         : np.array(var_yss)
    }

def bias_2(mean_ys, hs):
    return np.mean((hs - mean_ys) ** 2)

def variance(var_ys):
    return np.mean(var_ys)

# ----------------------------------------------------------------------

import matplotlib.pyplot as plt
import os.path

plot_xs = np.linspace(0, 1, 100+1)

def plot_setting(h, sd, phis, plot_xs=plot_xs):
    hs = h(plot_xs)
    upper_hs = hs + sd
    lower_hs = hs - sd
    phis_xs = phis(plot_xs)
    return {
        'plot_xs' : plot_xs,
        'hs'      : hs,
        'upper_hs': upper_hs,
        'lower_hs': lower_hs,
        'phis_xs' : phis_xs
    }

def plot_setting_from_dataset(dataset):
    h    = dataset['setting']['h']
    sd   = dataset['setting']['sd']
    phis = dataset['basis_funcs']
    return plot_setting(h, sd, phis)

def prepare_graph_1():
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.set_xlabel("x")
    ax.set_ylabel("t")
    ax.axis([0, 1, -1.5, 1.5])
    return (fig, ax)

def prepare_graph_2():
    fig = plt.figure()
    ax1 = fig.add_subplot(211)
    ax2 = fig.add_subplot(212)
    for ax in (ax1, ax2):
        ax.set_xlabel("x")
        ax.set_ylabel("t")
        ax.axis([0, 1, -1.5, 1.5])
    return (fig, ax1, ax2)

def add_title(ax, sample_no, basis_funcs_name, lamb):
    return ax.set_title(
        r"sample #%d, basis functions: %s, ln$\lambda=$ %f" %
        (sample_no, basis_funcs_name, np.log(lamb)))

def add_hrange(ax, xs, hs_upper, hs_lower):
    return ax.fill_between(xs, hs_upper, hs_lower,
                           color='pink', alpha=0.2)

def add_hs(ax, xs, hs): return ax.plot(xs, hs, 'g-')

def add_points(ax, xs, ts): return ax.plot(xs, ts, 'bo')

def add_w_phiss(ax, xs, w_phiss): return ax.plot(xs, w_phiss.T, 'y--')

def add_ys(ax, xs, ys, fmt='b-'): return ax.plot(xs, ys, fmt)

def standard_graph_fig(sample_no, basis_funcs_name, lamb,
                       xs, upper_hs, lower_hs, hs,
                       sample_xs, ts,
                       w_phiss, ys):
    fig, ax = prepare_graph_1()
    add_title(ax, sample_no, basis_funcs_name, lamb)
    add_hrange(ax, xs, upper_hs, lower_hs)
    add_hs(ax, xs, hs)
    add_points(ax, sample_xs, ts)
    add_w_phiss(ax, xs, w_phiss)
    add_ys(ax, xs, ys)
    return (fig, ax)

def add_yss(ax, xs, yss): return ax.plot(xs, yss.T, 'r-')

def add_errorbars(ax, xs, ys, bars):
    return ax.errorbar(xs, ys, yerr=bars, fmt='r.')
    
def bias_variance_fig(basis_funcs_name, lamb,
                      xs, hs, yss, plot_mean_ys,
                      sample_xs, mean_ys, var_ys):
    fig, ax1, ax2 = prepare_graph_2()
    fig.suptitle(r"basis functions: %s, ln$\lambda=$ %f" %
                 (basis_funcs_name, np.log(lamb)))
    add_yss(ax1, xs, yss)
    add_hs(ax2, xs, hs)
    add_ys(ax2, xs, plot_mean_ys, fmt='r-')
    add_errorbars(ax2, sample_xs, mean_ys, np.sqrt(var_ys))
    return (fig, ax1, ax2)

default_dir = "/dev/shm/graph"

def plot_bias_variance(dataset, plot_setting, dir=default_dir):
    basis_funcs_name = dataset['basis_funcs_name']
    setting          = dataset['setting']
    wsss             = dataset['wsss']
    mean_yss         = dataset['mean_yss']
    var_yss          = dataset['var_yss']
    
    sample_xs = setting['sample_xs']
    lambs     = setting['lambs']

    xs       = plot_setting['plot_xs']
    hs       = plot_setting['hs']
    phis_xs  = plot_setting['phis_xs']

    wsss20 = wsss[:, :20, :]
    mean_wss = np.mean(wsss, axis=1)
    
    for i in xrange(np.alen(lambs)):
        lamb = lambs[i]
        mean_ys = mean_yss[i, :]
        var_ys = var_yss[i, :]

        wss = wsss20[i, :, :]
        yss = np.dot(wss, phis_xs)

        mean_ws = mean_wss[i, :]
        plot_mean_ys = np.dot(phis_xs.T, mean_ws)

        bias_variance_fig(basis_funcs_name, lamb,
                          xs, hs, yss, plot_mean_ys,
                          sample_xs, mean_ys, var_ys)
        plt.savefig(os.path.join(dir, "lambda%.5f.png" % lamb))
    plt.close('all')

def plot_all_sample(dataset, i, plot_setting, dir=default_dir):
    basis_funcs_name = dataset['basis_funcs_name']
    setting          = dataset['setting']
    tss              = dataset['tss']
    wss              = dataset['wsss'][i, :, :]
    
    sample_xs = setting['sample_xs']
    lamb      = setting['lambs'][i]

    xs       = plot_setting['plot_xs']
    hs       = plot_setting['hs']
    upper_hs = plot_setting['upper_hs']
    lower_hs = plot_setting['lower_hs']
    phis_xs  = plot_setting['phis_xs']

    w_phisss = wss[:, :, np.newaxis] * phis_xs[np.newaxis, :, :]
    yss = np.sum(w_phisss, axis=1)

    for sample_no in xrange(np.alen(tss)):
        ts = tss[sample_no, :]
        w_phiss = w_phisss[sample_no, :, :]
        ys = yss[sample_no, :]
        standard_graph_fig(sample_no, basis_funcs_name, lamb,
                           xs, upper_hs, lower_hs, hs,
                           sample_xs, ts,
                           w_phiss, ys)
        plt.savefig(os.path.join(dir, "sample%03d.png" % sample_no))
    plt.close('all')


def plot_through_one_sample(dataset, sample_no, plot_setting,
                            dir=default_dir):
    basis_funcs_name = dataset['basis_funcs_name']
    setting          = dataset['setting']
    ts               = dataset['tss'][sample_no, :]
    wsss             = dataset['wsss']
    
    sample_xs = setting['sample_xs']
    lambs     = setting['lambs']

    xs       = plot_setting['plot_xs']
    hs       = plot_setting['hs']
    upper_hs = plot_setting['upper_hs']
    lower_hs = plot_setting['lower_hs']
    phis_xs  = plot_setting['phis_xs']

    for i in xrange(np.alen(lambs)):
        lamb = lambs[i]
        ws = wsss[i, sample_no, :]
        w_phiss = ws[:, np.newaxis] * phis_xs
        ys = np.sum(w_phiss, axis=0)
        
        standard_graph_fig(sample_no, basis_funcs_name, lamb,
                           xs, upper_hs, lower_hs, hs,
                           sample_xs, ts,
                           w_phiss, ys)
        plt.savefig(os.path.join(dir, "lambda%.5f.png" % lamb))
    plt.close('all')
    
def plot_testerror(dataset, plot_setting, dir=default_dir):
    basis_funcs_name = dataset['basis_funcs_name']
    lnlambs = np.log(dataset['setting']['lambs'])
    mean_yss = dataset['mean_yss']
    var_yss = dataset['var_yss']
    hs = dataset['setting']['h'](dataset['setting']['sample_xs'])

    bias_2s = np.apply_along_axis(bias_2, 1, mean_yss, hs)
    variances = np.apply_along_axis(variance, 1, var_yss)

    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.set_title("Test error, basis functions: %s" % basis_funcs_name)
    ax.set_xlabel(r"$ln\lambda$")
    ax.plot(lnlambs, bias_2s  , 'b-', label=r"$(bias)^2$")
    ax.plot(lnlambs, variances, 'r-', label=r"$variance$")
    ax.plot(lnlambs, bias_2s + variances, 'm-',
            label="$(bias)^2 + variance$")
    ax.legend()

    plt.savefig(os.path.join(dir, "test_error.png"))
    plt.close()

# example ipython session

# from __future__ import division
# import numpy as np
# import bv

# ds_gauss005 = bv.dataset(bv.gauss005, "Gaussian(s=0.05)")
# pset = bv.plot_setting_from_dataset(ds_gauss005)

# bv.plot_bias_variance(ds_gauss005, pset)
# bv.plot_all_sample(ds_gauss005, 10, pset)
# bv.plot_through_one_sample(ds_gauss005, 50, pset)
# bv.plot_testerror(ds_gauss005, pset)