%TLASSO_FINDMINTAU2 Determine the minimum regularisation parameter to use
% Users should refer to the 'tlasso' function to use this software.
function [tau2, b0_null, sigma2_null, L_null] = tlasso_FindMinTau2(X, y, nu)

%% Find the (b0,sigma2) 'null' model
options = optimoptions('fminunc','Display', 'off', 'GradObj', 'off', 'Algorithm', 'quasi-newton');
[x_null, L_null] = fminunc(@(Z) tnegll(y, Z(1), exp(Z(2)), nu), [mean(y), log(var(y))], options);

b0_null = x_null(1);
sigma2_null = exp(x_null(2));

%% Calculate the gradients at the null model w.r.t. to the beta's
e = y - b0_null;
g = -(nu+1)/nu/sigma2_null * (e ./ (1+e.^2/nu/sigma2_null) );
g = (g'*X)';

%% Use the gradients to calculate the minimum value of tau^2 needed to zero out all regression coefficients
gmax = max(abs(g));
tau2 = (sqrt(2)/gmax/sqrt(sigma2_null))^2 * 0.95;

end

%% Compute negative log-likelihood
function L = tnegll(y, b0, sigma2, nu)

k = 1;
n = length(y);
mu = b0;
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));

end