function [fp,s2p_sk,lower95_sk,upper95_sk] = posterior(s2,l,ns,xs,np,xp,fs,TSigma)

r = correlationfunc(l,xs,xp);
s = s2*r;

F = ones(ns,1);
x = TSigma'\fs;
M = TSigma'\F;
b = (M'*M)\(M'*x);  % b = (F'*Sigma^(-1)*F)^(-1)*F'*Sigma^(-1)*fs;
z = x - M * b;
zz = z'*z; % zz = (fs-F*b)'*Sigma^(-1)*(fs-F*b);

f = ones(np,1);
trend = f*b; % trend
complement = (TSigma'\s)' * z; % complement = s'*Sigma^(-1)*(fs-F*b); local departure
fp = trend + complement;

s21 = sum((TSigma'\s) .* (TSigma'\s), 1); % compute s'*Sigma^(-1)*s 
s2p_sk = s2 - s21';
% T_M = chol(M' * M);   % equivalently : qrR <- qr.R(qr(M))
% m = (T_M)' \ (f - (TSigma'\s)' * M)';
% s22 = sum(m .* m,1);
% s2p_uk= s2 - s21' + s22';
s2p_sk(s2p_sk<0) = 0.0;

q95 = norminv(0.975);
lower95_sk = fp - q95*sqrt(s2p_sk);
upper95_sk = fp + q95*sqrt(s2p_sk);
% lower95_uk = fp - q95*sqrt(s2p_uk);
% upper95_uk = fp + q95*sqrt(s2p_uk);