function [ b ] = Threshold_Gamma( Name_Procedure, T, d, Pre_beta, Gamma, w,s )
%THRESHOLD_NORMAL_MEAN Summary of this function goes here
%   This function aims to compute the threshold such that the ARL of the
%   procedure is almost Gamma. Here we use the approximation that the
%   threshold b is almost the (exp(-T/Gamma)) quantile of the maximum of
%   the mc Monte Carlo trials

% Number of Monte Carlo trials of running the procedure on a pre-change
% sequence
mc = 100;


switch Name_Procedure
    % CUSUM
    case 'CUSUM'
        % To construct the CUSUM procedure we need an exactly specified
        % post-change mean vector
        Post_beta_CUSUM = 2 * ones(d,1);
        theta0 = -Pre_beta;
        theta1 = -Post_beta_CUSUM;
        % 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 = gamrnd(1,1/Pre_beta, [1 T]);
            
            % initialization: compute stat(1)
            x = Signal(:,1);
            stat(1) = max(stat0 + NLL_Gamma(x, theta0) - ...
                NLL_Gamma(x, theta1), 0);
            for t = 2 : T
                % Acquire a new observation
                x = Signal(:,t);
                
                % compute stat(t) from stat(2)
                stat(t) = max(stat(t-1) + NLL_Gamma(x, theta0) - ...
                    NLL_Gamma(x, theta1), 0);
            end
            stat_max(i) = max(stat);
            
        end
        b = quantile(stat_max, exp(-T/Gamma));
        
    % WLGLR method without any programming optimization
    case 'WLGLR_V'
        stat_max = zeros(mc,1);
        theta_0 = -Pre_beta;
        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 = gamrnd(1,1/Pre_beta, [1 T]);
            
            for t = 1 : T
                for k = max(t-w,1) : t
                    MLE = -1/mean(Signal(:,k:t),2);
                    S(k) = 0;
                    for l = k : t
                        x = Signal(:,l);
                        S(k) = S(k) + NLL_Gamma(x, theta_0) - ...
                            NLL_Gamma(x, MLE);
                    end
                end
                stat(t) = max(S(max(t-w,1):t));
                    
            end
            stat_max = max(stat);

        end
        b = quantile(stat_max, exp(-T/Gamma));         
        
    % WLACM
    case 'WLACM'
        % Maximal statistics of the trials, used to compute b by setting b
        % as a quantile
        stat_max = zeros(mc,1);
        theta0 = -Pre_beta;
        for i = 1 : mc
            
            if (mod(i,100) == 0)
                disp(i);
            end
            % reset the statistics
            stat = zeros(T,1);
            theta = zeros(d,w);
            S = zeros(1,w); eta = zeros(1,w);
            
            % generate pre-change sequence with length T
            Signal = gamrnd(1,1/Pre_beta, [1 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_Gamma(x,theta0) - ...
                        NLL_Gamma(x,theta(:,index));
                    eta(1,index) = eta(1,index)+1;
                    mu = -1./theta(:,index);
                    mu = mu - 1/eta(1,index) * (mu - x);
                    theta(:,index) = -1./mu;
                end
                stat(t) = max(S);
            end
            stat_max(i) = max(stat);
            
        end
        b = quantile(stat_max, exp(-T/Gamma));
        
    % WLSRRS
    case 'WLSRRS'
        % Maximal statistics of the trials, used to compute b by setting b
        % as a quantile
        theta0 = -Pre_beta;        
        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);
            Lambda = zeros(1,w); 
            eta = zeros(1,w); % step size
            
            % generate pre-change sequence with length T
            Signal = gamrnd(1,1/Pre_beta, [1 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_Gamma(x,theta0) - ...
                        NLL_Gamma(x,theta(:,index)));
                    eta(1,index) = eta(1,index)+1;
                    mu = -1./theta(:,index);
                    mu = mu - 1/eta(1,index) * (mu - x);
                    theta(:,index) = -1./mu;
                end
                stat(t) = log(sum(Lambda));
            end
            stat_max(i) = max(stat);
            
        end
        b = quantile(stat_max, exp(-T/Gamma)); 
        
    % LORDEN 
    case 'LORDEN'
        stat_max = zeros(mc,1);
        theta_0 = -Pre_beta;
        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 = gamrnd(1,1/Pre_beta, [1 T]);
            
            for t = 1 : T
                for k = max(t-w,1) : t
                    MOM = length(Signal(:,k:t))*mean(Signal(:,k:t),2)/...
                        var(Signal(:,k:t))/(length(Signal(:,k:t))-1);
                    S(k) = 0;
                    for l = k : t
                        x = Signal(:,l);
                        S(k) = S(k) + NLL_Gamma(x, theta_0) - ...
                            NLL_Gamma(x, -MOM);
                    end
                end
                stat(t) = max(S(max(t-w,1):t));
                    
            end
            stat_max = max(stat);

        end
        b = quantile(stat_max, exp(-T/Gamma));        
end



end

