% Script file name: Simulate_MLCM.m
% Purpose:
%   Simulate data base on MLCM.
% 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 Response Probability
%   Output
%         : `DataGNgJ' (G \times Ng \times J) data set
%         : save in MLCMM%dG%dNg%dJ%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:');
%%  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 P(M|L): Dimension: (M \times L) Each column needs to sum 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;

% CondResProd
disp('Specify matrix of Conditional Response Prob:')
CondResProd=input('\n 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);
    for g=1:G,
        for ng=1:Ng,
            LatClass(g,ng)=sample_discrete(PXgivenH(:,LatCluster(g)));
        end
    end
    True_TauGivenH=PXgivenH;
    True_LatClass=LatClass;
    %
    for g=1:G
        for ng=1:Ng
            for j=1:J
                DataGNgJ(g,ng,j)=binornd(1,CondResProd(LatClass(g,ng),j));
            end
        end
    end
    %
    fname=sprintf('MLCMM%dG%dNg%dJ%dData%d',M,G,Ng,J,DataID);
    save(fname,'DataGNgJ','True_PH_g','True_TauGivenH','True_CondResProd',...
        'True_LatCluster','True_LatClass','True_M','True_L','J','G','Ng')
end