Example of regressing with a student t likelihood using Laplace

def student_t_approx(optimize=True, plot=True):
    """
    Example of regressing with a student t likelihood using Laplace
    """
    real_std = 0.1
    #Start a function, any function
    X = np.linspace(0.0, np.pi*2, 100)[:, None]
    Y = np.sin(X) + np.random.randn(*X.shape)*real_std
    Y = Y/Y.max()
    Yc = Y.copy()

    X_full = np.linspace(0.0, np.pi*2, 500)[:, None]
    Y_full = np.sin(X_full)
    Y_full = Y_full/Y_full.max()

    #Slightly noisy data
    Yc[75:80] += 1

    #Very noisy data
    #Yc[10] += 100
    #Yc[25] += 10
    #Yc[23] += 10
    #Yc[26] += 1000
    #Yc[24] += 10
    #Yc = Yc/Yc.max()

    #Add student t random noise to datapoints
    deg_free = 1
    print("Real noise: ", real_std)
    initial_var_guess = 0.5
    edited_real_sd = initial_var_guess

    # Kernel object
    kernel1 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
    kernel2 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
    kernel3 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
    kernel4 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])

    #Gaussian GP model on clean data
    m1 = GPy.models.GPRegression(X, Y.copy(), kernel=kernel1)
    # optimize
    m1['.*white'].constrain_fixed(1e-5)
    m1.randomize()

    #Gaussian GP model on corrupt data
    m2 = GPy.models.GPRegression(X, Yc.copy(), kernel=kernel2)
    m2['.*white'].constrain_fixed(1e-5)
    m2.randomize()

    #Student t GP model on clean data
    t_distribution = GPy.likelihoods.StudentT(deg_free=deg_free, sigma2=edited_real_sd)
    laplace_inf = GPy.inference.latent_function_inference.Laplace()
    m3 = GPy.core.GP(X, Y.copy(), kernel3, likelihood=t_distribution, inference_method=laplace_inf)
    m3['.*t_scale2'].constrain_bounded(1e-6, 10.)
    m3['.*white'].constrain_fixed(1e-5)
    m3.randomize()

    #Student t GP model on corrupt data
    t_distribution = GPy.likelihoods.StudentT(deg_free=deg_free, sigma2=edited_real_sd)
    laplace_inf = GPy.inference.latent_function_inference.Laplace()
    m4 = GPy.core.GP(X, Yc.copy(), kernel4, likelihood=t_distribution, inference_method=laplace_inf)
    m4['.*t_scale2'].constrain_bounded(1e-6, 10.)
    m4['.*white'].constrain_fixed(1e-5)
    m4.randomize()
    print(m4)
    debug=True
    if debug:
        m4.optimize(messages=1)
        from matplotlib import pyplot as pb
        pb.plot(m4.X, m4.inference_method.f_hat)
        pb.plot(m4.X, m4.Y, 'rx')
        m4.plot()
        print(m4)
        return m4

    if optimize:
        optimizer='scg'
        print("Clean Gaussian")
        m1.optimize(optimizer, messages=1)
        print("Corrupt Gaussian")
        m2.optimize(optimizer, messages=1)
        print("Clean student t")
        m3.optimize(optimizer, messages=1)
        print("Corrupt student t")
        m4.optimize(optimizer, messages=1)

    if plot:
        plt.figure(1)
        plt.suptitle('Gaussian likelihood')
        ax = plt.subplot(211)
        m1.plot(ax=ax)
        plt.plot(X_full, Y_full)
        plt.ylim(-1.5, 1.5)
        plt.title('Gaussian clean')

        ax = plt.subplot(212)
        m2.plot(ax=ax)
        plt.plot(X_full, Y_full)
        plt.ylim(-1.5, 1.5)
        plt.title('Gaussian corrupt')

        plt.figure(2)
        plt.suptitle('Student-t likelihood')
        ax = plt.subplot(211)
        m3.plot(ax=ax)
        plt.plot(X_full, Y_full)
        plt.ylim(-1.5, 1.5)
        plt.title('Student-t rasm clean')

        ax = plt.subplot(212)
        m4.plot(ax=ax)
        plt.plot(X_full, Y_full)
        plt.ylim(-1.5, 1.5)
        plt.title('Student-t rasm corrupt')

    return m1, m2, m3, m4