function [ b ] = Threshold_ERgraph( Name_Procedure, T, Num_vertice, Pre_p, Gamma, w )
%THRESHOLD_ERGRAPH Summary of this function goes here
%   Detailed explanation goes here

mc = 500;

d = round(Num_vertice*(Num_vertice-1)/2);

epsilon = 0.01;

switch Name_Procedure
    % CUSUM
    case 'CUSUM'
        % To construct the CUSUM procedure we need an exactly specified
        % post-change mean vector
        Post_p_CUSUM = 0.5;
        % Maximal statistics of the trials, used to compute b by setting b
        % as a quantile
        stat_max = zeros(mc,1);
        for i = 1 : mc
            
            if (mod(i,100) == 0)
                disp(i);
            end
            % reset the statistics
            stat = zeros(T,1); stat0 = 0;
            
            % generate pre-change sequence with length T
            Signal = binornd(1,Pre_p,d,T);
            
            % initialization: compute stat(1)
            theta0 = log(Pre_p/(1-Pre_p)) * ones(d,1);
            theta1 = log(Post_p_CUSUM/(1-Post_p_CUSUM)) * ones(d,1);
            
            x = Signal(:,1);
            stat(1) = max(stat0 + NLL_bernoulli(x, theta0) - ...
                NLL_bernoulli(x, theta1), -1e08);
            for t = 2 : T
                % Acquire a new observation
                x = Signal(:,t);
                
                % compute stat(t) from stat(2)
                stat(t) = max(stat(t-1) + NLL_bernoulli(x, theta0) - ...
                    NLL_bernoulli(x, theta1), -1e08);
            end
            stat_max(i) = max(stat);
            
        end
        b = quantile(stat_max, exp(-T/Gamma));

    case 'WLGLR'
        % Maximal statistics of the trials, used to compute b by setting b
        % as a quantile
        stat_max = zeros(mc,1);
        
        
        for i = 1 : mc
            if (mod(i,100) == 0)
                disp(i);
            end
            % reset the statistics
            stat = -1e08*ones(T,1);
            theta0 = log(Pre_p/(1-Pre_p))*ones(d,1);
            mu = zeros(d,w);
            S = -1e08*ones(1,w); eta = zeros(1,w);
            PartialSum = zeros(T,1);
           
            % generate pre-change sequence with length T
            Signal = binornd(1,Pre_p,d,T);
            for t = 1 : T
                % Acquire a new observation
                x = Signal(:,t);
                if t == 1
                    PartialSum(t) = NLL_bernoulli(x,theta0);
                else
                    PartialSum(t) = PartialSum(t-1) + NLL_bernoulli(x,theta0);
                end
                
                % Using the idea of circle array
                index = mod(t-1,w)+1;
                % reset the "index"th column
                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;
                    eta(1,index) = eta(1,index)+1;
                    mu(:,index) = mu(:,index) - ...
                        1/eta(1,index) * (mu(:,index) - x);
                    % projection to avoid log(0) and log(inf)
                    %mu(:,index) = max(epsilon, mu(:,index)); 
                    %mu(:,index) = min(1-epsilon,mu(:,index));
                    
                    
                    %Compute t*Phi_star
                    Phi_star = 0;
                    for l = 1 : d
                        n = mu(l,index);
                        if (n ~= 0) && (n ~= 1)
                            Phi_star = Phi_star + ...
                                n*log(n) + (1-n)*log(1-n);
                        end
                    end
                    
                    if j == 1
                        S(1,index) = PartialSum(t) + t*Phi_star;
                    else
                        S(1,index) = PartialSum(t) - PartialSum(j) +...
                            (t-j+1)*Phi_star;
                    end
                     
                end
                stat(t) = max(S);
            end
            stat_max(i) = max(stat);
            
        end
        b = quantile(stat_max, exp(-T/Gamma));
        
    % WLGLR_V
    
    case 'WLGLR_V'
        stat_max = zeros(mc,1);
        for i = 1 : mc
            if (mod(i,100) == 0)
                disp(i);
            end        
            % 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
            theta0 = log(Pre_p/(1-Pre_p))*ones(d,1);
            Signal = binornd(1,Pre_p,d,T);
            
            for t = 1 : T
                for k = max(t-w,1) : t
                    MLE = mean(Signal(:,k:t),2);
                    MLE = max(MLE,1e-08); MLE = min(MLE, 1-1e-08);
                    MLE = log(MLE./(1-MLE));
                    S(k) = 0;
                    for l = k : t
                        x = Signal(:,l);
                        S(k) = S(k) + NLL_bernoulli(x, theta0) - ...
                            NLL_bernoulli(x, MLE);
                    end
                end
                stat(t) = max(S(max(t-w,1):t));               
            end
            stat_max(i) = max(stat);
        end
        b = quantile(stat_max, exp(-T/Gamma));        
        
    case 'WLACM'
        % Maximal statistics of the trials, used to compute b by setting b
        % as a quantile
        stat_max = zeros(mc,1);
        
        
        for i = 1 : mc
            if (mod(i,100) == 0)
                disp(i);
            end
            % reset the statistics
            stat = -1e08*ones(T,1);
            theta = zeros(d,w);
            theta0 = log(Pre_p/(1-Pre_p))*ones(d,1);
            mu = zeros(d,w);
            S = -1e08*ones(1,w); eta = zeros(1,w);
            
           
            % generate pre-change sequence with length T
            Signal = binornd(1,Pre_p,d,T);
            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_bernoulli(x,theta0) - ...
                        NLL_bernoulli(x,theta(:,index));
                    mu(:,index) = exp(theta(:,index))./(1+exp(theta(:,index)));
                    eta(1,index) = eta(1,index)+1;
                    mu(:,index) = mu(:,index) - ...
                        1/eta(1,index) * (mu(:,index) - x);
                    % projection to avoid log(0) and log(inf)
                    mu(:,index) = max(epsilon, mu(:,index)); 
                    mu(:,index) = min(1-epsilon,mu(:,index));
                    
                    theta(:,index) = log(mu(:,index)./(1-mu(:,index)));
                end
                stat(t) = max(S);
            end
            stat_max(i) = max(stat);
            
        end
        b = quantile(stat_max, exp(-T/Gamma));


    case 'WLSRRS'
        % Maximal statistics of the trials, used to compute b by setting b
        % as a quantile
        stat_max = zeros(mc,1);
        
        
        for i = 1 : mc
            if (mod(i,100) == 0)
                disp(i);
            end
            % reset the statistics
            stat = zeros(T,1);
            theta = zeros(d,w);
            theta0 = log(Pre_p/(1-Pre_p))*ones(d,1);
            mu = zeros(d,w);
            Lambda = zeros(1,w); eta = zeros(1,w);
            
           
            % generate pre-change sequence with length T
            Signal = binornd(1,Pre_p,d,T);
            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; 
                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_bernoulli(x,theta0) - ...
                        NLL_bernoulli(x,theta(:,index)));
                    mu(:,index) = exp(theta(:,index))./(1+exp(theta(:,index)));
                    eta(1,index) = eta(1,index)+1;
                    mu(:,index) = mu(:,index) - ...
                        1/eta(1,index) * (mu(:,index) - x);
                    % projection to avoid log(0) and log(inf)
                    mu(:,index) = max(epsilon, mu(:,index)); 
                    mu(:,index) = min(1-epsilon,mu(:,index));
                    
                    theta(:,index) = log(mu(:,index)./(1-mu(:,index)));
                end
                stat(t) = log(sum(Lambda));
            end
            stat_max(i) = max(stat);
            
        end
        b = quantile(stat_max, exp(-T/Gamma));

end


end

