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
def student_t_approx ( optimize = True , plot = True ) : real_std = 0.1 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 ( ) Yc [ 75 : 80 ] += 1 deg_free = 1 print ( "Real noise: " , real_std ) initial_var_guess = 0.5 edited_real_sd = initial_var_guess 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 ] ) m1 = GPy . models . GPRegression ( X , Y . copy ( ) , kernel = kernel1 ) m1 [ '.*white' ] . constrain_fixed ( 1e-5 ) m1 . randomize ( ) m2 = GPy . models . GPRegression ( X , Yc . copy ( ) , kernel = kernel2 ) m2 [ '.*white' ] . constrain_fixed ( 1e-5 ) m2 . randomize ( ) 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 ( ) 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
Example of regressing with a student t likelihood using Laplace