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')
Kernel Density Estimator. | |
Create or alter KDE OPTIONS structure. | |
Random numbers from the Multivariate Kernel Function. | |
Transforms process X and up to four derivatives | |
Create object or return object class. | |
Covariance matrix. | |
Deal inputs to outputs. | |
Determinant. | |
Display message and abort function. | |
True for cell array. | |
True if arrays are numerically equal. | |
Convert string to lowercase. | |
Average or mean value. | |
Convert numeric array into cell array. | |
Matrix square root. | |
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
Comments or corrections to the WAFO group