%?Calculates?for?xstar?x-coordinates?the?bestEstimate?and?the?95%
%?confidence?interval?（bounds）?for?a?Gaussian?process?regression.
%?Mark?Ebden?July?2008
%
%?Inputs:
%??k?-?a?handle?to?the?covariance?function
%??Xy?-?the?training?data?x-?and?y-coordinates
%??theta?-?the?GP?parameters?to?be?passed?to?the?covariance?function
%??????????（e.g.?length?scale?periodicity?etc.
%??sigma_n?-?the?noise?（optional?--?default?0）
%????NOTE:?sigma_n?should?be?left?as?zero?if?k?*already*?includes?the?noise?term?as?I?usually?do.
%??xstar?-?specify?the?x-coordinates?to?solve?for?（optional?--?default?is
%??????????to?let?the?algorithm?choose?them）
%??graphing?-?produce?a?plot?of?results?（optional）
%?????????????Default?is?to?plot?the?results?unless?the?function?outputs?were?asked?for
%
%?Outputs:
%??xstar?-?x-coordinates?（default?100?are?chosen）
%??K?-?covariance?matrix?among?training?data
%??V?-?variance?at?each?xstar

function?[xstar?bestEstimate?bounds?K?V]?=?calcGP?（k?X?y?theta?sigma_n?xstar?graphing）

if?nargin?????if?nargout?==?0
????????graphing?=?1;?
????else
????????graphing?=?0;
????end
end
if?nargin?????r2?=（max（X）-min（X））;
????N?=?1e3;?%?How?many?points?to?use
????z?=?.5;?%?How?much?to?extend?the?time?series?on?either?side?（e.g.?.5?is?50%?extension?on?left?and?on?right）
????xstar?=?（0:N）/N?*?（z*2+1）*r2?+?min（X）-r2*z;
end
if?nargin?????sigma_n?=?0;
end
if?nargin?????theta?=?[1?0?1?0];
end

%?a）?Initializations
meany?=?mean（y）;?y?=?y?-?meany;
n?=?size（X1）;?lx?=?length（xstar）;
K?=?zeros（n）;?kstar?=?zeros（n1）;
for?i?=?1:n
????for?j?=?1:n
????????K（ij）?=?k?（X（i）X（j）theta）;
????end
end
fstarbar?=?zeros（lx1）;?V?=?zeros（lx1）;

%?b）?One-off?calculations
diags?=?max（1e3*eps?sigma_n^2）;?%?beef?up?the?diagonal?if?sigma_n?=?0
L?=?chol?（K?+?diags*eye（n）‘lower‘）;?
alpha?=?L‘\（L\y）;
logpyX?=?-y‘*alpha/2?-?sum（log（diag（L）））?-?n*log（2*pi）/2;?%?Log?marginal?likelihood

%?c）?xstar?loop
for?q?=?1:lx
????for?i?=?1:n
????????kstar（i）?=?k?（X（i）xstar（q）theta）;
????end
????kstar
????%?Mean?of?prediction
????fstarbar（q）?=?kstar‘?*?alpha;
????%?Variance?of?prediction
????v?=?L\kstar;
????ystar_noise?=?sigma_n^2;?%?recall?f*?has?no?noise?itself?（Algorithm?2.1?on?p?19?of?H152）
????V（q）?=?k?（xstar（q）xstar（q）theta）?-?v‘*v?+?ystar_noise;
end
bounds?=?[fstarbar+1.96*sqrt（V）?fstarbar-1.96*sqrt（V）]+meany;
bestEstimate?=?fstarbar+meany;

%?d）?Output
if?graphing?==?1
????figure?hold?on
????stdRegion?（xstarbounds）
????plot?（xstarbestEstimate‘k‘）
????plot?（Xy+meany‘b+‘）
end

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????2919??2012-10-27?12:15??gaussian\GPtutorial.m

?????文件???????2043??2011-03-31?11:51??gaussian\GPtutorialFcn.m

?????文件????????509??2011-02-03?12:23??gaussian\hypSample.m

?????文件???????1436??2010-11-09?18:03??gaussian\k_GP.m

?????文件????????599??2012-10-27?12:16??gaussian\stdRegion.m

?????文件???????2551??2012-04-27?14:23??gaussian\calcGP.m

?????目錄??????????0??2014-03-13?15:35??gaussian

-----------?---------??----------?-----??----

????????????????10057????????????????????7