function f = naplasso_est(x,y);
% implements the non-adaptive plasso in linear regression with conjugate prior.
% Inputs:
% x:    design matrix without the column 1s
% y:    response vector
% Output: plasso estimation of coefficient vector and noise level
% References: "The predictive Lasso" by Tran, Nott and Leng.
% Written by Minh-Ngoc Tran (ngocmath@gmail.com) with help of the LARS implementation
% of Karl Skoglund

steps = 100;
phi = 2.85;
a = 0.72;
d = 2.58;
e = 1e-10;
[n,p] = size(x);

x1 = x;
y1 = y;
mx = mean(x);
my = mean(y);
x = center(x);
y = center(y);

xx = [ones(n,1),x];
if n>p
    m = xx'*xx\xx'*y;
    sig2 = norm(y-xx*m)^2/(n-p-1);
else
    m = alasso(center(x),center(y));
    m = [mean(y)-mean(x)*m;m];
    sig2 = norm(y-xx*m)^2/n;
end
m(2:(p+1)) = 0;
% determine matrix V inverse
Vinv = diag([1/var(y),var(x)./phi^2]);
betatilde = inv(xx'*xx+Vinv)*(Vinv*m+xx'*y);
Vhat = inv(Vinv+xx'*xx);
s2 = (a+m'*Vinv*m+y'*y-(Vinv*m+xx'*y)'*Vhat*(Vinv*m+xx'*y))/(n+d-2);

ytilde = xx*betatilde;
%wei = abs(betatilde(2:p+1));
%X=x.*(ones(n,1)*wei');
%[XX,mm,ss]=normalize(X);
yy = center(ytilde);
step = 1/(steps - 1);
b_interpolated = zeros(steps, p);
b = lars(x, yy,'lars');
t = sum(abs(b),2);
s = (t - min(t))/max(t - min(t));
[sm s_idx] = unique(s, 'rows');
b_interpolated = interp1q(s(s_idx), b(s_idx, :), (0:step:1)');
b_original = zeros(size(b_interpolated,1),p+1);
for i = 1:size(b_interpolated,1)
    b = b_interpolated(i,:);
    b = b';
    b = [mean(ytilde)-mean(x)*b;b];
    b_original(i,:) = b';
    Cp(i) = norm(ytilde-xx*b_original(i,:)')^2/sig2+ 2*sum(abs(b_original(i,:))>e);
end
[minvalue,index] = min(Cp);
b = b_original(index,:)';
b(1) = my-mx*b(2:length(b));

% find plasso estimate of noise level

xx = [ones(n,1),x1];
y = y1;
Vinv = diag([1/var(y),var(x)./phi^2]);
betatilde = inv(xx'*xx+Vinv)*(Vinv*m+xx'*y);
Vhat = inv(Vinv+xx'*xx);
s2 = (a+m'*Vinv*m+y'*y-(Vinv*m+xx'*y)'*Vhat*(Vinv*m+xx'*y))/(n+d-2);
ytilde = xx*betatilde;
s = 0;
for i = 1:n
    s = s+xx(i,:)*Vhat*xx(i,:)';
end
sig2plasso = (s2*(n+s)+(ytilde-xx*b)'*(ytilde-xx*b))/n;
obj.beta = b;
obj.sig2 = sig2plasso;
f = obj;

end % end of function
