高斯过程是一种非参数模型估计方法。不像最小二乘,需要知道模型的参数,如:y=ax+b,我们就需要知道a和b来对模型进行估计。
高斯过程要设置一个核函数,来给不同观测值确定关系。这里我们需要设置核函数的超参数,比如下面的alpha和beta。
下面是几种常见的计算不同观测关系的核函数:
设置好核函数后我们需要计算三组矩阵,分别是自观测矩阵K,观测预测矩阵Ks和预测矩阵Kss。
计算方法就是将观测自变量和预测自变量代入设定的核函数即可,具体可以参考后面代码,不是很麻烦。
得到三个矩阵后只需要根据下面两个公式就能得到预测数据的均值和方差。
其中Y是观测因变量,sigma是观测可能的不确定性。
代码如下:
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);
结果如下: