%TLASSO_PATH Compute the maximum likelihood estimates for a t-regression
% Users should refer to the 'tlasso' function to use this software.
function [b, b0, sigma2, L, iter] = tlasso_ML(X, y, nu)

%% Setup
seps = sqrt(eps);
convcrit = 1e-6;
[~, p] = size(X);
n = length(y);
X = [ones(n,1) X];

%% Initial parameter estimates
k = 1;
beta_old = zeros(p + 1, 1);
b = X\y;
b = b(2:end);
b0 = 0;
sigma2 = 1;

%% EM algorithm
done = false;
iter = 0;
while(~done)
    
    %% E-step
    mu = b0 + X(:,2:end)*b;
    e2 = (y - mu).^2 / sigma2;
    w = (nu + 1) ./ (nu + e2);
    
    %% M-step
    beta = lscov(X, y, w);
    b0 = beta(1);
    b = beta(2:end);
    sigma2 = sum(w .* (y - mu).^2) / n;    

    %% Compute negative log-likelihood
    mu = b0 + X(:,2:end)*b;
    e2 = (y - mu).^2 / sigma2;    
    L = n/2*log(sigma2)  - n*gammaln((nu+k)/2) + n*k*gammaln(0.5) + n*gammaln(nu/2) + n*k/2*log(nu);
    L = L + (nu + k)/2 * sum(log(1 + e2/nu));
    
    %% Check convergence
    if(~any(abs(beta - beta_old) > convcrit * max(seps, abs(beta_old))))
        done = true;
    end    
    
    beta_old = beta;    
    iter = iter + 1;
end


end