n = 1000;
sigma = 5;
wn = normrnd(0,sigma,n,1);
  figure(1)
plot(wn)
  figure(2)
plot(xcov(wn))

bm = cumsum(wn);
  figure(1)
plot(bm)
  figure(2)
plot(xcov(bm))


data = load('./data/Kit2_2.dat'); 

V = log(data(:,2)+0.01);   x = data(:,3);  y = data(:,4);    % log TCE concentrations  
n = length(V);

plot3(x,y,V,'.');

X = [ones(n,1) x y  x.*y  x.^2  y.^2]; 

[xx,yy] = meshgrid(0:10:370,-90:5:-45);

betahat = (X' * X)^(-1) * X' * V;
pred = betahat(1) + betahat(2)*xx + betahat(3)* yy  +  betahat(4)* xx.*yy + betahat(5)* xx.^2  + betahat(6)* yy.^2;
  hold on
surf(xx,yy,pred)
 hold off

Vhat = X * betahat;    % predicted values 
resid = V - Vhat;

g =  variogram(x, y, resid,0);

%%%%%%%%%%%%%%%%%% 
% now fitting Gaussian variogram

global pts vario np;        % we use them later for the function wsq

sz = size(g);
nzones = sz(2);
 pts = g(2:end,3);           % bin centers
 vario = g(2:end,4);         % values of gamma
 np =  g(2:end,5);           % number of points in each variogram bin

 
 
% we will grid-search in the vicinity of range = 10, sill = 2 

 range = 3:0.5:9;
 sill = 1.5:0.1:2.5;
 crit = zeros(length(range), length(sill));
 for i = 1:length(range), 
    for j = 1:length(sill),  
        crit(i,j) = wsq([range(i) sill(j)]);     % , pts, vario, np);       
    end
 end    
 crit
 imagesc(-crit);   colormap gray;
  
 [pars,fval] = fminsearch(@wsq,[6, 2])
 
 
 plot(pts, vario,'*r');
   hold on
 xc = 0:max(pts);
 yc = pars(2) * (1 - exp(-(xc/pars(1)).^2));
 plot(xc, yc)
   hold off
   
 
 
 
 
   
 
 
 
 






% generating random fields from Exponetial, Gaussian and Spherical models

randomfield   % this initiates a pretty nice graphical interface 

%==================================================
% #5 in Lab 2

data = load('c:/local/586/data/hi_plain.txt');

Y = data(:,9);   East = data(:,6);  North = data(:,7);
n = length(Y);

X = [ones(n,1) East North]; 

plot3(East, North, Y, '.');

betahat = (X' * X)^(-1) * X' * Y;
Yhat = X * betahat;    % predicted values 
resid = Y - Yhat;

  % there is a slight curvature going on but the  linear model is probably acceptable

g =  variogram(East, North, resid,1);

   % for example, I entered 10 lags, lag step = 3, vario. angle = 54. There
   % seems to be a mistake in this code: not all points are good for aniso.  model