%% РпупфйкпрпЯзуз ЧщсйкЮт УхнЬцейбт уе УзмейбкЬ ДедпмЭнб
% *ЦбЯдщн КхсйбкЯдзт*
%
% ФмЮмб ГещгсбцЯбт РбнерйуфзмЯпх БйгбЯпх

%% ЕйубгщгЮ дедпмЭнщн уе узмеЯб кбй рлзспцпсЯет
load BayAreaData.mat
who

%% РспкбфбскфйкЬ
Data = bayprecip
% 1з уфЮлз -> бэопнфбт бсйимьт вспчпмефсйкпэ уфбимпэ
% 2з уфЮлз -> гещгсбцйкь мЮкпт уфбимпэ
% 3з уфЮлз -> гещгсбцйкь рлЬфпт уфбимпэ
% 4з уфЮлз -> вспчьрфщуз уе mm
% 5з уфЮлз -> вспчьрфщуз (уе mm) рпх кбфбгсЬцзке брь Энб Ьллп ьсгбнп
nDat = size(bayprecip,1);
idStation = num2str(Data(:,1));
% ухмвплпуейсЬ ме фпн бэопнфб бсйимь фпх кЬие уфбимпэ
plotPars = {'color', 'vbar'}; 
% color -> Эгчсщмз брейкьнйуз, vbar -> ме кЬиефп хрьмнзмб чсюмбфпт

%% Брейкьнйуз
bulletSize = 20;
cLims = [1 15]; % чсщмбфйкЬ ьсйб
scattermap(bayprecip,[2 3],bulletSize,4,'cont',cLims,plotPars);
% [2 3] -> уфЮлет ме ухнфефбгмЭнет уфбимюн
% 4 -> уфЮлз ме мефсЮуейт уфбимюн
title('Bay Area precipitation data set');
hold on
plot(baycoast(:,1),baycoast(:,2),'k-');
hold off

%% Емрейсйкь змйвбсйьгсбммб ме фзн енфплЮ points2variogram
help spassociation
%dirSpecs = [];  % ден ендйбцЭсей з дйеэихнуз (рбн-дйехихнуйбкь)
LagMidTol = []; % бхфьмбфз ерйлпгЮ дйбуфзмЬфщн брьуфбузт кбй бнпчЮт
G = spassociation(bayprecip,[2 3],4,LagMidTol,'variogram')
%G.assoc

%% Брейкьнйуз змйвбсйпгсЬммбфпт ме фзн енфплЮ plotvariogram
help plotspassociation
plotspassociation(G)
dir2plot = []; % Брейкьнйуз ьлщн фщн вбсйпгсбммЬфщн
plotspassociation(G,[],'r*','pairs');

%% ЕобгщгЮ дедпмЭнщн гйб рбсЬдейгмб
lagObs = G.assoc(:,1); % dist  -> ч брпуфЬуейт
gamObs = G.assoc(:,2); % gam   -> ш фймЭт вбсйпгсЬммбфпт
nPairs = G.assoc(:,3); % pairs -> бсйимьт жехгбсйюн гйб кЬие брьуфбуз
% брблейцЮ 1пх жеэгпхт (0,0)
lagObs(1) = [];
gamObs(1) = [];
nPairs(1) = [];
plot(lagObs,gamObs,'*'); ylim([0 12]);
xlabel('distance: lag');
ylabel('semivariogram: gamObs');
title('A scatter plot of y (dependent) vs x (independent variable)');

%% Function handle for Exponential semivariogram model
% Create a function handle for the exponential semivariogram model as:
 funGamExp = @(b,x) b(1)*(1-exp(-3*x/b(2)))
%% 
% where |x| stands for an input distance and |b = [b(1) b(2)]| 
% for a parameter vector; type:
whos funGamExp 
%%
% to see some properties of this data type. 
%
% You can now construct a finely discretized distance vector |lagFine|, 
% say: 
lag = (0.01:0.01:1.5);
disp(lag(1:5));
%% 
% and compute the associated semivariogram values for an 
% exponential model with sill 9 and range 1.2, i.e., |b = [10 1]|, as:
b = [10 1];
gamMod = funGamExp(b,lag);
plot(lag,gamMod,'.'); 
title('Exponential semivariogram model: sill = 10, range = 1');
hold on
plot(lagObs,gamObs,'*');
gamModObs = funGamExp(b,lagObs); % model evaluated at data distances
plot(lagObs,gamModObs,'ro');
hold off
SSE = sum((gamObs-gamModObs).^2)

%% 
% For a sill of 9 and a range of 1.2, we get:
b = [9 1.2];
gamMod = funGamExp(b,lag); 
plot(lag,gamMod,'.');
title('Exponential model: sill = 9, range = 1.2');
hold on
plot(lagObs,gamObs,'*');
gamModObs = funGamExp(b,lagObs); % model evaluated at data distances
plot(lagObs,gamModObs,'ro');
hold off
SSE = sum((gamObs-gamModObs).^2)

%% Nugget+Exponential semivariogram model
% You can add a nugget effect contribution to the above model, by 
% creating a new function handle as: 
funGamNugExp = @(b,x) b(1)*(x~=0) + b(2)*(1-exp(-3*x/b(3))) 
%% 
% where now |b(1)| is the nugget sill, |b(2)| is the exponential sill,
% and |b(3)| is the range; note that the total sill for this model is:
% |b(1)+b(2)|. For the same distance vector |lagFine|, compute and display
% the corresponding semivariogram values for |b = [0.5 8.5 1.5]| as:
b = [0.5 8.5 1.5];
gamMod = funGamNugExp(b,lag);
plot(lag,gamMod,'.');
title('Exponential model: sill=8.5, range=1.5 + nugget with sill=0.5');
hold on
plot(lagObs,gamObs,'*');
gamModObs = funGamNugExp(b,lagObs); % model evaluated at data distances
plot(lagObs,gamModObs,'ro');
hold off
SSE = sum((gamObs-gamModObs).^2)

%% Fitting an Exponential variogram model with function |nlinfit|
% You can use function |nlinfit| and the function handle |funGamExp| to 
% estimate the parameter vector |bHat| of an exponential semivariogram 
% model that best fits the sample semivariogram in |gamObs|, as:
funGamExp = @(b,x) b(1)*(1-exp(-3*x/b(2)))
bInit = [5 1];  % initial guess for the parameter vector
bHat = nlinfit(lagObs,gamObs,funGamExp,bInit)
%%
% Compute the fitted model semivariogram values at the 30 distances 
% in vector |lag|, as: 
gamHat = funGamExp(bHat,lagObs);
disp([lagObs gamObs gamHat]);
%%
% and overlay the results on the sample semivariogram plot as: 
plot(lagObs,gamObs,'*',lagObs,gamHat,'ro'); ylim([0 12]);
legend('Sample','Model','Location','SouthEast');
title('Sample semivariogram of precip and fitted model');
%%
% Compute the model residuals
rHat = gamObs - gamHat;
disp( mean(rHat) );
%%
% the sum of squared errors (SSE):
SSE = sum(rHat.^2)
%%
% the mean squared error (MSE):
nLag = numel(lagObs);
nPar = numel(bHat);
MSE = SSE/(nLag - nPar) % degrees of freedom = (nLag-nPar)
%%
% and the Akaike information criterion (AIC), as:
AIC = nLag*(log(2*pi*SSE/nLag)+1) + 2*(nPar+1)
%%
% We can now verify that the MSE value computed above is the same as the
% one obtained from function |nlinfit| as:
[b,r,Covb,J,mse2] = nlinfit(lagObs,gamObs,funGamExp,bInit);
disp(mse2);

%% Fitting a Nugget+Exponential variogram model with |nlinfit|
funGamNuggExp = @(b,x) b(1)*(x~=0) + b(2)*(1-exp(-3*x/b(3)))
bInit3 = [0.1 8 0.5];  % initial guess for the parameter vector
[bHat3,r3,Covb3,J3,mse3] = nlinfit(lagObs,gamObs,funGamNuggExp,bInit3);
disp(mse3);
%% 
% Compute the fitted model semivariogram values at the 30 distances 
% in vector |lag|, as: 
gamHat3 = funGamNuggExp(bHat3,lagObs); 
%%
% and overlay the results on the sample semivariogram plot as: 
plot(lagObs,gamObs,'*',lagObs,gamHat3,'ro'); ylim([0 12]);
legend('Sample','Model','Location','SouthEast');
title('Sample semivariogram of precip and fitted model (w/ nugget)');

%% Physical validity (constraints) of fitted parameters
% Function |nlinfit| outputs a negative nugget parameter; this 
% is not physically possible, but nowhere did we instruct the function that
% this was the case. To handle this issue we need to do one or more of the 
% following: (i) change the initial guess parameter vector, hoping to
% find a combination of initial guesses that does not lead to non-sensical
% parameter estimates, (ii) change our model structure, or (iii) keep the
% same model structure and resort to constrained non-linear minimization; 
% see function |fmincon| from the Optimization Toolbox.

%% Recap
% We can conclude that the simpler Exponential-only model should
% be our selected model. Note, however, that in many cases the nugget 
% component is included, if only to account for small-scale attribute 
% variability and measurement error. It is therefore advisable to 
% consider a Nugget+Exponential model in your subsequent analyses.

%% Handling equally plausible (in terms of fit) variogram models
% When multiple variogram models (possibly with different functional forms)
% fit equally well a particular sample semivariogram, one needs to 
% use additional information to resolve the ambiguity. Such information
% could stem from spatial theories dictating the smoothness of the
% particular attribute under study. For example, if an attribute has a
% diffusive-type spatial variation, the appropriate variogram model might
% be the Gaussian. 

%% Weighted least squares for model fitting
% Consider a weights vector with nLag entries defined as: 
w = nPairs./(lagObs.^2); 
% standardize so that weights sum to nLags
w = w./mean(w);
%%
% this weighting scheme gives more weight to small lag distances 
% associated with many data pairs; other possible weighting schemes 
% exist, too. Note that the actual magnitude of the weights does 
% not matter, just their relative magnitude.

%% Weighted fit for the exponential model
% We can create a weighted version of the observations, i.e., a weighted 
% sample semivariogram as:
gamObsW = sqrt(w).*gamObs;
% as well as weighted model predictions as:
funGamExpW = @(b,x) sqrt(w).*funGamExp(b,x); 
%%
% The corresponding parameter vector can be obtained as:
bInit = [5 1];  % initial guess for the parameter vector
bHatW = nlinfit(lagObs,gamObsW,funGamExpW,bInit);
disp(bHatW);
%%
% Note the difference in estimated model parameters between the
% equally and unequally weighted case:
disp('Estimated sill parameters (non-weighted & weighted)');
disp([bHat(1) bHatW(1)]);
disp('Sample variance');
disp(var(bayprecip(:,4),1));
disp('Estimated range parameters (non-weighted & weighted)');
disp([bHat(2) bHatW(2)]);
%%
% Compute the fitted model semivariogram values at the nLag distances 
% in vector |lagObs|, as: 
gamHatW = funGamExp(bHatW,lagObs);
%% 
% *note* the use of |funGamExp| not |funGamExpW| – and overlay the
% results on the sample semivariogram plot as: 
plot(lagObs,gamObs,'*',lagObs,gamHat,'ro',lagObs,gamHatW,'k+'); ylim([0 12]);
legend('Sample','Fitted','Fitted-W','Location','SouthEast');
title('Sample semivariogram & model fit (with & without weights)');
%%
% Compute the weighted residuals:
rHatW = sqrt(w).*(gamObs - gamHatW);
%%
% the weighted sum of squared errors (WSSE):
WSSE = sum(rHatW.^2)
%%
% the mean squared error (MSE):
WMSE = WSSE/(sum(w) - numel(bHat))

%% Variogram model fitting with function fitvariogramisowls
help fitvariogramisowls
iCell = 1; % cell with directional variogram
ModSpecsInit = {'exponential',[5 1]};
wType = 1; % equal weigthing
Out = fitvariogramisowls(G,iCell,ModSpecsInit,wType,'plot');
disp(Out);
%%
Out.specs % array with model type and estimated model parameters
Out.specs{1}
% Model type -> this is FIXED, i.e., given, not estimated
Out.specs{2}
% 1st entry = estimated 1st model parameter (sill)
% 2nd entry = estimated 2nd model parameter (range)
%%
Out.specsdim
% Model specification in 2 dimensions
% 1st entry -> model type (here 2 -> exponential)
% Type vmodtype for a list of availabel models
vmodtype
% Function vmodelspecs translates the entries of Out.ModSpecsDim to 
% those of Out.ModSpecs
vmodelspecs(Out.specsdim)

%% 2D Gaussian simulation of attribute fields with a given variogram model
% You will now generate simulated attribute realizations from a Gaussian 
% spatial process with a spatially constant mean (no first-order effects)
% and the above variogram model as its spatial smoothness descriptor 
% (2nd-order stationarity). Realizations will be generated at the grid 
% nodes of a (300 x 360) regular raster, whose specification are stored 
% in array |GridSpecs|, using custom function |mafftsim|. 
% 
% First, let's see the contenct of array |GridSpecs| with the 
% grid specifications:
disp(GridSpecs);

%%
% Next, create a vector |simPars| with simulation parameters as:
simPars = [10 1234 0];
%%
% corresponding to |nSim = 10| realizations with a random number seed of 
% |1234| and random (|0|) sampling. 
%
% Last, specify the target mean of the simulations as that of the 
% sample precipitation data:
targMean = mean(bayprecip(:,4))
%%
% You can now generate the simulated realizations as: 
vModelSpecs = Out.specsdim;
vModelSpecs(:,2) = 1; % we need a stanardized (unit sill) model
Sim = mafftsim(GridSpecs,[],vModelSpecs,targMean,simPars);
%%
% Note that the output from function |mafftsim| contains |(nx*ny)| rows
% and |nSim| columns; each column is a simulated attribute image in this
% case, stored as a long 1D vector not as a 2D matrix.
disp(size(Sim));
disp([GridSpecs(:,1)' prod(GridSpecs(:,1))]);
%%
% You can display, say, the 1st simulated realization as:
rastermap(GridSpecs,Sim,1,'cont',[],{'color','vbar'});
title('Realization #1');
%%
nBins = 50;
hist(Sim(:,1),nBins);
title('Histogram of realization #1');
%%
% and, say, the 8-th realization as:
rastermap(GridSpecs,Sim,8,'cont',[],{'color','vbar'});
title('Realization #8');
%%
hist(Sim(:,8),nBins);
title('Histogram of realization #8');

%% 2D realizations with non-Gaussian histograms
DataForHist = bayprecip;
icol = 4;
[min(DataForHist(:,4)) max(DataForHist(:,4))]
minmax = [1 15];
Sim = mafftsim(GridSpecs,[],vModelSpecs,targMean,simPars,DataForHist,4,minmax);
%%
% Note that the output from function |mafftsim| contains |(nx*ny)| rows
% and |nSim| columns; each column is a simulated attribute image in this
% case, stored as a long 1D vector not as a 2D matrix.
disp(size(Sim));
disp([GridSpecs(:,1)' prod(GridSpecs(:,1))]);
%%
% You can display, say, the 1st simulated realization as:
rastermap(GridSpecs,Sim,1,'cont',[],{'color','vbar'});
title('Realization #1');
%%
nBins = 50;
hist(Sim(:,1),nBins);
title('Histogram of realization #1');
figure
hist(bayprecip(:,4),30);
title('Data histogram');
%%
% and, say, the 8-th realization as:
rastermap(GridSpecs,Sim,8,'cont',[],{'color','vbar'});
title('Realization #8');
%%
hist(Sim(:,8),nBins);
title('Histogram of realization #8');