首页 文章

高斯过程回归

提问于
浏览
2

我正在编写高斯过程回归算法 . 这是代码:

% Data generating function

fh = @(x)(2*cos(2*pi*x/10).*x);

% range

x = -5:0.01:5;
N = length(x);

% Sampled data points from the generating function

M = 50;
selection = boolean(zeros(N,1));
j = randsample(N, M);

% mark them

selection(j) = 1;
Xa = x(j);

% compute the function and extract mean

f = fh(Xa) - mean(fh(Xa));
sigma2 = 1;

% computing the interpolation using all x's
% It is expected that for points used to build the GP cov. matrix, the
% uncertainty is reduced...

K = squareform(pdist(x'));
K = exp(-(0.5*K.^2)/sigma2);

% upper left corner of K

Kaa = K(selection,selection);

% lower right corner of K

Kbb = K(~selection,~selection);

% upper right corner of K

Kab = K(selection,~selection);

% mean of posterior

m = Kab'*inv(Kaa+0.001*eye(M))*f';

% cov. matrix of posterior

D = Kbb - Kab'*inv(Kaa + 0.001*eye(M))*Kab;

% sampling M functions from from GP

[A,B,C] = svd(Kaa);
F0 = A*sqrt(B)*randn(M,M);
% mean from GP using sampled points

F0m = mean(F0,2);
F0d = std(F0,0,2);

%%
% put together data and estimation

F = zeros(N,1);
S = zeros(N,1);
F(selection) = f' + F0m;
S(selection) = F0d;

% sampling M function from posterior

[A,B,C] = svd(D);
a = A*sqrt(B)*randn(N-M,M);
% mean from posterior GPs

Fm = m + mean(a,2);
Fmd = std(a,0,2);
F(~selection) = Fm;
S(~selection) = Fmd;

%%

figure;
% show what we got...

plot(x, F, ':r', x, F-2*S, ':b', x, F+2*S, ':b'), grid on;
hold on;
% show points we got

plot(Xa, f, 'Ok');
% show the whole curve

plot(x, fh(x)-mean(fh(x)), 'k');
grid on;

我希望得到一些不错的数据,其中未知数据点的不确定性会很大,并且采样数据点周围很小 . 我得到了一个奇怪的数字甚至更奇怪的是,采样数据点周围的不确定性大于其余数据点 . 有人可以向我解释我做错了什么吗?谢谢!!

1 回答

  • 3

    您的代码存在一些问题 . 以下是最重要的一点:

    • 使一切出错的主要错误是 f 的索引 . 您正在定义 Xa = x(j) ,但实际上应该执行 Xa = x(selection) ,以便索引与您在内核矩阵 K 上使用的索引一致 .

    • 减去样本均值 f = fh(Xa) - mean(fh(Xa)) 不起任何作用,并使得图中的圆圈与实际函数不同 . (如果你选择减去某个东西,它应该是一个固定的数字或函数,而不是取决于随机抽样的观察结果 . )

    • 你应该直接从 mD 计算后验均值和方差;无需从后验采样,然后获得那些样本估计 .

    这是修改后的脚本版本,修复了以上几点 .

    %% Init
    % Data generating function
    fh = @(x)(2*cos(2*pi*x/10).*x);
    % range
    x = -5:0.01:5;
    N = length(x);
    % Sampled data points from the generating function
    M = 5;
    selection = boolean(zeros(N,1));
    j = randsample(N, M);
    % mark them
    selection(j) = 1;
    Xa = x(selection);
    
    %% GP computations
    % compute the function and extract mean
    f = fh(Xa);
    sigma2 = 2;
    sigma_noise = 0.01;
    var_kernel = 10;
    % computing the interpolation using all x's
    % It is expected that for points used to build the GP cov. matrix, the
    % uncertainty is reduced...
    K = squareform(pdist(x'));
    K = var_kernel*exp(-(0.5*K.^2)/sigma2);
    % upper left corner of K
    Kaa = K(selection,selection);
    % lower right corner of K
    Kbb = K(~selection,~selection);
    % upper right corner of K
    Kab = K(selection,~selection);
    % mean of posterior
    m = Kab'/(Kaa + sigma_noise*eye(M))*f';
    % cov. matrix of posterior
    D = Kbb - Kab'/(Kaa + sigma_noise*eye(M))*Kab;
    
    %% Plot
    figure;
    grid on;
    hold on;
    % GP estimates
    plot(x(~selection), m);
    plot(x(~selection), m + 2*sqrt(diag(D)), 'g-');
    plot(x(~selection), m - 2*sqrt(diag(D)), 'g-');
    % Observations
    plot(Xa, f, 'Ok');
    % True function
    plot(x, fh(x), 'k');
    

    由此产生的结果具有5个随机选择的观察结果,其中真实函数以黑色显示,后验平均值以蓝色显示,置信区间以绿色显示 .

    GP estimates

相关问题