% eda10_04
% simple example of kriging

% the model parameters are the values of a sinusoidal
% function evenly spaced along the t axis
M=101;
Dt=1.0;
tmin=0;
t = tmin+Dt*[0:M-1]';
tmin=t(1);
tmax=t(M);

% the data are the vales the function measured on
% only a subset of points along the t axis
N=50;
% randomly choose N of the M points to have data
rowindex=sort(unidrnd(M,N,1));
% determine their t positions
tobs=tmin+Dt*rowindex;
% evaluate a sinusoidal function at these points to
% simulate the data
freq=0.02;
dobs = sin(2*pi*freq*tobs);

% estimate data at these points
test=t;
dest=zeros(M,1);

% plot the observed data with circles
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis([tmin tmax -2 2]);
plot(tobs,dobs,'ko','LineWidth',2);
xlabel('t');
ylabel('d');

% assume the autocorrelation function of the data
% is a a Normal function with variance L^2.  The normalization
% constant cancels, so just set to unity

% set the constant L arbitrarily
L2 = 10;

% basics equations are Aw=v(i), dest(i)=dobs'*w(i);

% set up the kriging matrix, A
A = exp(-abs(tobs*ones(N,1)'-ones(N,1)*tobs').^2 /(2*L2));

% set up kriging vector, v, but with vector for each interplation
% point a column of the matrix, V
V = exp(-abs(tobs*ones(M,1)'-ones(N,1)* test').^2 /(2*L2));

% solve, but add some damping, just in case A
% is near singular.  This selects the smallest w
% that solves the equation
dest=dobs'*((A+1e-6*eye(N))\V);
    
% plot estimated data with line
plot(test,dest,'k-','LineWidth',2);

