
%toolbox fdaM of Matlab should be installed and put in the path (in example> mathtest are the functions used)
%part of the code from mathtest.m

%  ----------------------------------------------------------------
%  ----------------------  enter ACT male data  -------------------
%  ----------------------------------------------------------------

nit =   60;
nex = 2115;

%  Read in the data.  The first record contains the key,
%    and subsequent records the responses for each examinee

fid = fopen('actm.txt','rt');
ACTmtest = reshape(fscanf(fid, '%s'), [nit,nex])';

%  set up limits on quadrature points

thetamax = floor(20*(1 + log10(nex))/3)/10;
thetamin = -thetamax;
thetarng = [thetamin, thetamax];

%  Smooth the data using the TestGraf algorithm that
%  estimates item response functions by kernel smoothing
%  The objective is to get an initial estimate of the
%  item response functions

nbin = 51;  %  number of bin boundaries
thetavec = linspace(thetamin, thetamax, nbin)';  %  bin boundaries

%  bin the data

Pbin0  = FirstStep(ACTmtest, thetavec);

%  set up smoothing matrix used in TestGraf

Smat = smoothmat(thetavec, nex);

%  smooth the binned data

Pbin = Smat*Pbin0;

%  compute log odds of smoothed data

Wbin  = log(Pbin./(1-Pbin)); 

%   set up B-spline basis  

nbasis = 11;
norder = 6;
wbasis = create_bspline_basis(thetarng, nbasis, norder);

%  set up a functional data object for initial W

Wfd0  = data2fd(Wbin, thetavec, wbasis); 

%  compute asymptotic standard error

Wvec  = eval_fd(Wfd0,0);
DWvec = eval_fd(Wfd0,1);
Pvec  = exp(Wvec)./(1+exp(Wvec));
serr  = 1/sqrt(sum(DWvec.^2.*Pvec.*(1-Pvec)));

%  compute number of quadrature points

nq = floor((thetamax - thetamin)/serr);  % nq = 26

%  set up equally-spaced trait values theta_q and thetafine

thetaq    = linspace(thetamin, thetamax, nq)';  %  quadrature points
thetafine = linspace(thetamin, thetamax, 101)'; %  plotting   points

%  least squares weights for B-spline test functions

norder = 4;
nbasis = nq + norder - 2;
wgtq   = gausswgt(thetaq, nbasis, norder);

% Matrix phimat has basis function values for each thetaq value

phimat = getbasismatrix(thetaq, wbasis);

%  set up linear differential operator

Ldata = zeros(101,3);
theta = linspace(thetamin,thetamax,101)';
Ldata(:,3) = (exp(theta)-1)./(exp(theta)+1);
wfd = data2fd(Ldata,theta,wbasis);
Lfdobj = Lfd(3,fd2cell(wfd));

% Matrix K is the penalty matrix for penalizing roughness;

Kmat = eval_penalty(wbasis, Lfdobj);  %Penalize size of 2nd derivative

%  ------------------------------------------------------------------
%                         The EM Algorithm
%  ------------------------------------------------------------------

%  Function FirstStep computes the proportions of examinees
%    associated with each value of thetaq that pass items

Pbin  = FirstStep(ACTmtest, thetaq);
Wbin  = log(Pbin./(1-Pbin)); % matrix of logit transforms of probabilities
Wfd0  = data2fd(Wbin, thetaq, wbasis); % functions corresponding to these
coef0 = getcoef(Wfd0);         %  coefficients for the basis expansion
W0    = eval_fd(Wfd0, thetaq); % evaluate functions at thetaq values
P0    = 1./(1+exp(-W0));    %  compute corresponding probabilities
P     = P0; Q = 1-P;        % starting values for probabilities
coef  = coef0;              % starting values for coefficients
Wfd   = Wfd0;               % starting functions w

%  ------------  initialize the EM algorithm  --------

lambda   = 1e0;              % smoothing parameter
penmat   = lambda.*Kmat;      % penalty matrix times smoothing parameter
iter     = 0;                 % initialize iteration number
convtest = 1e-3;              % convergence criterion
itermax  = 100;                % maximum number of iterations
F        = 1e10;              % initialize new and old function values
Fold     = F + 2*convtest;

%  ---- iterate through EM algorithm  calling functions Estep and Mstep --------

while Fold - F > convtest & iter < itermax
   if iter == 0, disp('It.    -log L   penalty    F     Fold - F'); end
   iter = iter + 1;
   Fold = F;

  %  --------  E step  ------------

  [N, CN, CP, L, CL] = Estep(ACTmtest, P, wgtq);
   
  %  compute penalized negative sum of marginal likelihoods
  logL = sum(log(L));
  pen  = sum(diag(coef' * penmat * coef));
  F = -logL + pen;
  %  print out F, which is being minimized
  fprintf('%g  ', [iter, -logL, pen, F, Fold - F]);  fprintf('\n');
   
  %   -------  Mstep  -------------
  
  coef = Mstep(CN, N, P, coef, phimat, penmat);
  % update Q by n matrix of probabilities 
  P = 1./(1+exp(-phimat * coef));  
  Q = 1 - P;
  
end

Wfd = putcoef(Wfd, coef);

%  --------------------  end of EM algorithm ---------------------

%  set results for fine mesh of theta values

Wfdfine = eval_fd(thetafine, Wfd);
Pfine   = exp(Wfdfine)./(1 + exp(Wfdfine));
Qfine   = 1 - Pfine;

%Figure
subplot(1,2,1)

plot(thetafine,Pfine)
axis([thetamin, thetamax, 0, 1])
xlabel('\fontsize{16} \theta')
ylabel('\fontsize{16} P(\theta)')



subplot(1,2,2)
%plot(Wfd)
%plot(Wfd,0,1,0)
plot(thetafine,Wfdfine)
axis([thetamin, thetamax, -4.6, 7.3])
xlabel('\fontsize{16} \theta')
ylabel('\fontsize{16} W(\theta)')




%%%PCA


nharm = 4;
Wpcastr = pca_fd(Wfd, nharm);
Wpcastr = varmx_pca(Wpcastr);

Wharmmat = eval_fd(Wpcastr.harmfd, thetaq);

Wfdmean = mean(Wfd);
Wmean = eval_fd(Wfdmean, thetaq);
Wmeanmat = Wmean*ones(1,4);

Pmean = exp(Wmean)./(1 + exp(Wmean));


Wconst = ones(nq,1)*[1, 1.5, .5, .5]; 
Wmatp = Wmeanmat + Wconst.*Wharmmat;
Wmatm = Wmeanmat - Wconst.*Wharmmat;
Pharp = exp(Wmatp)./(1+exp(Wmatp));
Pharm = exp(Wmatm)./(1+exp(Wmatm));

%Wpcastr

%varprop: [0.4281 0.1965 0.1006 0.2743] 0.9995

%Figure
titlestr = {' I: 43%'; ' II: 20%'; 'III: 10%'; ' IV: 27%'}; 
for j=1:nharm
subplot(2,2,j)
plot(thetaq, Pmean, '--')
axis([thetamin, thetamax, 0, 1])
text(thetaq-.1, Pharp(:,j), '+')
text(thetaq-.1, Pharm(:,j), '-')
title(['\fontsize{12} Harmonic ',titlestr{j}])
end


%Figure
subplot(1,2,1)

plot(Wpcastr.harmscr(:,1),Wpcastr.harmscr(:,2),'.')
 text(Wpcastr.harmscr(:,1)+.1,Wpcastr.harmscr(:,2)+.1,cellstr(num2str([1:60]')))
xlabel('\fontsize{16} Component 1 Score')
ylabel('\fontsize{16} Component 2 Score')


subplot(1,2,2)

plot(Wpcastr.harmscr(:,3),Wpcastr.harmscr(:,4),'.')
 text(Wpcastr.harmscr(:,3)+.1,Wpcastr.harmscr(:,4)+.1,cellstr(num2str([1:60]')))
xlabel('\fontsize{16} Component 3 Score')
ylabel('\fontsize{16} Component 4 Score')


%%%%%
%For FADA in R, 4 de Wfd

%save('Wfd.mat','Wfd')

%Figure 
%FADA

plot(thetafine,Pfine(:,[2]),'-g')
hold on
plot(thetafine,Pfine(:,[18]),'b:')
hold on
plot(thetafine,Pfine(:,[28]),'k-.')
hold on
plot(thetafine,Pfine(:,[60]),'r--')

%hold on 
%plot(thetaq,Pmean,'y')
axis([thetamin, thetamax, 0, 1])
xlabel('\fontsize{16} \theta')
ylabel('\fontsize{16} P(\theta)')

hleg1 = legend('2','18','28','60');

set(hleg1,'Location','SouthEast')
