function [ EDD,mean1,std1 ] = EDD_normal_mean( Name_Procedure, T, d, Pre_mean, Post_mean, b, w)
%EDD_NORMAL_MEAN Summary of this function goes here
%   Compute the simulated EDD for certain procedure

% Number of Monte Carlo trials
mc = 2000;

EDD = zeros(mc,1);

switch Name_Procedure
    % CUSUM
    case 'CUSUM'
        for i = 1 : mc
            % reset the statistics
            stat = zeros(T,1); stat0 = 0;
            
            % generate post-change sequence
            Signal = bsxfun(@plus,randn(d,T),Post_mean);

            % To construct the CUSUM procedure we need an exactly specified
            % post-change mean vector
            Post_mean_CUSUM = 1 * ones(d,1);

            % initialization: compute stat(1) and test if over threshold
            x = Signal(:,1);
            stat(1) = max(stat0 + NLL_normal_mean(x, Pre_mean) - ...
                NLL_normal_mean(x, Post_mean_CUSUM), 0);
            if stat(1) > b
                EDD(i) = 1; continue;
            end
            
            for t = 2 : T
                % Acquire a new observation
                x = Signal(:,t);

                % compute stat(t) from stat(2) and test if over threshold
                stat(t) = max(stat(t-1) + NLL_normal_mean(x, Pre_mean) - ...
                    NLL_normal_mean(x, Post_mean_CUSUM), 0);
                if stat(t)>b
                    break;
                end
            end
            EDD(i) = t;
        end
        
        
    % WLGLR
    case 'WLGLR'
        
        for i = 1 : mc
            
            % reset the statistics
            stat = zeros(T,1);
            S = zeros(1,w); % The statistic inside the max (over window-size)
            PartialSum = zeros(d,T); % sum of the recent t observations
            
            % generate pre-change sequence with length T
            Signal = bsxfun(@plus, randn(d,T), Post_mean);
            
            % compute stat(1)
            x = Signal(:,1); PartialSum(:,1) = x;
            S(1,1) = x'*x/2; stat(1) = S(1,1);
            
            for t = 2 : T
                % Acquire a new observation
                x = Signal(:,t);
                PartialSum(:,t) = PartialSum(:,t-1) + x;
                
                % compute stat(t) from stat(2)
                for j = max(1,t-w+1) : t
                    % Using the idea of circle array
                    index = mod(j-1,w)+1;
                    if j == 1
                        y = PartialSum(:,t); % set PartialSum(:,0)=0;
                    else
                        y = PartialSum(:,t) - PartialSum(:,j-1);
                    end
                    S(1,index) = y'*y / (2*(t-j+1));
                end
                stat(t) = max(S);
                if (stat(t) > b)
                    break;
                end
            end
            EDD(i) = t;
            
        end
        
        
        
        
    %WLACM
    case 'WLACM'       
        for i = 1 : mc
            
            % reset the statistics
            stat = zeros(T,1);
            theta = zeros(d,w);
            theta0 = Pre_mean;
            S = zeros(1,w); eta = zeros(1,w);

            % generate post-change sequence
            Signal = bsxfun(@plus,randn(d,T),Post_mean);
            
            for t = 1 : T
                % Acquire a new observation
                x = Signal(:,t);
                
                % Using the idea of circle array
                index = mod(t-1,w)+1;
                % reset the "index"th column
                theta(:,index) = theta0; S(1,index) = 0; eta(1,index) = 0;
                
                
                % compute stat(t) from stat(2)
                for j = max(1,t-w+1) : t
                    index = mod(j-1,w)+1;
                    S(1,index) = S(1,index) + NLL_normal_mean(x,Pre_mean) - ...
                        NLL_normal_mean(x,theta(:,index));
                    eta(1,index) = eta(1,index)+1;
                    theta(:,index) = theta(:,index) - ...
                        1/eta(1,index) * (theta(:,index) - x);
                end
                stat(t) = max(S);
                if stat(t) > b
                    break;
                end
          end
            EDD(i) = t;            
        end        

    % WLSRRS
    case 'WLSRRS'
        for i = 1 : mc
            
            % reset the statistics
            stat = zeros(T,1);
            theta = zeros(d,w);
            Lambda = zeros(1,w); 
            eta = zeros(1,w); % step size
            
            % generate pre-change sequence with length T
            Signal = bsxfun(@plus,randn(d,T),Post_mean);
            
            for t = 1 : T
                % Acquire a new observation
                x = Signal(:,t);
                
                % Using the idea of circle array
                index = mod(t-1,w)+1;
                % reset the "index"th column
                theta(:,index) = Pre_mean; Lambda(1,index) = 1; eta(1,index) = 0;
                
                
                % compute stat(t) from stat(2)
                for j = max(1,t-w+1) : t
                    index = mod(j-1,w)+1;
                    Lambda(1,index) = Lambda(1,index) * exp(NLL_normal_mean(x,Pre_mean) - ...
                        NLL_normal_mean(x,theta(:,index)));
                    eta(1,index) = eta(1,index)+1;
                    theta(:,index) = theta(:,index) - ...
                        1/eta(1,index) * (theta(:,index) - x);
                end
                stat(t) = log(sum(Lambda));
                if (stat(t)>b)
                    break;
                end
            end
            EDD(i) = t;
        end

        
    %WLACM with projection onto L1-ball
    case 'WLACM_L1'
        s = 5;
        for i = 1 : mc
            
            % reset the statistics
            stat = zeros(T,1);
            theta = zeros(d,w);
            theta0 = Pre_mean;
            S = zeros(1,w); eta = zeros(1,w);

            % generate post-change sequence
            Signal = bsxfun(@plus,randn(d,T),Post_mean);
            
            for t = 1 : T
                % Acquire a new observation
                x = Signal(:,t);
                
                % Using the idea of circle array
                index = mod(t-1,w)+1;
                % reset the "index"th column
                theta(:,index) = theta0; S(1,index) = 0; eta(1,index) = 0;
                
                
                % compute stat(t) from stat(2)
                for j = max(1,t-w+1) : t
                    index = mod(j-1,w)+1;
                    S(1,index) = S(1,index) + NLL_normal_mean(x,Pre_mean) - ...
                        NLL_normal_mean(x,theta(:,index));
                    eta(1,index) = eta(1,index)+1;
                    theta(:,index) = theta(:,index) - ...
                        1/eta(1,index) * (theta(:,index) - x);
                    theta(:,index) = ...
                        perform_l1ball_projection(theta(:,index), s);                    
                end
                stat(t) = max(S);
                if stat(t) > b
                    break;
                end
                   
            end
            EDD(i) = t;
            
        end        

    % WLSRRS with projection onto L1-ball
    case 'WLSRRS_L1'
        s = 5;
        for i = 1 : mc
            
            % reset the statistics
            stat = zeros(T,1);
            theta = zeros(d,w);
            Lambda = zeros(1,w); 
            eta = zeros(1,w); % step size
            
            % generate pre-change sequence with length T
            Signal = bsxfun(@plus,randn(d,T),Post_mean);
            
            for t = 1 : T
                % Acquire a new observation
                x = Signal(:,t);                
                % Using the idea of circle array
                index = mod(t-1,w)+1;
                % reset the "index"th column
                theta(:,index) = Pre_mean; Lambda(1,index) = 1; eta(1,index) = 0;                                
                % compute stat(t) from stat(2)
                for j = max(1,t-w+1) : t
                    index = mod(j-1,w)+1;
                    Lambda(1,index) = Lambda(1,index) * exp(NLL_normal_mean(x,Pre_mean) - ...
                        NLL_normal_mean(x,theta(:,index)));
                    eta(1,index) = eta(1,index)+1;
                    theta(:,index) = theta(:,index) - ...
                        1/eta(1,index) * (theta(:,index) - x);
                    theta(:,index) = ...
                        perform_l1ball_projection(theta(:,index), s);                    
                end
                stat(t) = log(sum(Lambda));
                if (stat(t)>b)
                    break;
                end
            end
            EDD(i) = t;
        end
        
    % Mei
%     case 'Mei'        
%         for i = 1 : mc
%             %linear shrinkage factor a
%             a = 0.22;
%             % set b=c=0
%             %hard-threshold omega
%             omega = 0.35;
%             % reset the statistics
%             stat = zeros(T,1);
%             %S = zeros(1,w); % The statistic inside the max (over window-size)
%             Lambda = zeros(T,T);
%             mu_hat = zeros(d,T,T);
%             %PartialSum = zeros(d,T); % sum of the recent t observations
%             
%             % generate pre-change sequence with length T
%             Signal = bsxfun(@plus, randn(d,T), Post_mean);
%             
%             % compute the shrinkage estimators
%             for t = 1:T-1
%                 for tt = t+1:T
%                     temp = mean(Signal(:,t:tt-1),2);
%                     for k = 1:d
%                         if abs(temp(k)) >= omega
%                             mu_hat(k,t,tt) = a*temp(k);
%                         else
%                             mu_hat(k,t,tt) = 0;
%                         end
%                     end
%                 end
%             end
%             % compute the Lambda statistics
%             Lambda(T,T) = 1;
%             for t = 1:T-1
%                 Lambda(t,t) = 1;
%                 for tt = t+1:T
%                     Lambda(tt,t) = Lambda(tt-1,t)*exp(NLL_normal_mean(Signal(:,tt),Pre_mean) - ...
%                         NLL_normal_mean(Signal(:,tt),mu_hat(:,t,tt)));
%                 end
%             end
%             % compute R(n)
%             for t = 1:T
%                 stat(t) = sum(Lambda(t,1:t));            
% %                 if (stat(t) > b)
% %                     break;
% %                 end
%             end
%             EDD(i) = t;
%         end 
    case 'MEI'
        a = 1;
        omega = 0.1;
        for i = 1 : mc
            % reset the statistics
            stat = zeros(T,1);
            S = zeros(T,1); % The statistic inside the max (over window-size)
            % generate pre-change sequence with length T
            Signal = bsxfun(@plus, randn(d,T), Post_mean);
            
            for t = 1 : T
                for k = max(t-w,1) : t
                    %MLE = 1/mean(Signal(:,k:t),2);
                    temp = mean(Signal(:,k:t),2);
                    mu_hat = zeros(d,1);
                    for p = 1:d
                        if abs(temp(p)) >= omega
                            mu_hat(p) = temp(p);
                        else
                            mu_hat(p) = 0;
                        end
                    end
                    S(k) = 0;
                    for l = k : t
                        x = Signal(:,l);
                        S(k) = S(k) + NLL_normal_mean(x,Pre_mean) - ...
                         NLL_normal_mean(x,mu_hat);
                    end
                end
                stat(t) = max(S(max(t-w,1):t));
                if stat(t) > b
                    break;
                end            
            end
            EDD(i) = t;

        end
        
%         for i = 1 : mc
%             
%             % reset the statistics
%             stat = zeros(T,1);
%             S = zeros(1,w); % The statistic inside the max (over window-size)
%             PartialSum = zeros(d,T); % sum of the recent t observations
%             
%             % generate pre-change sequence with length T
%             Signal = bsxfun(@plus, randn(d,T), Post_mean);
%             
%             % compute stat(1)
%             x = Signal(:,1); PartialSum(:,1) = x;
%             S(1,1) = x'*x/2; stat(1) = S(1,1);
%             
%             for t = 2 : T
%                 % Acquire a new observation
%                 x = Signal(:,t);
%                 PartialSum(:,t) = PartialSum(:,t-1) + x;
%                 
%                 % compute stat(t) from stat(2)
%                 for j = max(1,t-w+1) : t
%                     % Using the idea of circle array
%                     index = mod(j-1,w)+1;
%                     if j == 1
%                         y = PartialSum(:,t); % set PartialSum(:,0)=0;
%                     else
%                         y = PartialSum(:,t) - PartialSum(:,j-1);
%                     end
%                     S(1,index) = y'*y / (2*(t-j+1));
%                 end
%                 stat(t) = max(S);
%                 if (stat(t) > b)
%                     break;
%                 end
%             end
%             EDD(i) = t;
%             
%         end
    end

mean1 = mean(EDD);
std1 = std(EDD);

end

