Home > wafo > kdetools > ssample.m

ssample

PURPOSE ^

Random sampling from a smoothed empirical distribution

SYNOPSIS ^

s=ssample(A,m,varargin)

DESCRIPTION ^

 SSAMPLE  Random sampling from a smoothed empirical distribution 
  
  CALL: s = ssample(data,m,options) 
    
    s      = sampled selection from data,  size m x D 
    data   = data matrix, size N x D (D = # dimensions) 
    m      = sampling size  
   options = kdeoptions-structure or cellvector of named parameters with 
             corresponding values, see kdeoptset for details.   
   
   SSAMPLE(DATA,M) selects a random sample of M data points from the 
   multivariate data-set in the matrix DATA. 
  
  Example: 
      data=wnormrnd(0,1,500,1); 
      s = ssample(data,100,'kernel','gauss'); 
      f = kdebin(s);   
      pdfplot(f),hold on, 
      x = linspace(-5,5); 
      plot(x,wnormpdf(x),'r')

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function s=ssample(A,m,varargin) 
002 %SSAMPLE  Random sampling from a smoothed empirical distribution 
003 % 
004 % CALL: s = ssample(data,m,options) 
005 %   
006 %   s      = sampled selection from data,  size m x D 
007 %   data   = data matrix, size N x D (D = # dimensions) 
008 %   m      = sampling size  
009 %  options = kdeoptions-structure or cellvector of named parameters with 
010 %            corresponding values, see kdeoptset for details.   
011 %  
012 %  SSAMPLE(DATA,M) selects a random sample of M data points from the 
013 %  multivariate data-set in the matrix DATA. 
014 % 
015 % Example: 
016 %     data=wnormrnd(0,1,500,1); 
017 %     s = ssample(data,100,'kernel','gauss'); 
018 %     f = kdebin(s);   
019 %     pdfplot(f),hold on, 
020 %     x = linspace(-5,5); 
021 %     plot(x,wnormpdf(x),'r')   
022    
023 %  Reference:   
024 %  B. W. Silverman (1986)  
025 % 'Density estimation for statistics and data analysis'   
026 %  Chapman and Hall, pp. 143  
027 %   
028 %  Wand, M. P. and Jones, M. C. (1995)  
029 % 'Density estimation for statistics and data analysis'   
030 %  Chapman and Hall,  
031  
032 % History: 
033 % revised pab dec2003   
034 % Revised by gl 13.07.2000 
035 %    changed ind generation to avoid stats 
036 % by pab 10.12.1999 
037  
038 % TODO % Needs further checking  
039    
040    defaultoptions = kdeoptset('kernel','gauss');   
041 % If just 'defaults' passed in, return the default options in g 
042 if ((nargin==1) & (nargout <= 1) &  isequal(A,'defaults')), 
043   s = defaultoptions; 
044   return 
045 end    
046  
047 error(nargchk(2,inf, nargin)) 
048 [n d]=size(A); 
049  
050 if (nargin<3) 
051   options  = defaultoptions; 
052 else 
053   switch lower(class(varargin{1})) 
054    case {'char','struct'}, 
055     options = kdeoptset(defaultoptions,varargin{:}); 
056    case {'cell'} 
057     
058       options = kdeoptset(defaultoptions,varargin{1}{:}); 
059    otherwise 
060     error('Invalid options') 
061   end 
062 end 
063 kernel   = options.kernel; 
064 h        = options.hs; 
065 alpha    = options.alpha; 
066 L2       = options.L2; 
067 hsMethod = options.hsMethod; 
068 if isempty(h)  
069   h=zeros(1,d); 
070 end 
071  
072 k3 = []; 
073 if isempty(L2) 
074   L2=ones(1,d); % default no transformation 
075   elseif iscell(L2)   % cellarray of non-parametric and parametric transformations 
076   Nl2 = length(L2); 
077   if ~(Nl2==1|Nl2==d), error('Wrong size of L2'), end 
078   [L22{1:d}] = deal(L2{1:min(Nl2,d)}); 
079   L2 = ones(1,d); % default no transformation 
080   for ix=1:d, 
081     if length(L22{ix})>1, 
082       k3=[k3 ix];       % Non-parametric transformation 
083     else  
084      L2(ix) = L22{ix};  % Parameter to the Box-Cox transformation 
085     end 
086   end 
087 elseif length(L2)==1 
088   L2=L2(:,ones(1,d)); 
089 end 
090  
091 amin=min(A); 
092 if any((amin(L2~=1)<=0))  , 
093   error('DATA cannot be negative or zero when L2~=1') 
094 end 
095  
096 %new call 
097 lA = A;   
098  
099 k1=find(L2==0); % logaritmic transformation 
100 if any(k1) 
101   lA(:,k1)=log(A(:,k1)); 
102 end 
103 k2=find(L2~=0 & L2~=1); % power transformation 
104 if any(k2) 
105   lA(:,k2)=sign(L2(ones(n,1),k2)).*A(:,k2).^L2(ones(n,1),k2); 
106 end 
107 for ix = k3, 
108   lA(:,ix) = tranproc(A(:,ix),L22{ix}); 
109   lX(:,ix) = tranproc(X(:,ix),L22{ix}); 
110 end 
111  
112 hsiz=size(h); 
113 if (min(hsiz)==1)|(d==1) 
114   if max(hsiz)==1, 
115     h=h*ones(1,d); 
116   else 
117     h=reshape(h,[1 d]); % make sure it has the correct shape 
118   end 
119  
120   ind=find(h<=0); 
121   if any(ind)      % If no value of h has been specified by the user then  
122     h(ind)=feval(hsMethod,lA(:,ind),kernel); % calculate automatic values.   
123   end 
124   hmat=diag(h); 
125 else 
126   deth=det(h); 
127   if deth<=0, 
128     error('bandwidth matrix h must be positive definit') 
129   end 
130   hmat=h; 
131 end 
132  
133 ma=mean(lA); 
134 sa=cov(lA);  
135  
136 E=mkernelrnd(kernel,m,d); 
137 % ind = unidrnd(n,m,1);   
138 % new generation of random sample 
139 ind = floor(n*rand(m,1))+1; 
140 if alpha>0, 
141   tmp=num2cell(lA,1); 
142   opt1 = kdeoptset('kernel',kernel,'hs',h,'alpha',alpha,'L2',L2); 
143   [f, hs,lambda]= kdefun(lA,opt1,tmp{:}); 
144  E=E.*lambda(ind,ones(1,d)); 
145 end 
146  
147  switch lower(kernel(1:4)) 
148   case 'epa1', sk=eye(d)/5; 
149   case {'norm','gaus'}, sk=eye(d); 
150   end 
151  
152 s=ma(ones(m,1),:)+(lA(ind,:)-ma(ones(m,1),:)+ E*hmat)/sqrtm(eye(d)+hmat^2*sk/sa); 
153  
154 % transforming back 
155 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
156 if any(k1), % L2=0 i.e. logaritmic transformation 
157  s(:,k1)=exp(s(:,k1)); 
158 end 
159 if any(k2) % L2~=0 i.e. power transformation 
160   for ix=k2 
161     s(:,ix)=sign(L2(ix)).*(s(:,ix)).^(1/L2(ix)); 
162   end 
163 end 
164

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003