Wolfgang Schwanghart

@Sébastien: Thanks for pointing at the error when using a nugget variogram.

I am currently abroad, so I probably won't find much time to update the function. So, in the meantime please adopt the changes proposed by Sébastien.

The problem with the gaussian variogram is most likely due to negative weights (see http://dx.doi.org/10.1016/0098-3004(96)00005-2 ). I will try to implement the procedure suggested by Deutsch (1996) in a future update.

Jiro, thanks for the quick reply and fix!

@Wolfgang: I really appreciate your interactivity and your pertinence in all your comments. I just coded the Deutsch method to solve the problem of the negative weights. Thanks to your doi link. Thus, the gaussian variogram is now well working :)

Here are the lines code to insert within the "for" loop, between the paragraphs "solve system" and "estimate zi".

However, I'm not an excellent developper. I'll let you find some better/faster way to code this.

% solve system
lambda = A*b;

% correct negative weights (Deutsch 1996)
if sum(lambda(:)<0) % if negative weigths are found
% compute the "average absolute magnitude of the negative weights"
lambdaLim=lambda(1:end-1,:);
PosInd=lambdaLim>=0;
lambdaLim(PosInd)=NaN;
lambdaLim=nanmean(abs(lambdaLim),1);
% compute the mean covariance between "negative" locations and the IX ones
bNoOnes=b(1:end-1,:);
bLim=bNoOnes;
bLim(PosInd)=NaN;
bLim=nanmean(bLim,1);
% correct the weights
lambdaC=lambda(1:end-1,:); % step 1 of Deutsch 1996
lambdaC(lambdaC<0)=0; % step 2 of Deutsch 1996
lambdaC(PosInd & bNoOnes-repmat(bLim,numobs,1)<0 & lambda(1:end-1,:)-repmat(lambdaLim,numobs,1)<0)=0; % step 3 of Deutsch 1996
% restandardization of the weights
lambda(1:end-1,:)=lambdaC./repmat(sum(lambdaC,1),numobs,1);
end

% estimate zi
zi(IX) = lambda'*z;