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
Bandwidth selection for kernel smoothing of a cycle matrix. | |
Display message and abort function. | |
X and Y arrays for 3-D plots. | |
Extract upper triangular part. |
Calculates the cycle count matrix from a cycle count. | |
Script to computer exercises 2 | |
Script to computer exercises 4 | |
Rainflow matrix for Switching Markov Chains of Turning Points. | |
Extrapolates a rainflow matrix. | |
Quick test of the routines in module 'cycles' |
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
Comments or corrections to the WAFO group