Sampling from Posterior Distribution of GPR Model from fitrgp()

Question

Sterling Baird 2021-1-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/709828-sampling-from-posterior-distribution-of-gpr-model-from-fitrgp

评论： Jitin Malhotra 2021-6-25

I can get a gprMdl via:

rng(0,'twister');
N = 100;
x = linspace(-10,10,N)';
y = 1 + x*5e-2 + sin(x)./x + 0.2*randn(N,1);
gprMdl = fitrgp(x,y,'FitMethod','Exact','PredictMethod','Exact');

and the posterior covariance matrix via:

kfcn = gprMdl.Impl.Kernel.makeKernelAsFunctionOfXNXM(gprMdl.Impl.ThetaHat)
K = kfcn(x(1:5,:),x(1:7,:))

(credit: Gautam Pendse in MATLAB Answers)

How can I go about sampling a model from the posterior distribution of gprMdl?

To clarify, I mean like the curves from the 2nd tile of:

Source: https://commons.wikimedia.org/wiki/File:Gaussian_Process_Regression.png

Being able to do this is a matter of being able to finish up a graduate program, so any help is much appreciated!

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Sterling Baird 2021-1-6

编辑：Sterling Baird 2021-1-6

Perhaps all I need is to use the mean from predict and the covariance matrix from the undocumented function along with mvnrnd, similar to what's described in: https://www.mathworks.com/matlabcentral/answers/442492-sampling-based-on-a-correlation-matrix-for-multiple-parameters

Though I'm not sure if this is sampling from the posterior..

Thoughts?

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Gautam Pendse 2021-1-7

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/709828-sampling-from-posterior-distribution-of-gpr-model-from-fitrgp#answer_592493

在 MATLAB Online 中打开

Hi Sterling,

Here's an example illustrating how to sample from the posterior distribution of a GPR model. The code uses an undocumented function predictExactWithCov. If you have categorical predictors, you would need to convert them to dummy variables before using that function.

Hope this helps,

Gautam

function testPosteriorCovarianceGP()
% 1. Some dummy data.
rng(0,'twister');
X = [-10,-8,-5,-1,1,3,7,10]';
N = numel(X);
y = 1 - X*5e-2 + sin(X)./X + 1e-4*randn(N,1);
M = 1000;
Xnew = linspace(-10,10,M)';
% 2. Fit a GPR model.
gpr = fitrgp(X,y,'Sigma',0.01,'SigmaLowerBound',1e-6);
% 3. Access the posterior covariance matrix by calling the undocumented
% predictExactWithCov method on the gpr.Impl object. Vectors ynew can be
% simulated from a Normal distribution with mean pred and covariance covmat
% using cholcov or svd and randn.
alpha = 0.05;
[pred,covmat,ci] = predictExactWithCov(gpr.Impl,Xnew,alpha);
figure(1);
plot(Xnew,pred,'b');
hold on;
plot(Xnew,ci(:,1),'m');
plot(Xnew,ci(:,2),'g');
plot(X,y,'ko','MarkerFaceColor','k');
xlabel('x');
ylabel('y');
figure(2);
T = cholcov(covmat);
numSamples = 50;
ynew = pred + T'*randn(M,numSamples);
plot(Xnew,ynew);
hold on;
plot(X,y,'ko','MarkerFaceColor','k');
xlabel('x');
ylabel('y');
end

2 个评论
显示无隐藏无

Sterling Baird 2021-1-14

编辑：Sterling Baird 2021-1-14

What is alpha? (Edit: I just realized it controls the confidence intervals)

Jitin Malhotra 2021-6-25

@Gautam Pendse, I din't understand what exactly you did and how i should use it on my data. If possible can you please provide more details about this undocumented function? predictExactWithCov(gpr.Impl,Xnew,alpha);

My main motive is to calculate CI_width and CI_var, which is basically the 95% confidence interval width and variance acc. to given formulae