netcdf/ncdiffvars.m

   1 function [relerr,ssq,p]=ncdiffvars(s,v1,v2)
   2 % [relerr,ssq,p]=ncdiffvars(s,v1,v2)
   3 % compare two array values v1 and v2 using a description s
   4 % relerr = max relative difference
   5 % ssq = square mean relative difference
   6 % p = p-value, close to +1 -1 means errors significantly in one direction
   7 % for rounding errors expect relerr=const*eps, p small
   8 scale=max(big(v1),big(v2));
   9 relerr=big(v2(:)-v1(:))/(scale+realmin);
  10 ssq=norm(v2(:)-v1(:))/(max(norm(v1(:)),norm(v2(:)))+realmin);
  11 d=(v2(:)-v1(:))/(scale+realmin); % scaled diff
  12 avgdiff=mean(d);
  13 stdev=std(d);
  14 n=length(v1(:));
  15 t=sqrt(n)*avgdiff/(stdev+realmin);
  16 p=erf(t);
  17
  18 fprintf('%s max abs %g relative diff %g max %g stdev %g avg %g t-stats %g p-value %g\n',...
  19     s,scale,relerr,max(d),stdev,avgdiff,t,p)
  20 end