function [nukSmth,nukSmthU,nukSmthL,nukFilt,nukFiltU,nukFiltL] = NoncentChi2FiltSmooth(Devk,Md,sz0,rho,NN,Nem)

% This function is a simplified version of what is available in the AGC
% toolbox that has been modified for the deviance tests for coordinated
% spiking. URL: https://github.com/Arsha89/AGC_Analysis

    K = length(Devk);
    nukFilt = zeros(K,1);
    nukFiltU = zeros(K,1); nukFiltL = zeros(K,1);
    nukSmth = zeros(K,1);
    nukSmthU = zeros(K,1); nukSmthL = zeros(K,1);
    eta = norminv(0.975); % the 95% confidence intervals

        Dnk = Devk;
        zk = zeros(K,1); s2k = zeros(K,1); 
        nuk = zeros(K,1); nukU = zeros(K,1); nukL = zeros(K,1);
        sz = sz0;

        % EM Algorithm on sz parameter
        for ll = 1:Nem
            zk(1) = 0; 
            s2k(1) = 1;
            for k = 1:K
                zl = rho*zk(k);
                for l = 1:NN
                    nul = exp(zl);
                    ul2 = nul*max(Dnk(k),0);
                    fl = sqrt(ul2 + (Md-1)^2/4) - (Md-1)/2;
                    gl = zl - rho*zk(k) + (rho^2*s2k(k)+sz)*(nul - fl)/2;
                    gpl = 1 + (rho^2*s2k(k)+sz)*(nul - ul2/(2*(fl+(Md-1)/2)))/2;
                    zl = zl - gl/gpl;
                end
                zk(k+1) = zl;
                nuk(k) = exp(zk(k+1));
                uzk = sqrt(nuk(k)*max(Dnk(k),0));
                fuk = sqrt(uzk^2 + (Md-1)^2/4) - (Md-1)/2;
                s2k(k+1) = 1/( 1/(rho^2*s2k(k) + sz) + nuk(k)/2 - (uzk^2 - (Md-2)*fuk - fuk^2)/4 );
                nukU(k+1) = exp(zk(k+1)+ eta*sqrt(s2k(k+1)));
                nukL(k+1) = exp(zk(k+1)- eta*sqrt(s2k(k+1)));
            end

            %%% Smoothing the Deviance Samples
            zkK = zk; s2kK = s2k;
            s2kp1K = zeros(K,1);
            for k = K:-1:1
                sk = rho*s2k(k+1)/(rho^2*s2k(k+1)+ sz);
                zkK(k) = zk(k+1)+ sk*(zkK(k+1) - rho*zk(k+1));
                s2kK(k) = s2k(k+1)+ (sk^2)*(s2kK(k+1) - (rho^2*s2k(k+1)+sz));
                s2kp1K(k) = sk*s2kK(k+1);
            end
            
            nukKU = nukU; nukKL = nukL;
            for k = 1:K+1
                nukKU(k) = exp(zkK(k)+ eta*sqrt(s2kK(k)));
                nukKL(k) = exp(zkK(k)- eta*sqrt(s2kK(k)));
            end

            alpha = -1; beta = 0;
            % EM update rule for sz based on smoothed values
            Eterm = rho^2*zkK(1)^2 + (1+rho^2)*sum(zkK(2:K).^2)+ zkK(K+1).^2 + rho^2*s2kK(1)+ (1+rho^2)*sum(s2kK(2:K)) + s2kK(K+1) -2*rho*(sum(s2kp1K)+ sum(zkK(2:K+1).*zkK(1:K)));
            szp =  (Eterm/2+beta)/(K/2+alpha+1);
            sz = szp;
        end
        nukFilt(:) = nuk;
        nukFiltU(:) = nukU(2:end); nukFiltL(:) = nukL(2:end);
        nukSmth(:) = exp(zkK(2:end));
        nukSmthU(:) = nukKU(2:end); nukSmthL(:) = nukKL(2:end);

end