gaussian process

gp_regression.m

kernel = 2;
switch kernel
    case 1; k = @(x, y) 1 * (x'*y + 1);
    case 2; k = @(x, y) exp(-1 * (x-y)' * (x-y));
    case 3; k = @(x, y) 1 * (x'*y + 1) * (x'*y + 1) + x'*y+1
    case 4; k = @(x, y) 1 * min(x, y)
end

x_test = (0:.005:3)';
n_test = length(x_test);

n_train = 30;
x_train = x_test(randperm(n_test)(1:n_train));
f = @(x) 1 + cos(x) + sin(x);
y_train = f(x_train) + randn(size(x_train)) * sqrt(0.05);

n = n_test + n_train;
x = [x_test; x_train];

C = zeros(n,n);
for i = 1 : n
    for j = 1 : n
        C(i,j) = k(x(i), x(j));
    end
end
C = C + 1 * eye(n);

Caa = C(1:n_test, 1:n_test);
Cab = C(1:n_test, n_test+1 : n);
Cba = C(n_test+1:n, 1:n_test);
Cbb = C(n_test+1:n, n_test+1:n);

m = Cab * inv(Cbb) * (y_train);
D = Caa - Cab * inv(Cbb)*Cba;

u = randn(n_test, 1);
[A, S, B] = svd(D);
z = m + A * sqrt(S) * u;

figure(2); hold on; clf
fitted_curve = plot(x_test, m, 'b-');
hold on;
true_curve = plot(x_test, f(x_test), 'g-');
hold on;
train_points = plot(x_train, y_train, 'xr');
m1 = "fitted curve";
m2 = "true curve";
m3 = "training points";
legend([fitted_curve; true_curve; train_points], [m1; m2; m3]);
title("linear kernel");

gp_sample.m

kernel = 4;
switch kernel
    case 1; k = @(x, y) 1 * x'*y;
    case 2; k = @(x, y) exp(-10 * (x-y)' * (x - y));
    case 3; k = @(x, y) 1 * (x'*y) * (x'*y) + x'*y;
    case 4; k = @(x, y) 1 * min(x, y);
    case 5; k = @(x, y) exp(sin(5*pi*(x-y))^2);
end

x = (0:.005:1);
n = length(x);

C = zeros(n,n);

for i = 1 : n
    for j = 1 : n
        C(i,j) = k(x(i), x(j));
    end
end

u = randn(n, 1);
[A, S, B] = svd(C);
z = A * sqrt(S) * u;

figure(2); hold on; %clf
plot(x, z, 'g-');