%generating function
f = @(x) (sin(x*10)*10+(x+3).*(x+2).*(x+1).*x.*(x-1).*(x-2).*(x-3)).*(1-tanh((abs(x)-3)*2.5));

%generate training data: noisy samples from true function
se = 50; %observation noise (std)
ntrain = 100; %number of observations
Xtrain = rand(ntrain,1)*10-5; %random training stimuli from interval [-5,5]
Ytrain = f(Xtrain)+randn(ntrain,1)*se; %observed outcomes

%gaussian process model
a = .5; %width of kernel
l = .1; %penalty parameter
k = @(x,y) exp(-((x-y)/a).^2); %gaussian kernel
K = @(x,y) k(repmat(x,1,length(y)),repmat(y,1,length(x))'); %creates kernel matrix between any two sets of stimuli (x and y)
fhat = @(Xtest,Xtrain,Ytrain) K(Xtest,Xtrain)*(K(Xtrain,Xtrain)+l*eye(length(Xtrain)))^(-1)*Ytrain; %estimated function values (i.e., posterior means) for any given test stimuli Xtest, after learning from Xtrain and Ytrain (see first paragraph of section 3.2 in Griffiths et al., 2009)

%do inference and plot
figure(1),clf,hold on
Xtest = (-5:.01:5)'; %test points densely covering stimulus space, to show model's predictions for all possible stimuli
plot(Xtrain,Ytrain,'r.','MarkerSize',20) %plot training data as red points
plot(Xtest,f(Xtest),'k-') %plot true function in black
plot(Xtest,fhat(Xtest,Xtrain,Ytrain),'g-') %plot estimated function (for all test points) in green