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:
 kdefun Kernel Density Estimator. kdeoptset Create or alter KDE OPTIONS structure. mkernelrnd Random numbers from the Multivariate Kernel Function. tranproc Transforms process X and up to four derivatives class Create object or return object class. cov Covariance matrix. deal Deal inputs to outputs. det Determinant. error Display message and abort function. iscell True for cell array. isequal True if arrays are numerically equal. lower Convert string to lowercase. mean Average or mean value. num2cell Convert numeric array into cell array. sqrtm Matrix square root. x
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