/[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.1 - (hide annotations) (download)
Thu Jan 6 23:16:55 2011 UTC (14 years, 6 months ago) by gforget
Branch: MAIN
- new directory dedicated to statistical analyses of MITprof files.

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     if nargin>4; choiceStat=varargin{1}; else; choiceStat='pdf'; end;
18     if nargin>5; z=varargin{2}; else; z=[]; end;
19    
20     %take care of input dimensions:
21     if size(x,1)==1; x=x'; end; if size(y,1)==1; y=y'; end;
22     %
23     if size(y,2)==1; %x is array => make sure y is so
24     if size(x,2)==size(y,1);
25     y=ones(size(x,1),1)*y';
26     else;
27     y=y*ones(1,size(x,2));
28     end;
29     end;
30     %
31     if size(x,2)==1; %y is array => make sure x is so
32     if size(y,2)==size(x,1);
33     x=ones(size(y,1),1)*x';
34     else;
35     x=x*ones(1,size(y,2));
36     end;
37     end;
38     %
39     if size(z,2)>1;%z is array => make sure x/y are so
40     if size(x,2)==1&size(x,1)==size(z,1); x=x*ones(1,size(z,2)); end;
41     if size(x,2)==1&size(x,1)==size(z,2); x=ones(size(z,1),1)*x'; end;
42     %
43     if size(y,2)==1&size(y,1)==size(z,1); y=y*ones(1,size(z,2)); end;
44     if size(y,2)==1&size(y,1)==size(z,2); y=ones(size(z,1),1)*y'; end;
45     end;
46     if isempty(z); z=NaN*x; end;%to facilitate vector length tests
47     %switch to vector form:
48     x=x(:); y=y(:); z=z(:);
49     %check vectors length:
50     if length(x)~=length(y)|length(x)~=length(z)|length(y)~=length(z);
51     error('inconsistent diemnsions of [x,y[,z]]');
52     end;
53    
54     %prepare output arrays:
55     if size(xBinCenters,1)==1; xh=xBinCenters'; else; xh=xBinCenters; end;
56     if size(yBinCenters,1)==1; yh=yBinCenters'; else; yh=yBinCenters; end;
57     if size(xh,2)>1|size(yh,2)>1; error('inconsistent diemnsions of [xBinCenters yBinCenters]'); end;
58     xh=xh*ones(1,size(yh,1));
59     yh=ones(size(xh,1),1)*yh';
60    
61     %build the histogram
62     TRI=delaunay(xh,yh); nxy = prod(size(xh));
63     Stri=sparse(TRI(:,[1 1 2 2 3 3]),TRI(:,[2 3 1 3 1 2]),1,nxy,nxy);
64    
65     %compute grid point vector associated with lon/lat vectors
66     if nargin>5; ii=find(~isnan(x.*y.*z)); else; ii=find(~isnan(x.*y)); end;
67     x=x(ii); y=y(ii); z=z(ii);
68     ik=dsearch(xh,yh,TRI,x,y,Stri);
69    
70     %compute bin sums:
71     nh=zeros(size(xh)); if nargin>5; zh=nh; zh2=zh; end;
72     for k=1:length(ik)
73     nh(ik(k))=nh(ik(k))+1;
74     if nargin>5;
75     zh(ik(k))=zh(ik(k))+z(k);
76     zh2(ik(k))=zh2(ik(k))+z(k)^2;
77     end;
78     end % k=1:length(ik)
79     ii=find(nh(:)==0); nh(ii)=NaN;
80     if nargin>5; zh(ii)=NaN; zh2(ii)=NaN; end;
81    
82     %compute output:
83     if strcmp(choiceStat,'pdf'); %hist + normalization along second dimension
84     zh=nh; zh(isnan(zh))=0;
85     %normalize by total number of samples
86     zh=zh./(nansum(zh,2)*ones(1,size(zh,2)));
87     %normalize by delta(yh) => pdf
88     tmp2=diff(yh,1,2);
89     tmp2=0.5*(tmp2(:,1:end-1)+tmp2(:,2:end));
90     tmp2=[tmp2(:,1) tmp2 tmp2(:,end)];
91     zh=zh./tmp2;
92     elseif strcmp(choiceStat,'hist');
93     zh=nh;
94     elseif strcmp(choiceStat,'sum');
95     zh=zh;
96     elseif strcmp(choiceStat,'mean');
97     zh=zh./nh;
98     elseif strcmp(choiceStat,'var');
99     zh=zh2./nh-(zh./nh).^2;
100     end;
101    
102     if nargout>3;
103     varargout{1}=nh;
104     else;
105     varargout={};
106     end;
107    
108     %for debugging purposes:
109     if 0; figure; pcolor(xh,yh,zh); colorbar; end;
110    
111     warning('on','MATLAB:dsearch:DeprecatedFunction');
112    

  ViewVC Help
Powered by ViewVC 1.1.22