/[MITgcm]/MITgcm_contrib/gael/profilesMatlabProcessing/profiles_stats/MITprof_stats.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/profilesMatlabProcessing/profiles_stats/MITprof_stats.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.3 - (hide annotations) (download)
Wed Dec 12 00:49:20 2012 UTC (12 years, 7 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint65x, checkpoint65r, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, HEAD
Changes since 1.2: +18 -3 lines
- switch to DelaunayTri for recent matlab versions.

1 gforget 1.1 function [xh,yh,zh,varargout]=MITprof_pdf(x,xBinCenters,y,yBinCenters,varargin);
2     %input:
3     % x/y position vectors or 2D arrays
4     % xBinCenters position bin centers
5     % yBinCenters position bin centers
6     % choiceStat [optional] type of statistic (pdf, hist, sum, mean, or var)
7     % z [optional] data vector or 2D array
8     %
9     %note: all inputs are assumed to have be nan-masked accordingly
10     %
11     %output:
12     % xh/yh are the position arrays for the historgam bin centers
13     % zh is the 2D histogram
14    
15     warning('off','MATLAB:dsearch:DeprecatedFunction');
16    
17 gforget 1.3 gcmfaces_global;
18    
19     if ~isfield(myenv,'useDelaunayTri');
20     myenv.useDelaunayTri=~isempty(which('DelaunayTri'));
21     end;
22    
23 gforget 1.1 if nargin>4; choiceStat=varargin{1}; else; choiceStat='pdf'; end;
24     if nargin>5; z=varargin{2}; else; z=[]; end;
25    
26     %take care of input dimensions:
27     if size(x,1)==1; x=x'; end; if size(y,1)==1; y=y'; end;
28     %
29     if size(y,2)==1; %x is array => make sure y is so
30     if size(x,2)==size(y,1);
31     y=ones(size(x,1),1)*y';
32     else;
33     y=y*ones(1,size(x,2));
34     end;
35     end;
36     %
37     if size(x,2)==1; %y is array => make sure x is so
38     if size(y,2)==size(x,1);
39     x=ones(size(y,1),1)*x';
40     else;
41     x=x*ones(1,size(y,2));
42     end;
43     end;
44     %
45     if size(z,2)>1;%z is array => make sure x/y are so
46     if size(x,2)==1&size(x,1)==size(z,1); x=x*ones(1,size(z,2)); end;
47     if size(x,2)==1&size(x,1)==size(z,2); x=ones(size(z,1),1)*x'; end;
48     %
49     if size(y,2)==1&size(y,1)==size(z,1); y=y*ones(1,size(z,2)); end;
50     if size(y,2)==1&size(y,1)==size(z,2); y=ones(size(z,1),1)*y'; end;
51     end;
52     if isempty(z); z=NaN*x; end;%to facilitate vector length tests
53     %switch to vector form:
54     x=x(:); y=y(:); z=z(:);
55     %check vectors length:
56     if length(x)~=length(y)|length(x)~=length(z)|length(y)~=length(z);
57     error('inconsistent diemnsions of [x,y[,z]]');
58     end;
59    
60     %prepare output arrays:
61     if size(xBinCenters,1)==1; xh=xBinCenters'; else; xh=xBinCenters; end;
62     if size(yBinCenters,1)==1; yh=yBinCenters'; else; yh=yBinCenters; end;
63     if size(xh,2)>1|size(yh,2)>1; error('inconsistent diemnsions of [xBinCenters yBinCenters]'); end;
64     xh=xh*ones(1,size(yh,1));
65     yh=ones(size(xh,1),1)*yh';
66    
67     %build the histogram
68 gforget 1.3 if myenv.useDelaunayTri;
69     %the new way, using DelaunayTri&nearestNeighbor
70     mytri.TRI=DelaunayTri(xh(:),yh(:));
71     else;
72     TRI=delaunay(xh,yh,{'QJ'}); nxy = prod(size(xh));
73     Stri=sparse(TRI(:,[1 1 2 2 3 3]),TRI(:,[2 3 1 3 1 2]),1,nxy,nxy);
74     end;
75 gforget 1.1
76     %compute grid point vector associated with lon/lat vectors
77     if nargin>5; ii=find(~isnan(x.*y.*z)); else; ii=find(~isnan(x.*y)); end;
78     x=x(ii); y=y(ii); z=z(ii);
79 gforget 1.3 if myenv.useDelaunayTri;
80     ik = mytri.TRI.nearestNeighbor(x,y);
81     else;
82     ik = dsearch(xh,yh,TRI,x,y,Stri);
83     end;
84 gforget 1.1
85     %compute bin sums:
86     nh=zeros(size(xh)); if nargin>5; zh=nh; zh2=zh; end;
87     for k=1:length(ik)
88     nh(ik(k))=nh(ik(k))+1;
89     if nargin>5;
90     zh(ik(k))=zh(ik(k))+z(k);
91     zh2(ik(k))=zh2(ik(k))+z(k)^2;
92     end;
93     end % k=1:length(ik)
94     ii=find(nh(:)==0); nh(ii)=NaN;
95     if nargin>5; zh(ii)=NaN; zh2(ii)=NaN; end;
96    
97     %compute output:
98     if strcmp(choiceStat,'pdf'); %hist + normalization along second dimension
99     zh=nh; zh(isnan(zh))=0;
100     %normalize by total number of samples
101     zh=zh./(nansum(zh,2)*ones(1,size(zh,2)));
102     %normalize by delta(yh) => pdf
103     tmp2=diff(yh,1,2);
104     tmp2=0.5*(tmp2(:,1:end-1)+tmp2(:,2:end));
105     tmp2=[tmp2(:,1) tmp2 tmp2(:,end)];
106     zh=zh./tmp2;
107     elseif strcmp(choiceStat,'hist');
108     zh=nh;
109     elseif strcmp(choiceStat,'sum');
110     zh=zh;
111     elseif strcmp(choiceStat,'mean');
112     zh=zh./nh;
113     elseif strcmp(choiceStat,'var');
114     zh=zh2./nh-(zh./nh).^2;
115     end;
116    
117     if nargout>3;
118     varargout{1}=nh;
119     else;
120     varargout={};
121     end;
122    
123     %for debugging purposes:
124     if 0; figure; pcolor(xh,yh,zh); colorbar; end;
125    
126     warning('on','MATLAB:dsearch:DeprecatedFunction');
127    

  ViewVC Help
Powered by ViewVC 1.1.22