Gauss Code for Simulation Program
/* Generates random data for a two-class latent class model. Useful for conducting a parametric bootstrap. The latent class proportion for the first class is Thc; The conditional probabilities for the variables are Alpc[.,1] and Alpc[.,2] for
the first and second latent classes. Intermediate values for a given simulation prior to convergence are in Th, Thn, Alp and Alpn; after convergence of the EM computations the estimates are in ethn and ealpn. Note that computations are based on raw data in matrix, y. Simulator works well for reasonable numbers of items and sample size. */
t1=time;
nlc=2;
nconct=0;
/*                                       */
/* Enter the number of replications of the simulation */
/*                                       */
sim=500;
/*                                       */
/* Enter the number of variables, nk         */
/*                                       */
nk=4;
/*                                       */
/* Enter the sample size, nn          */
/*                                       */
nn=319;
ethn=zeros(sim,2);
ealpn=zeros(sim,2*nk);
p=zeros(nlc,nn);
/*                                       */
/* Enter parametric values, Thc and alpc         */
/*                                       */
Thc=.1606|.8394;
Alpc=.5769~.5890~.2160~.3764|.0166~.0292~.0371~.1820;
cut = cumsumc(Thc);
cat = cumsumc(ones(2,1));
/* Max iterations and convergence criterion for EM computations */
kter=300;
critc=.0001;
s=1;
do while s <= sim;
start:
/*  This statement actually generates random sample of nn cases and stores in y  */
y=1-round(rndu(nn,nk)-Alpc[subscat(rndu(nn,1),cut,cat),.]+.5);
print "step " s;
Th=Thc; Alp=Alpc;
iter=1;
pr = 1/nn;
/* Procedure to calculate probs for complete data */
Proc liki(Th,Alp);
local l,j,Tha,lik;
Tha=Th|(1-Th);
l = 1;
do while l <= nlc;
j=1;
do while j <= nn;
  p[l,j] = Tha[l]*prodc((Alp[l,.]^y[j,.])')*prodc(((1-Alp[l,.])^(1-y[j,.]))');
  j = j+1;
endo;
l = l + 1;
endo;
/* lik=(sumc(ln(sumc(p))));   */
retp(p);
endp;
/* Start of EM loop */
begin:
lik=liki(Th,Alp);
p =  (p'./sumc(p))';
Thn = sumc(p')*pr;
Alpn = Alp*0;
l = 1;
do until l == nlc+1;
j = 1;
do until j == nk+1;
       Alpn[l,j] = (p[l,.]*(pr.*y[.,j]))/Thn[l];
  j = j + 1;
endo;
l = l + 1;
endo;
crit = sumc(sumc(abs(Th-Thn)))+sumc(sumc(abs(Alp-Alpn)));
if crit < critc;
  goto end;
endif;
Th=Thn;
Alp=Alpn;
iter = iter+1;
  if iter == kter+1;
/*    print "convergence not attained"; */
    nconct=nconct+1;
    goto start;
  endif;
/*print crit;*/
goto begin;
end:
ethn[s,.]=Thn';
ealpn[s,1:nk]=Alpn[1,.];
ealpn[s,nk+1:2*nk]=Alpn[2,.];
s=s+1;
endo;
/*print ethn ealpn nconct;*/
print time-t1;
print "Did not converge" nconct "times";
print "Mean Theta" meanc(ethn)';
print "Mean Alpha" meanc(ealpn)';
print "Standard errors for Thetas" stdc(ethn)';
print "Standard errors for Alphas" stdc(ealpn)';
end