function [theta,sigma2] = EMnormal0(id,t,x,a,b,nk,order)
%function [theta,sigma2] = EMnormal0(id,t,x,a,b,nk,order)
%Semiparametric Normal-model 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)
%   nk:     Number of knots             (scalar)
%             (Equispaced knots in [a,b] will be used)
%   order:  Spline order                (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;

A1 = zeros(p,p);
A2 = zeros(p,1);
N = 0;
for i = 1:n
    iab = find(id==indiv(i) & t>=a & t<=b);
    tt = t(iab);
    xx = x(iab);
    m = length(tt);
    B = bspl(tt,order,knots,0);
    A1 = A1 + B'*B;
    A2 = A2 + B'*xx;
    N = N + m;
end
theta = inv(A1)*A2;
ss = 0;
for i = 1:n
    iab = find(id==indiv(i) & t>=a & t<=b);
    tt = t(iab);
    xx = x(iab);
    B = bspl(tt,order,knots,0);
    ss = ss + norm(xx-B*theta)^2;
end
sigma2 = ss/N;