% Example 9.1
fLF= @(x) 0.5*sin(3*x); fHF= @(x) sin(3*x) + 0.25*sin(0.3*x);
xtrainLF=[0.2,0.3,0.4,0.5,0.6,0.7,0.8];
ytrainLF=fLF(xtrainLF);                               % LF training samples
xtrainHF=[0.2,0.4,0.6];
ytrainHF=fHF(xtrainHF);                               % HF training samples
xplot=0:0.01:1; yplotLF=fLF(xplot); yplotHF=fHF(xplot);
% Low-fidelity surrogate (quadratic PRS)
X_LF=[ones(length(xtrainLF),1) xtrainLF' xtrainLF.^2'];
[B_LF,BINT_LF,R_LF,RINT_LF,STATS_LF]=regress(ytrainLF',X_LF);
y_LF=[ones(length(xplot),1) xplot' xplot.^2']*B_LF;
% High-fidelity surrogate (linear PRS)
X_HF=[ones(length(xtrainHF),1) xtrainHF'];
[B_HF,BINT_HF,R_HF,RINT_HF,STATS_HF]=regress(ytrainHF',X_HF);
y_HF=[ones(length(xplot),1) xplot']*B_HF;
%
% Multi-fidelity surrogate (linear discrepancy PRS)
ydiff=fHF(xtrainHF)-fLF(xtrainHF);
X_diff=[ones(length(xtrainHF),1) xtrainHF'];
[B_diff,BINT_diff,R_diff,RINT_diff,STATS_diff]=regress(ydiff',X_diff);
y_MF=[ones(length(xplot),1) xplot' xplot.^2']*B_LF ...
    +[ones(length(xplot),1) xplot']*B_diff;
%
plot(xplot,yplotLF,'b',xplot,yplotHF,'r',xplot,y_LF,'b--',...
     xplot,y_HF,'r--',xplot,y_MF,'k--',xtrainLF,ytrainLF,'bo',...
     xtrainHF,ytrainHF,'r*');
legend('LF function','HF function','LF prediction',...
       'HF prediction', 'MF prediction');
%-----------------------------------------------------------------------------
%
% Example 9.2
X_comp=[[ones(length(xtrainHF),1) xtrainHF' xtrainHF.^2']*B_LF ...
        ones(length(xtrainHF),1) xtrainHF'];
[B_comp,BINT_comp,R_comp,RINT_comp,STATS_comp]=regress(ytrainHF',X_comp);
y_MF_comp=[[ones(length(xplot),1) xplot' xplot.^2']*B_LF ...
           ones(length(xplot),1) xplot']*B_comp;
%-----------------------------------------------------------------------------
%
% Example 9.6
lowfun = @(x) 0.5*((6*x-2).^2.*sin(12*x-4))+10*(x-0.5)+5;
highfun = @(x) (6*x-2).^2.*sin(12*x-4);
xs=[0.5862,0.3638,0.8198]';nHy=size(xs,1);     % HF sample locations & size
yLs=lowfun(xs);yHs=highfun(xs);                           % LF & HF samples
x=(0:0.01:1)';N=size(x,1);                              % Prediction points
yLtrue=lowfun(x);yHtrue=highfun(x);                % True LF & HF functions
[yDp,yDsig2,rho,nbD] = MFSLF(xs,yHs,yLs,x,0,0);      % fitting Bayesian MFS
%
yHp=rho*yLtrue+yDp;                            % Multi-fidelity predictions
yHsig=rho*sqrt(yDsig2);                                    % Prediction STD
PI = yHp + tinv([0.025 0.975],nHy-nbD).*yHsig;       % Prediction intervals
figure(1); hold on;                          % Plotting result MF surrogate
xShaded=[x;sort(x,'descend')]; yShaded=[PI(:,1);PI(end:-1:1,2)];
patch(xShaded,yShaded,[.24 .73 .89],'EdgeColor','none');
plot(x,yHp,'r-',xs,yHs,'k+',x,yHtrue,'k-','LineWidth',2,'Markersize',10);
legend('95% PI','y_M(x) prediction','y_H(x) samples','y_H(x) true');
xlabel('x'); set(gca,'Fontsize',14)
%
yDsig=sqrt(yDsig2); yDtrue=yHtrue-rho*yLtrue;% Discrepancy STD & prediction
PI = yDp + tinv([0.025 0.975],nHy-nbD).*yDsig;       % Prediction intervals
figure(2); hold on;                        % Plotting discrepancy surrogate
xShaded=[x;sort(x,'descend')]; yShaded=[PI(:,1);PI(end:-1:1,2)];
patch(xShaded,yShaded,[.24 .73 .89],'EdgeColor','none');
plot(x,yDp,'r-',xs,yHs-rho*yLs,'k+',x,yDtrue,'k-','LineWidth',2,'Markersize',10)
legend('95% PI','y_D(x) prediction','y_D(x) samples','y_D(x) true');
xlabel('x'); set(gca,'Fontsize',14)
%
eRMS=sqrt(sum((yHtrue-yHp).^2)/N);             % RMS error at N grid points
eav=sum(abs(yHtrue-yHp))/N;              % Average error at NxN grid points
emax=max(abs(yHtrue-yHp));               % Maximum error at NxN grid points
fprintf('RMS error:%8.3g\nAve error:%8.3g\nMax error:%8.3g\n',eRMS,eav,emax);
%-----------------------------------------------------------------------------
%
% Example 9.7
lowfun = @(x) 0.5*((6*x-2).^2.*sin(12*x-4))+10*(x-0.5)+5;
highfun = @(x) (6*x-2).^2.*sin(12*x-4);
xLs=[0.1548,0.7148,0.4934,0.9932,0.5862,0.0007,0.3638,0.8198]';     % LF inputs
xHs=[0.5862,0.0007,0.3638,0.8198]';                                 % HF inputs
yLs=lowfun(xLs);                                         % low-fidelity outputs
yHs=highfun(xHs);                                       % high-fidelity outputs
trL=0;                                      % LF trend function order (0, 1, 2)
trD=0;                             % Discrepancy trend function order (0, 1, 2)
x=(0:0.01:1)';
nHy=size(yHs,1);nLy=size(yLs,1);
% fitting Bayesian MFS
[yLp,yLsig2,yDp,yDsig2,rho,nbL,nbH] = MFS(xLs,yLs,xHs,yHs,x,trL,trD,0);
yHp=rho*yLp+yDp;
yHsig=sqrt(rho^2*yLsig2+yDsig2);
cilb = yHp + tinv(0.025,nHy-nbH)*yHsig;                    % Lower-bounds of PI
ciub = yHp + tinv(0.975,nHy-nbH)*yHsig;                    % Upper-bounds of PI
yHtrue=highfun(x);
% plotting result 
figure(1); hold on;
xShaded=[x;sort(x,'descend')]; yShaded=[cilb;ciub(end:-1:1)];
patch(xShaded,yShaded,[.24 .73 .89],'EdgeColor','none');
plot(x,yHp,'r-','LineWidth',2)
plot(xHs,yHs,'k+','LineWidth',2,'Markersize',10)
plot(x,yHtrue,'k-','lineWidth',2)
legend('95% CI','y_M(x) prediction','y_H(x) samples','y_H(x) true');
xlabel('x'); set(gca,'Fontsize',14)
%
yLtrue=lowfun(x); yLsig=sqrt(yLsig2);
cilb = yLp + tinv(0.025,nLy-nbL)*yLsig;                    % Lower-bounds of PI
ciub = yLp + tinv(0.975,nLy-nbL)*yLsig;                    % Upper-bounds of PI
figure(2); hold on;
xShaded=[x;sort(x,'descend')]; yShaded=[cilb;ciub(end:-1:1)];
patch(xShaded,yShaded,[.24 .73 .89],'EdgeColor','none');
plot(x,yLp,'r-','LineWidth',2)
plot(xLs,yLs,'k+','LineWidth',2,'Markersize',10)
plot(x,yLtrue,'k-','lineWidth',2)
legend('95% CI','y_L(x) prediction','y_L(x) samples','y_L(x) true');
xlabel('x'); set(gca,'Fontsize',14)
%
figure(3); hold on;
yDsig=sqrt(yDsig2); yDtrue=yHtrue-rho*yLtrue;
cilb = yDp + tinv(0.025,nHy-nbH)*yDsig;                    % Lower-bounds of PI
ciub = yDp + tinv(0.975,nHy-nbH)*yDsig;                    % Upper-bounds of PI
xShaded=[x;sort(x,'descend')]; yShaded=[cilb;ciub(end:-1:1)];
patch(xShaded,yShaded,[.24 .73 .89],'EdgeColor','none');
plot(x,yDp,'r-','LineWidth',2)
plot(xHs,yHs-rho*lowfun(xHs),'k+','LineWidth',2,'Markersize',10)
plot(x,yDtrue,'k-','lineWidth',2)
legend('95% CI','y_D(x) prediction','y_D(x) samples','y_D(x) true');
xlabel('x'); set(gca,'Fontsize',14)
%-----------------------------------------------------------------------------