Home > wafo > cycles > smoothcmat.m

smoothcmat

PURPOSE ^

Smooth a cycle matrix using (adaptive) kernel smoothing

SYNOPSIS ^

[Fsmooth,h] = smthcmat(F,method,h,NOsubzero,alpha)

DESCRIPTION ^

 SMOOTHCMAT Smooth a cycle matrix using (adaptive) kernel smoothing
 
  CALL:  Fsmooth = smoothcmat(F,method);
         Fsmooth = smoothcmat(F,method,[],NOsubzero);
         Fsmooth = smoothcmat(F,2,h,NOsubzero,alpha);
 
  Input:
         F       = Cycle matrix.           [nxn]
         method  = 1: Kernel estimator (constant bandwidth). (Default)
                   2: Adaptiv kernel estimator (local bandwidth). 
         h       = Bandwidth (Optional, Default='automatic choice')
       NOsubzero = Number of subdiagonals that are zero
                   (Optional, Default = 0, only the diagonal is zero)
         alpha   = Parameter for method (2) (Optional, Default=0.5).
                   A number between 0 and 1.
                   alpha=0 implies constant bandwidth (method 1).
                   alpha=1 implies most varying bandwidth.
 
  Output:
  F       = Smoothed cycle matrix.   [nxn]
  h       = Selected bandwidth.
 
  See also  cc2cmat, tp2rfc, tp2mm, dat2tp

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [Fsmooth,h] = smthcmat(F,method,h,NOsubzero,alpha)
002 %SMOOTHCMAT Smooth a cycle matrix using (adaptive) kernel smoothing
003 %
004 % CALL:  Fsmooth = smoothcmat(F,method);
005 %        Fsmooth = smoothcmat(F,method,[],NOsubzero);
006 %        Fsmooth = smoothcmat(F,2,h,NOsubzero,alpha);
007 %
008 % Input:
009 %        F       = Cycle matrix.           [nxn]
010 %        method  = 1: Kernel estimator (constant bandwidth). (Default)
011 %                  2: Adaptiv kernel estimator (local bandwidth). 
012 %        h       = Bandwidth (Optional, Default='automatic choice')
013 %      NOsubzero = Number of subdiagonals that are zero
014 %                  (Optional, Default = 0, only the diagonal is zero)
015 %        alpha   = Parameter for method (2) (Optional, Default=0.5).
016 %                  A number between 0 and 1.
017 %                  alpha=0 implies constant bandwidth (method 1).
018 %                  alpha=1 implies most varying bandwidth.
019 %
020 % Output:
021 % F       = Smoothed cycle matrix.   [nxn]
022 % h       = Selected bandwidth.
023 %
024 % See also  cc2cmat, tp2rfc, tp2mm, dat2tp
025 
026 % Tested on Matlab 5.3
027 %
028 % History:
029 % Revised by PJ 18-May-2000
030 %   Updated help text.
031 % Revised by PJ  01-Nov-1999
032 %   updated for WAFO
033 % Created by PJ (Pär Johannesson) 1997
034 %   from 'Toolbox: Rainflow Cycles for Switching Processes V.1.0'
035 
036 % Check input arguments
037 
038 ni = nargin;
039 no = nargout;
040 error(nargchk(1,5,ni));
041 
042 if ni<2, method=1; end
043 if ni<3, h=[]; end
044 if ni<4, NOsubzero=[]; end
045 if ni<5, alpha=[]; end
046 
047 if method == 1 | method == 2  % Kernel estimator
048   if isempty(h), aut_h=1; else aut_h=0; end
049   if isempty(NOsubzero), NOsubzero=0; end
050 else
051   error('Input argument "method" should be 1 or 2');
052 end
053 
054 if method == 2   % Adaptive Kernel estimator
055   if isempty(alpha), alpha=0.5; end
056 end
057 
058 n = length(F);    % Size of matrix
059 N = sum(sum(F));  % Total number of cycles
060 
061 Fsmooth = zeros(n,n);
062 
063 if method == 1 | method == 2 % Kernel estimator
064 
065   d = 2;   % 2-dim
066   [I,J] = meshgrid(1:n,1:n);
067 
068   % Choosing bandwidth
069   % This choice is optimal if the sample is from a normal distr.
070   % The normal bandwidth usualy oversmooths,
071   % therefore we choose a slightly smaller bandwidth
072   
073   if aut_h == 1
074     h_norm = smoothcmat_hnorm(F,NOsubzero);
075     h = 0.7*h_norm;         % Don't oversmooth 
076     
077     %h0 = N^(-1/(d+4));
078     %FF = F+F';
079     %mean_F = sum(sum(FF).*(1:n))/N/2;
080     %s2 = sum(sum(FF).*((1:n)-mean_F).^2)/N/2;
081     %s = sqrt(s2);       % Mean of std in each direction
082     %h_norm = s*h0;      % Optimal for Normal distr.
083     %h = h_norm;         % Test
084   end
085 
086   % Calculating kernel estimate
087   % Kernel: 2-dim normal density function
088   
089   test=0;
090   if test==0
091     for i = 1:n-1
092       for j = i+1:n
093         if F(i,j) ~= 0
094           F1 = exp(-1/(2*h^2)*((I-i).^2+(J-j).^2));  % Gaussian kernel
095           F1 = F1+F1';                     % Mirror kernel in diagonal
096           F1 = triu(F1,1+NOsubzero);       % Set to zero below and on diagonal
097           F1 = F(i,j) * F1/sum(sum(F1));   % Normalize
098           Fsmooth = Fsmooth+F1;
099         end
100       end
101     end
102   else
103     
104     [II,JJ]=find(F>0);
105     for k = 1:length(II)
106       i=II(k); j=JJ(k);
107       F1 = exp(-1/(2*h^2)*((I-i).^2+(J-j).^2));  % Gaussian kernel
108       F1 = F1+F1';                     % Mirror kernel in diagonal
109       F1 = triu(F1,1+NOsubzero);       % Set to zero below and on diagonal
110       F1 = F(i,j) * F1/sum(sum(F1));   % Normalize
111       Fsmooth = Fsmooth+F1;
112     end
113     
114   end
115   
116   
117 end
118 
119 if method == 2
120   Fpilot = Fsmooth/N;
121   Fsmooth = zeros(n,n);
122   [I1,I2] = find(F>0);
123   logg = 0;
124   for i =1:length(I1)
125     logg = logg + F(I1(i),I2(i)) * log(Fpilot(I1(i),I2(i)));
126   end
127   g = exp(logg/N);
128   lambda = (Fpilot/g).^(-alpha);
129 
130   for i = 1:n-1
131     for j = i+1:n
132       if F(i,j) ~= 0
133         hi = h*lambda(i,j);
134         F1 = exp(-1/(2*hi^2)*((I-i).^2+(J-j).^2));  % Gaussian kernel
135         F1 = F1+F1';                     % Mirror kernel in diagonal
136         F1 = triu(F1,1+NOsubzero);       % Set to zero below and on diagonal
137         F1 = F(i,j) * F1/sum(sum(F1));   % Normalize
138         Fsmooth = Fsmooth+F1;
139       end
140     end
141   end
142 
143 end
144 
145

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