clear all;

% This is an example of showing the codes of "FSIM"
% it takes about 15 mins when n=100

n = 100; % sample size (number of subjects)

[y z t xi] = sim_FSIM(n); % generate data


% set parameters
bwmu = [0  0];         % bandwidths for mean (set [0 0] to apply CV)
bwxcov = 0;            % bandwidths for covariance (set [0 0] to apply CV)
mfc = 10;              % number of fold for CV (while estimating covariance)
FVE_threshold=0.80;    % threshold of fraction of variation explained 
                       % (to select the number of eigenfunctions)

error=1;             % measurement error with constant variance exists
                     %  0: no measurement errors
                     %  2: the variance of measurement errors varies with
                     %  time
                     
ngrid=51;            % number of grid for output

selection_k = 'FVE'; % methods of selecting the number of eigenfunctions
                     % other choices: AIC, BIC

tic     
[noeig, Esig, lambda, eigen, xi_est, beta, muout, bw_mu, ecov, bw_xcov, ...
           yfit, tout]...
        = FSIM(y, z, t, bwmu, bwxcov, mfc, selection_k, FVE_threshold, ...
          error, ngrid);
toc      
        
nt = length(tout);

bt=[2 1 0 3 0 -1]'./sqrt(15); % true beta


z2 = cell2mat(z)*beta;    % estimate index

z1 = min(z2):(range(z2)/(ngrid-1)):max(z2);

Tmu = zeros(ngrid,nt); % True Mean function

    for j=1:ngrid
        Tmu(j,:) = sin(z1(j))*sin(tout*pi)+(1-sin(z1(j)))*cos(tout*pi); 
    end    
 
% plots of true and estimated mean functions    
figure;
subplot(1,2,1)
surf(tout,z1,Tmu);
zlim([-3 3]);
xlabel('z')
ylabel('time')
title('True Mean Surface Function');
subplot(1,2,2)
surf(tout,z1,muout);
zlim([-3 3]);
xlabel('z')
ylabel('time')
title('Estimated Mean Surface Function');

fprintf(1,'True beta:\n');
bt
fprintf(1,'Estimated beta:\n');
beta

