kalmandetection.m

worm detect,基于传染病模型的蠕虫检测算法

%======================================================================
% This program is used for Code Red/Slammer early detection in our paper in CCS'03.

% In our paper in CCS'03, we only used Kalman filter based on epidemic model (Filter #1 in this program).
% Here we add two additional Kalman filters based on two exponential models (see our submitted journal paper on early detection on my webpage).
% The input file is from Code Red C simulation program. Writing Matlab program for data processing is very convienent and good for figure output.
% Updated Date: 02/21/2004:  
%               Add inline explanation for sharing with other researchers.
% Author: Cliff Changchun Zou.


clear all;
FilterEndTime = 310;  % arbitrary value, just used to show the stabilized estimation before the worm finishes propagation.
%BeginTime=0;
BeginTime=1; 
monitorBit = 20;  % monitored IP space 2^{20}.
p = 1/power(2, 32 - monitorBit);  % monitor 20 bit.  This value is used for "bias correction" (see our paper in CCS'03).
eta = 358; 
%eta = 4000;
%eta = 6000;
%eta = 100;
%eta = 1000;
N=360000; 
%N=720000; 
delta =  eta * p; alpha = eta * N / power(2,32);
scanRate = 358;
%scanRate = 4000;
%scanRate = 6000;
%scanRate = 100;
%scanRate = 1000;
q= 1 - power( 1 - 1/power(2,32-monitorBit), scanRate);

threshold = 29.5 * 2;
alarmTime = 3;  % when monitored Z_t is over alarm threshold for 3 consecutive times, activate the Kalman filter.

fid = fopen( 'codered-monitor20-100run.txt', 'r');   % consider noise. this is a output file of 100 simulation runs.
%fid = fopen( 'slammer-monitor20-100run.txt', 'r');
%fid = fopen( 'codered-monitor6000-100run.txt', 'r');
%fid = fopen( 'codered-monitor100-100run.txt', 'r');
%fid = fopen( 'codered-monitor1000-100run.txt', 'r');
simulationNum = 700;
sampleRun = 1;
Datatemp = fscanf( fid, '%ld', [4,inf]);
fclose(fid);

I_all = zeros(simulationNum, sampleRun); C_all = I_all; Z_all = I_all; 

for i=1:sampleRun,
    I_all(:,i) = Datatemp(2, 1+ (i-1)*simulationNum:i*simulationNum)';
    C_all(:,i) = Datatemp(4, 1+ (i-1)*simulationNum:i*simulationNum)';
    Z_all(:,i) = Datatemp(3, 1+ (i-1)*simulationNum:i*simulationNum)';
end

runNo = 1;
I_t = I_all(:, runNo); Z_t = Z_all(:,runNo);  C_t = C_all(:, runNo); % we can change the runNo to deal with different simulation runs. 
DataNum = length(I_t);
% ---------determine when should we begin to use Kalman filter. (BeginTime)
flag = 0;
for i=1:FilterEndTime,
    if Z_t(i) > threshold,
        flag = flag + 1;
    else
        flag = 0;
    end
    if flag >= alarmTime,
        BeginTime = i
        break;
    end
end
%BeginTime = 50;
Step = 1+(FilterEndTime-BeginTime);  %discrete-time Kalman Filter overall step.

% Using low-pass filter. (see our submitted Journal version paper on early detection)
% OldZ_t = Z_t;
% a = 0.3;
% for i=2:FilterEndTime,
%     Z_t(i) = a * OldZ_t(i) + (1-a)*Z_t(i-1);
%     %c1=5; c2=2-c1; Z_t(i) = ((OldZ_t(i)+OldZ_t(i-1))-c2*Z_t(i-1))/c1;
% end
%figure; hold on; plot(OldZ_t,'r--'); plot(Z_t,'b'); axis([0 250 0 max(OldZ_t(255))]);return; 
%testNum = 250;scale = I_t(testNum)./Z_t(testNum);
%figure; hold on; semilogy(I_t(2:testNum),'b'); semilogy(scale*Z_t(2:testNum),'r--'); return;


%--------- use Z_t to estimate ----------------------------------------

% Filter #1:  Epidemic model.

X = [ ones(Step, 1) zeros(Step, 1)]'; Pt_1 = eye(2);  Pt= Pt_1; H = ones(Step, 2); R = 1; K = X;
Z = Z_t(BeginTime:FilterEndTime);
for t=2:Step,
    H(t,:) = [ Z(t-1) Z(t-1)*Z(t-1)];
    K(:,t) = Pt_1 * H(t,:)' /  ( H(t,:)*Pt_1*H(t,:)' + R );
    Pt     = ( eye(2) - K(:,t)*H(t,:) ) * Pt_1;
    X(:,t) = X(:,t-1) + ( Z(t) - H(t,:)*X(:,t-1) )* K(:,t);
    Pt_1 = Pt;
end
AlphaEpidemic = X(1,2:Step)' - 1;

% Filter #2:  AR exponential model.y_t = exp(alpha) y_{t-1}.

Z=Z_t(BeginTime:FilterEndTime); P = 1;  X = ones(1, Step); R =1;K=X;
for t=2:Step - 1,
    R = 1/(t); %1/(sqrt(t));    %put more weight on later part data.
    H = Z(t-1);
    K = P * H / (H*H*P + R);
    P = ( 1 - K*H)*P;
    X(t) = X(t-1) + ( Z(t) - H*X(t-1) )* K;
end
AlphaAR = X(2:Step)'-1;

% Filter #3: Transformed linear model :  ln y_t = ln C + alpha t.

Z = log(Z_t(BeginTime:FilterEndTime)); P = eye(2); I = eye(2);  X = zeros(2, Step); R =1; K=X;
for t = 2:Step -1,
    R = 1/(t*t);    % put more weight on later part data. The reason is explained in our submitted journal version paper on early detection.
    H = [ 1 t ];
    K = P * H' / (H * P * H' + R);
    P = ( I - K * H ) * P;
    X(:, t) = X(:, t-1) +  ( Z(t) - H*X(:, t-1) ) * K;
end
AlphaLinear = X(2,2:Step)';

% %--------- use C_t after bias correction to estimate. ( use Filter #1: epidemic model for Kalman filter)
% X = [ ones(Step, 1) zeros(Step, 1)]'; Pt_1 = eye(2);  Pt= Pt_1; H = ones(Step, 2);R = 1; K = X;
% Infect = C_t;
% for i=BeginTime:length(C_t)-1,
%     Infect(i) = ( C_t(i+1) - (1-q)*C_t(i) ) / q;
% end
% 
% Ihat = Infect(BeginTime:FilterEndTime);
% %Ihat = C_t(BeginTime:FilterEndTime);  % no bias correction.
% for t=2:Step,
%     H(t,:) = [ Ihat(t-1) Ihat(t-1)*Ihat(t-1)];
%     K(:,t) = Pt_1 * H(t,:)' /  ( H(t,:)*Pt_1*H(t,:)' + R );
%     Pt     = ( eye(2) - K(:,t)*H(t,:) ) * Pt_1;
%     X(:,t) = X(:,t-1) + ( Ihat(t) - H(t,:)*X(:,t-1) )* K(:,t);
%     Pt_1 = Pt;
% end
% AlphaC = X(1,2:Step)' - 1; AlphaC(Step - 2);

%---------------- plot results.
RealAlpha = ones(Step-1, 1) * alpha;
RealN     = ones(Step-1, 1) * N;

% ----------- plot Kalman Filter estimation of \alpha, one for each Kalman filter.
figure; hold on;
xlabel('Time t (minute)'); ylabel('Estimated rate');
plot( BeginTime+1:FilterEndTime, RealAlpha, 'k:');
plot( BeginTime+1:FilterEndTime, AlphaEpidemic,'b' ); 
axis( [BeginTime FilterEndTime-10 0 1.2]);
legend('Real value of \alpha', 'Estimated value of \alpha');

%figure; hold on;
%xlabel('Time t (minute)'); ylabel('Estimated infection rate');
%plot( BeginTime+1:FilterEndTime, RealAlpha, 'k:');
%plot( BeginTime+1:FilterEndTime, AlphaAR,'b' ); 
%axis( [BeginTime FilterEndTime-10 0 0.2]);
%legend('Real value of \alpha', 'Estimated value of \alpha');

%figure; hold on;
%xlabel('Time t (minute)'); ylabel('Estimated infection rate');
%plot( BeginTime+1:FilterEndTime, RealAlpha, 'k:');
%plot( BeginTime+1:FilterEndTime, AlphaLinear,'b' ); 
%axis( [BeginTime FilterEndTime-10 0 0.2]);
%legend('Real value of \alpha', 'Estimated value of \alpha');

%------------ plot estimated population N ( see our paper in CCS'03 for this estimation formula).
EstimateNLinear = power(2,32).*AlphaLinear./eta;
EstimateNEpidemic = power(2,32).*AlphaEpidemic./eta;
EstimateNAR = power(2,32).*AlphaAR./eta;

figure; hold on;
xlabel('Time t (minute)'); ylabel('Estimated population N');
plot( BeginTime+1:FilterEndTime, RealN, 'k:');
plot( BeginTime+1:FilterEndTime, EstimateNLinear,'b' ); 
plot( BeginTime+1:FilterEndTime, EstimateNEpidemic,'r-.' ); 

axis( [BeginTime FilterEndTime-10 0 N*1.8]);
legend('Real population size N', 'Transformed linear model', 'Simple epidemic model');


return;