/*
******************************************************************
*   (C) Copyright 1999, Peter Lenk. All Rights Reserved.
*  GMIXREG.GSS
*		Generats data for HB Regression Model
*		Y_i 	= X_i*beta_i + epsilon_i
*		beta_i	= sum_{k=1}^K N(beta_i|theta_k,Lambda_k)
*		epsilon_i is N(0,sigma2*I)
*		Different Design Matrices
*****************************************************************
*/
new;
nsub	= 100;					@ Number of subjects 					@
yrows	= 10*ones(nsub,1);		@ Number of observations per subject @
@yrows	= 3 + ceil(7*rndu(nsub,1));@		@ Number of observations per subject		@
ntot	= sumc(yrows);
sigmat	= 5;					@ True error STD			@
rankx	= 2;					@ rank(X_i)		@
mmodt	= 3;					@ Three mixture componts @
@ thetat[.,j] is the mean for component j @
thetat  = {     0    -10     7,
                0      7     5  };
@ lambdat stacks 3 covariance matrics, each is a 2 by 2 matrix. @
lambdat = {  1    0,
             0    1,
             25   9,
             9    4,
             9   -5,
            -5    5   };



@ Get pointer into stacked xy matrices @
b	= cumsumc(yrows);
a	= 1|(1+b[1:nsub-1]);
iptxy	= a~b;

@ Get pointer into sacked lambdat @
b	= seqa(rankx,rankx,mmodt);
a	= 1|(1+b[1:mmodt-1]);
iptgp	= a~b;


xdim	= rankx - 1;
a	= seqa(1,1,xdim);
xnames	= 0 $+ "X " $+ ftocv(a,1,0);
xynames	= xnames|"Y ";

xdata	= rndn(ntot,xdim);
xmat	= ones(ntot,1)~xdata;
psit    = { 0.2,  0.3,  0.5 };		@ Mixture probabilities 	@
zprob   = ones(nsub,1).*.(psit');
z       = rndzmn(zprob);			@ Z gives class membership 	@

@ Compute Chol of lambdat @
lmbd12	= lambdat;
for i0 (1,mmodt,1); i = i0;
	lmbd12[iptgp[i,1]:iptgp[i,2],.] = chol(lambdat[iptgp[i,1]:iptgp[i,2],.]);
endfor;

@ Generate Beta and Ydata @
betat	= zeros(rankx,nsub);
ydata	= zeros(ntot,1);

for i0 (1,nsub,1); i = i0;
	betai	= thetat[.,z[i]] + lmbd12[iptgp[z[i],1]:iptgp[z[i],2],.]'rndn(rankx,1);
	betat[.,i]	= betai;
	xi	= xmat[iptxy[i,1]:iptxy[i,2],.];
	yi	= xi*betai + sigmat*rndn(yrows[i],1);
	ydata[iptxy[i,1]:iptxy[i,2],.] = yi;
endfor;
xydata	= xdata~ydata;
/*
**************************************************************************
* Create a Gauss file.  f1 is the file handle.  
* The Gauss file will be called "XYDATA."       
* The column will be named by the strings in the character array xyname.  
* ^xyname means use the names in the character string.	
* 0, 8 gives double precision real numbers.			
**************************************************************************
*/
create f1 = xydata with ^xynames, 0, 8;
/*
**************************************************************************
* Next read data into the Gauss file by using the writer command.	
* f1 is the file handle defined in previous command.			
* xydata is the data matrix that we just created.			
* writer returns the number of rows read to f1.			
* If it is not rows(xydata), something bad happended.		
**************************************************************************
*/
if writer(f1,xydata) /= rows(xydata);
		errorlog "Conversion of XYDATA to Gauss File did not work";
endif;
closeall f1;
save yrows		= yrows;
save classt		= z;		@ true classifications @
save sigmat		= sigmat;
save betat		= betat;
save thetat		= thetat;
save lambdat	= lambdat;
save psit		= psit;


@ Compute individual level MLEs @

bhat	= zeros(rankx,nsub);		
sse	= 0;
for i0 (1,nsub,1);	i = i0;
	xi	= xmat[iptxy[i,1]:iptxy[i,2],.];
	yi	= ydata[iptxy[i,1]:iptxy[i,2],.];
	xitxi	= xi'xi;
	xitxii	= invpd(xitxi);
	xityi	= xi'yi;
	bhat[.,i] = xitxii*xityi;
	resid	= yi - xi*bhat[.,i];
	sse	= sse + resid'resid;
endfor;
s2hat	= sse/ntot;		@ MLE of sigma2 @
sighat	= sqrt(s2hat);


@ Do some plots @
_plctrl = -1;
title("Y versus X");
xy(xydata[.,1],xydata[.,cols(xydata)]);
title("True Slope verus Intercept");
xy(betat[1,.]',betat[2,.]');
title("True & MLE Slope versus Intercept");
xy(betat[1,.]'~bhat[1,.]',betat[2,.]'~bhat[2,.]');
graphset;