pdf2d.m

Stanford的SRB实验室Quantitative Seismic Interpretation的免费MATLAB程序

%PDF2D 2D pdf estimation
%function [cpdftbl,nbin1,nbin2]=pdf2d(trdat,nbin1,nbin2,smfilt)
%Make 2D pdfs by smoothing histogram.
%trdat: 2-column training data [data1, data2]
%nbin1, nbin2: (optional) vectors defining bins for data1, data2.
%             if nbin1, nbin2 are scalars, they define the 
%             number of bins
%cpdftbl: discretized 2-D pdf table, Prob(data1, data1)
%smfilt (optional) smoothing filter parameters (size, std. deviation)
% for Gaussian filter kernel; uses 'fspecial('Gaussian')'.

%written by Tapan Mukerji

trdat(isnan(trdat(:,1)),:)=[]; trdat(isnan(trdat(:,2)),:)=[]; %% remove NaNs
%%% select good default bin size and filter size
optbin=floor(sqrt(length(trdat)/0.03));
if nargin==1, 
nbin1=optbin; nbin2=optbin;
end;

if nargin<4, smfilt=fspecial('gaussian',floor(optbin/3),floor(sqrt(optbin)/2));end;
if nargin==4, smfilt=fspecial('gaussian',smfilt(1),smfilt(2)); end;
%%%%


if prod(size(nbin1))==1, nbin1=linspace(min(trdat(:,1)),max(trdat(:,1)),nbin1); end;  
if prod(size(nbin2))==1, nbin2=linspace(min(trdat(:,2)),max(trdat(:,2)),nbin2); end;  
n1=length(nbin1); n2=length(nbin2); ndat=length(trdat(:,1));


pdftbl=zeros(n1,n2);  %% raw histogram
for k=1:ndat
 i1=sum(trdat(k,1)>=nbin1); i2=sum(trdat(k,2)>=nbin2); 
pdftbl(i1,i2)=pdftbl(i1,i2)+1;
end;

pdftbl=filter2(smfilt,pdftbl);
cpdftbl=pdftbl;
cpdftbl=cpdftbl./sum(sum(cpdftbl));

if nargout==0
  imagesc(nbin2,nbin1,cpdftbl); axis xy; colormap hot;colorbar;
 hold on; contour(nbin2,nbin1,cpdftbl,15,'-w'); 
 plot(trdat(:,2),trdat(:,1),'.c','markersize',1); hold off;
end;