function [theta,etas,lambdas,sigma2,y,xhat] = EMnormal(id,t,x,a,b,theta,etas,lambdas,sigma2,nk,order,maxiter)
%function [theta,etas,lambdas,sigma2,y,xhat] = EMnormal(id,t,x,a,b,theta,etas,lambdas,sigma2,nk,order,maxiter)
%Semiparametric Normal-model estimators of the mean and principal components (via EM algorithm)
%
%NOTE 1: 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
%
%NOTE 2: although the program accepts arbitrary initial estimators, it is 
% best used in a stepwise manner, i.e.: fit a mean + error model first 
% (using program EMnormal0), then fit a one-component model using the output THETA
% from EMnormal0 as initial THETA, and so on.
%
%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 mean coefficients   ((nk+order) x 1)
%            (See Note 2 above)
%   etas:   Initial PC coefficients     ((nk+order) x q)
%            (See Note 2 above; we suggest you take etas = [etas0, ones((nk+order),1)]
%            where etas0 is the output from the previous (q-1) component model)
%   lambdas: Initial eigenvalues        (q x 1)
%            (See Note 2 above; we suggest you take lambdas = [lambdas0; lambdas0(end)/2]
%            where lambdas0 is the output from the previous (q-1) component model)
%   sigma2: Initial variance estimator  (scalar)
%             (See Note 2 above)
%   nk:     Number of knots             (scalar)
%             (Equispaced knots in [a,b] will be used)
%   order:  Spline order                (scalar)
%   maxiter: Max. number of iterations  (scalar)
%
%OUTPUT:
%   theta:  Estimated mean coefficients ((nk+order) x 1)
%   etas:   Estimated PC coefficients   ((nk+order) x 1)
%   lambdas: Estimated eigenvalues      (q x 1)
%   sigma2: Estimated variance          (scalar)
%   y:      Predicted component scores  (q x n)
%   xhat:   Predicted values of X       (N x 1)
%
%External programs called: BSPL

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

knots = linspace(a,b,nk+2);
p = nk+order;
q = length(lambdas);
if q==0
    disp('For a mean + error model use EMt0')
    return
end
if size(etas,1)~=p || size(etas,2)~=q
    disp('Wrong dimension of ETAS')
    return
end
if length(theta)~=p
    disp('Wrong dimension of THETA')
    return
end

H = etas*diag(sqrt(lambdas));
grN = max(200,10*nk);
B0 = bspl(linspace(a,b,grN+1),order,knots,0);
J0 = B0'*B0*(b-a)/grN;
[U,S,V] = svd(H'*J0*H);
H = H*U;

y = zeros(q,n);
xhat = NaN(size(x));
err = 1;
iter = 0;
while err>1e-4 && iter<maxiter
    sigma20 = sigma2;
    iter = iter + 1;
    A1 = zeros(p,p);
    A2 = zeros(p,1);
    A3 = zeros(p*q,p*q);
    A4 = zeros(p*q,1);
    ss = 0;
    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);
        V0_1 = inv(eye(q) + H'*B'*B*H/sigma2);
        S0_1 = eye(m)/sigma2 - B*H*V0_1*H'*B'/sigma2^2;
        y(:,i) = H'*B'*S0_1*(xx-B*theta);
        Eu = 1;     % Eu = (nu+m)/(nu+d2) for t(nu)
        Ezz = V0_1 + Eu*y(:,i)*y(:,i)';
        Euz = Eu*y(:,i);
        A1 = A1 + Eu*B'*B;
        A2 = A2 + Eu*B'*(xx-B*H*y(:,i));
        A3 = A3 + kron(Ezz,B'*B);
        A4 = A4 + Eu*kron(y(:,i),B')*(xx-B*theta);
        ss = ss + Eu*norm(xx-B*theta)^2 + trace(B*H*Ezz*H'*B') ...
            - 2*Euz'*H'*B'*(xx-B*theta);
        N = N + m;        
        xhat(iab) = B*theta + B*H*y(:,i);
    end
    theta = inv(A1)*A2;
    H = reshape(inv(A3)*A4,p,q);
    sigma2 = ss/N;
    err = abs(sigma2/sigma20-1);
    disp(['Iter: ' num2str(iter) ', S^2 = ' num2str(sigma2) ', Error = ' num2str(err)])
end

[U,S,V] = svd(H'*J0*H);
lambdas = diag(S);
etas = H*U*diag(1./sqrt(lambdas));