% Script file name: Simulate_LMM.m
% Purpose:
%   Simulate data base on LMM.
% Define variables:
%   Input
%         : Number of subjects (N)
%         : Vector of Latent Class Probability
%         : Matrix of Conditional Response Probability
%         : Matrix of Transition Probability (Tau)
%         : Number of subjects (N)
%   Output
%         : `DataNJT' (N \times J \times T) data set
%         : save in LMMM%dN%dJ%dT%dData%d with true values
DataN = input('How many data sets:');
N= input('Number of subjects:');
M= input('Input the number of latent classes:');
PX=input('Specify latent class probability:\n(e.g., [0.3 0.7])\n');
disp('Specify matrix of Conditional Response Prob:')
CondResProd=input('(Matrix need to have the dimension of M X J)\n');
True_M=length(PX);
True_PX=PX;
Tau=input('Matrix of Transition Probabilities (M X M):\n');
True_Tau=Tau;
J=length(CondResProd(1,:));
T= input('Number of time points:');
True_T=T;
True_CondResProd=CondResProd;
for DataID=1:DataN,
    S = zeros(N,T);
    for i=1:N
        S(i, 1) = sample_discrete(PX);
        for t=2:T
            S(i, t) = sample_discrete(Tau(S(i,t-1),:));
        end
    end
    TrueState=S; % N by T
    DataNJT=zeros(N,J,T);
    for i=1:N
        for t=1:T
            for j=1:J
                DataNJT(i,j,t)=binornd(1,CondResProd(S(i,t),j));
            end
        end
    end
    Start_PX=True_PX;
    Start_Tau=reshape(repmat(Tau,1,T),M,M,T);
    Start_CondResProd=True_CondResProd;
    fname=sprintf('LMMM%dN%dJ%dT%dData%d',M,N,J,T,DataID);
    save(fname,'DataNJT','True_PX','True_CondResProd','TrueState',...
        'Start_PX','Start_Tau','Start_CondResProd','True_M','J','N');
end