function [theta,sigma2] = EMt0(id,t,x,a,b,theta,sigma2,nk,order,nu,maxit)
%function [theta,sigma2] = EMt0(id,t,x,a,b,theta,sigma2,nk,order,nu,maxit)
%Semiparametric t-estimators of the mean (via EM algorithm)
%
%NOTE: since the time grids may be sparse and irregular, the data is input 
% as a concatenated vector rather than a matrix; a vector of labels is used
% to identify the individuals
%
%INPUT:
%   id:     Individual labels           (N x 1)
%   t:      Time grid                   (N x 1)
%   x:      Observations                (N x 1)
%   a:      Time range, lower bound     (scalar)
%   b:      Time range, upper bound     (scalar)
%   theta:  Initial spline coefficients ((nk+order) x 1)
%             (Enter [] for default)
%   sigma2: Initial variance estimator  (scalar)
%             (Enter [] for default)
%   nk:     Number of knots             (scalar)
%             (Equispaced knots in [a,b] will be used)
%   order:  Spline order                (scalar)
%   nu:     t model degrees of freedom  (scalar)
%             (Use nu = 1 for Cauchy estimators)
%   maxit:  Max. number of iterations   (scalar)
%
%OUTPUT:
%   theta:  Estimated spline coefficients ((nk+order) x 1)
%   sigma2: Estimated variance            (scalar)
%
%External programs called: BSPL

indiv = unique(id);
n = length(indiv);

knots = linspace(a,b,nk+2);
p = nk+order;

if isempty(theta)
    theta = zeros(p,1);
end
if isempty(sigma2)
    sigma2 = var(x);
end

err = 1;
iter = 0;
while err>1e-3 && iter<maxit
    sigma20 = sigma2;
    iter = iter + 1;
    A1 = zeros(p,p);
    A2 = zeros(p,1);
    ss = 0;
    N = 0;
    for i = 1:n
        tt = t(id==indiv(i));
        xx = x(id==indiv(i));
        m = length(tt);
        B = bspl(tt,order,knots,0);
        d2 = (xx-B*theta)'*(xx-B*theta)/sigma2;
        Eu = (nu+m)/(nu+d2);
        A1 = A1 + Eu*B'*B;
        A2 = A2 + Eu*B'*xx;
        ss = ss + Eu*norm(xx-B*theta)^2;
        N = N + m;
    end
    theta = inv(A1)*A2;
    sigma2 = ss/N;
    err = abs(sigma2/sigma20-1);
%    disp(['Iter: ' num2str(iter) ', S^2 = ' num2str(sigma2) ', Error = ' num2str(err)])
end