condpdfs.m

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

%CONDPDFS Conditional pdfs P(data|class) and P(class|data)
%function [pdftbl,nbin1,nbin2,cpdftbl]=condpdfs(trdat,nbin1,nbin2,smfilt)
%Make conditional pdfs for different facies categories.
%nbin1, nbin2: vectors defining bins for data1, data2.
%             if nbin1, nbin2 are scalars, they define the 
%             number of bins
%trdat: 3-column training data [data1, data2, group]
%       can have any number of groups. If only one group
%       trdat can be a 2-column matrix. 
%pdftbl: discretized data-conditioned pdf table, Prob(group | data1, data2)
%cpdftbl: discretized class-conditioned pdf table, Prob(data1, data1 | group)
%Assumes uniform equally likely prior for each group.
%
%smfilt (optional) smoothing filter parameters (size, std. deviation)
% for Gaussian filter kernel; default [5, 3]. Uses 'fspecial('Gaussian')'.

%written by Tapan Mukerji

trdat(isnan(trdat(:,1)),:)=[]; trdat(isnan(trdat(:,2)),:)=[]; %% remove NaNs

if nargin<4, smfilt=fspecial('gaussian',5,3);end;
if nargin==4, smfilt=fspecial('gaussian',smfilt(1),smfilt(2)); end;

if size(trdat,2)==2, trdat=[trdat, ones(length(trdat),1)]; end;

catgr=unique(trdat(:,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); n3=length(catgr); ndat=length(trdat(:,1));

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

for k=1:n3,pdftbl(:,:,k)=filter2(smfilt,pdftbl(:,:,k));end;
cpdftbl=pdftbl;
for k=1:n3, cpdftbl(:,:,k)=cpdftbl(:,:,k)./sum(sum(cpdftbl(:,:,k))); end;

nrm=sum(pdftbl,3); nrm(nrm==0)=1; 
pdftbl=pdftbl./nrm(:,:,ones(1,n3));

if nargout==0
  for jj=1:n3, 
  subplot(n3,2,2*jj-1), imagesc(nbin2,nbin1,cpdftbl(:,:,jj)); axis xy; colorbar;
  subplot(n3,2,2*jj), imagesc(nbin2,nbin1,pdftbl(:,:,jj)); axis xy; colorbar;
  end;
end;