1function [IR, CC, D] = adim(U, UC, IR, CC, arg5); 2% ADIM adaptive information matrix. Estimates the inverse 3% correlation matrix in an adaptive way. 4% 5% [IR, CC] = adim(U, UC [, IR0 [, CC0]]); 6% U input data 7% UC update coefficient 0 < UC << 1 8% IR0 initial information matrix 9% CC0 initial correlation matrix 10% IR information matrix (inverse correlation matrix) 11% CC correlation matrix 12% 13% The algorithm uses the Matrix Inversion Lemma, also known as 14% "Woodbury's identity", to obtain a recursive algorithm. 15% IR*CC - UC*I should be approx. zero. 16% 17% Reference(s): 18% [1] S. Haykin. Adaptive Filter Theory, Prentice Hall, Upper Saddle River, NJ, USA 19% pp. 565-567, Equ. (13.16), 1996. 20 21 22% $Id$ 23% Copyright (C) 1998-2003 by Alois Schloegl <alois.schloegl@gmail.com> 24% 25% This program is free software: you can redistribute it and/or modify 26% it under the terms of the GNU General Public License as published by 27% the Free Software Foundation, either version 3 of the License, or 28% (at your option) any later version. 29% 30% This program is distributed in the hope that it will be useful, 31% but WITHOUT ANY WARRANTY; without even the implied warranty of 32% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 33% GNU General Public License for more details. 34% 35% You should have received a copy of the GNU General Public License 36% along with this program. If not, see <http://www.gnu.org/licenses/>. 37 38 39 40[ur,p] = size(U); 41 42Mode_E = 1; 43%if nargin < 4, 44% m = zeros(size(U)+[0,Mode_E]); 45%end; 46 47if nargin<2, 48 fprintf(2,'Error ADIM: missing update coefficient\n'); 49 return; 50else 51 if ~((UC > 0) & (UC <1)), 52 fprintf(2,'Error ADIM: update coefficient not within range [0,1]\n'); 53 return; 54 end; 55 if UC > 1/p, 56 fprintf(2,'Warning ADIM: update coefficient should be smaller than 1/number_of_dimensions\n'); 57 end; 58end; 59 60if nargin<3, 61 IR = []; 62end; 63if nargin<4, 64 CC = []; 65end; 66if nargin<5, 67 arg5 = 6; 68end; 69if isempty(IR), 70 IR = eye(p+Mode_E); 71end; 72if isempty(CC), 73 CC = eye(p+Mode_E); 74end; 75 76D = zeros(ur,(p+Mode_E)^2); 77%D = zeros(ur,1); 78W = eye(p+Mode_E)*UC/(p+Mode_E); 79W2 = eye(p+Mode_E)*UC*UC/(p+Mode_E); 80 81for k = 1:ur, 82 if ~Mode_E, 83 % w/o mean term 84 d = U(k,:); 85 else 86 % w/ mean term 87 d = [1,U(k,:)]; 88 end; 89 90 if ~any(isnan(d)), 91 CC = (1-UC)*CC + UC*(d'*d); 92 93 %K = (1+UC)*IR*d'/(1+(1+UC)*d*IR*d'); 94 v = IR*d'; 95 %K = v/(1-UC+d*v); 96 97 if arg5==0; % original version according to [1], can become unstable 98 IR = (1+UC)*IR - (1+UC)/(1-UC+d*v)*v*v'; 99 100 elseif arg5==1; % this provides IR*CC == I, this seems to be stable 101 IR = IR - (1+UC)/(1-UC+d*v)*v*v' + sum(diag(IR))*W; 102 103 elseif arg5==6; % DEFAULT: 104 IR = ((1+UC)/2)*(IR+IR') - ((1+UC)/(1-UC+d*v))*v*v'; 105 106 % more variants just for testing, do not use them. 107 elseif arg5==2; 108 IR = IR - (1+UC)/(1-UC+d*v)*v*v' + sum(diag(IR))*W2; 109 110 elseif arg5==3; 111 IR = IR - (1+UC)/(1-UC+d*v)*v*v' + eps*eye(p+Mode_E); 112 113 elseif arg5==4; 114 IR = (1+UC)*IR - (1+UC)/(1-UC+d*v)*v*v' + eps*eye(p+Mode_E); 115 116 elseif arg5==5; 117 IR = IR - (1+UC)/(1-UC+d*v)*v*v' + eps*eye(p+Mode_E); 118 119 end; 120 end; 121 122 %D(k) = det(IR); 123 D(k,:) = IR(:)'; 124end; 125 126IR = IR / UC; 127if Mode_E, 128 IR(1,1) = IR(1,1) - 1; 129end; 130 131 132