matlab练习程序（高斯过程回归）

clear all;
close all;
clc;

x=0:0.1:8;              %原始数据
y = x.*sin(x) + x/5;    
plot(x,y,‘r‘)

gc_x = x(1:end);        %观测数据加噪声
gc_y = y(1:end) + (rand(1,length(gc_x))-0.5)/3;
hold on;
plot(gc_x,gc_y,‘ro‘)

x=0:0.05:8;             %带预测数据
num = length(gc_x);
K = zeros(num,num);                 %根据核函数给出自观测矩阵
Ks = zeros(num,length(x));          %根据核函数给出观测预测矩阵
Kss = zeros(length(x),length(x));   %根据核函数给出预测矩阵

%根据核函数赋值
for i=1:num
    for j=1:num
        K(i,j) =exp(-((gc_x(i)-gc_x(j))*(gc_x(i)-gc_x(j))));        
    end
    for j=1:length(x)
        Ks(i,j) = exp(-((gc_x(i)-x(j))*(gc_x(i)-x(j))));        
    end
end

for i=1:length(x)
    for j=1:length(x)
        Kss(i,j) = exp(-((x(i)-x(j))*(x(i)-x(j))));
    end
end

%计算均值与方差
sigma = 0.001;  %sigma改成2可以明显看出方差范围    
mean_y = Ks‘*inv(K + sigma*eye(length(K)))*gc_y‘;
std_y = diag(Kss - Ks‘*inv(K + sigma*eye(length(K)))*Ks);
plot(x,mean_y,‘c‘);

%画出方差范围
max_y = mean_y + abs(std_y);
min_y = mean_y - abs(std_y);
yFill = [max_y‘, fliplr(min_y‘)];
xFill = [x, fliplr(x)];

plot(x,max_y,‘c‘);
plot(x,min_y,‘c‘);
fill(xFill, yFill, ‘c‘,‘edgealpha‘,0,‘facealpha‘,0.1);