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.

## CROSS-REFERENCE INFORMATION

This function calls:
 smoothcmat_hnorm Bandwidth selection for kernel smoothing of a cycle matrix. error Display message and abort function. meshgrid X and Y arrays for 3-D plots. triu Extract upper triangular part.
This function is called by:
 cc2cmat Calculates the cycle count matrix from a cycle count. itmkurs_lab2 Script to computer exercises 2 itmkurs_lab4 Script to computer exercises 4 rfcdemo2 Rainflow matrix for Switching Markov Chains of Turning Points. rfmextrapolate Extrapolates a rainflow matrix. test_cycles Quick test of the routines in module 'cycles'

## 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 %
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