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

Contents 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 - (show 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 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