% Script file name: Simulate_MLMM.m
% Purpose:
%   Simulate data base on MLMM.
% Define variables:
%   Input
%         : Number of groups (G)
%         : Number of subjects per group (Ng)
%         : Vector of Latent Class Probability
%         : Matrix of Conditional Latent Class Probability
%         : Matrix of Conditional Latent Transition Probability
%         : Matrix of Conditional Response Probability
%   Output
%         : `Data4D' (G \times Ng \times J \times T) data set
%         : save in MLMMM%dG%dNg%dJ%dT%dData%d with true values
DataN = input('How many data sets:');
G = input('Number of groups:');
Ng= input('Number of subjects per group:');
L= input('Input the number of latent clusters:');
M= input('Input the number of latent classes:');
T= input('Number of Time Points:');
%%  Specify parameters:
% $PH_g$
disp('Specify a vector a latent cluster probability:')
PH_g=input('Dimension: L \times 1 (e.g., [0.2; 0.8])\n');
True_PH_g=PH_g;

% PXgivenH $P(X_{ig}=m|H_{g}=l)$
disp('Specify a matrix of conditional latent class probability:')
disp('Dimension: (M \times L)\n Each column needs to sum up to 1')
PXgivenH=input('(e.g., [0.2 0.6;0.8 0.4])\n');
True_PXgivenH=PXgivenH;
[True_M True_L] = size(PXgivenH);

% PX_ig $P(X_{ig}=m)$
Px_ig=PXgivenH*PH_g;
True_PX_ig=Px_ig;

% Tau (Transition Matrix) For two time points
% TauGivenH (M M L)
disp('Specify the Transition Matrix by cluster')
disp('Dimension: (M \times M) \n Each row needs to sum up to 1')
disp('(e.g., [0.2 0.6;0.8 0.4])')
TauL1=input('Specify the transition matrix for the 1st cluster\n');
TauL2=input('Specify the transition matrix for the 2nd cluster\n');
TauGivenH=cat(3,TauL1,TauL2); % M M L
True_TauGivenH=TauGivenH; % M M L
TauGivenHT=reshape(repmat(TauGivenH,1,T),M,M,L,T);
% CondResProd
disp('Specify matrix of Conditional Response Prob:')
CondResProd=input('Matrix need to have the dimension of M X J:\n');
J=length(CondResProd(1,:));
True_CondResProd=CondResProd;

for DataID=1:DataN,
    %
    LatCluster=zeros(1,G);
    for g=1:G
        LatCluster(g) = sample_discrete(PH_g);
    end
    True_LatCluster=LatCluster;
    %
    LatClass=zeros(G,Ng,T);
    for g=1:G,
        for ng=1:Ng,
               LatClass(g,ng,1)=sample_discrete(PXgivenH(:,LatCluster(g)));
            for t=2:T,
          LatClass(g,ng,t)=...
            sample_discrete(TauGivenH(LatClass(g,ng,t-1),:,LatCluster(g)));
            end
        end
    end
    True_LatClass=LatClass;
    %
    for g=1:G,
        for ng=1:Ng,
            for j=1:J,
                for t=1:T,
                    Data4D(g,ng,j,t)=...
                        binornd(1,CondResProd(LatClass(g,ng,t),j));
                end
            end
        end
    end
    %
    fname=sprintf('MLMMM%dG%dNg%dJ%dT%dData%d',M,G,Ng,J,T,DataID);
    save(fname,'Data4D','True_PH_g','True_TauGivenH','True_CondResProd',...
        'True_PXgivenH','True_LatCluster','True_LatClass',...
        'True_M','True_L','J','T','G','Ng')
end