| 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 |
gcmfaces_global; |
| 18 |
|
| 19 |
if ~isfield(myenv,'useDelaunayTri'); |
| 20 |
myenv.useDelaunayTri=~isempty(which('DelaunayTri')); |
| 21 |
end; |
| 22 |
|
| 23 |
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 |
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 |
|
| 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 |
if myenv.useDelaunayTri; |
| 80 |
ik = mytri.TRI.nearestNeighbor(x,y); |
| 81 |
else; |
| 82 |
ik = dsearch(xh,yh,TRI,x,y,Stri); |
| 83 |
end; |
| 84 |
|
| 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 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph