TITLE1 'San Diego Homeless Data - Estimated Marginal Probabilities';
PROC IML;

/* Results from MIXNO analysis - output file SDHOUSE.OUT  */;

w0 = { 0 0 0 0 0 0 0,
       0 1 0 0 0 0 0,
       0 0 1 0 0 0 0,
       0 0 0 1 0 0 0};
w1 = { 1 0 0 0 0 0 0,
       1 1 0 0 1 0 0,
       1 0 1 0 0 1 0,
       1 0 0 1 0 0 1};

mua   = {-0.45189};
sda   = { 0.87071};

alphaa= { 0.52099,
          1.94152,
          2.81975,
          2.25949,
         -0.13527,
         -1.91710,
         -0.95229};

mub   = {-2.67540};
sdb   = { 2.33380};


alphab= { 0.78108,
          2.68209,
          4.08805,
          4.09897,
          2.00273,
          0.54838,
          0.30403};


/* Now get the estimated marginal probabilities */;
/* number of quadrature points, quadrature nodes & weights  */
nq = 10;
bq = {  -4.85946282833231,  -3.58182348355193,  -2.48432584163895,
        -1.46598909439116,  -0.48493570751550,   0.48493570751550,
         1.46598909439116,   2.48432584163895,   3.58182348355193,
         4.85946282833231};
aq = {   0.00000431065265,   0.00075807095698,   0.01911158107317,
         0.13548370704150,   0.34464234526294,   0.34464234526294,   
         0.13548370704150,   0.01911158107317,   0.00075807095698,
         0.00000431065265};


/* initialize to zero */
grp0  = J(4,1,0);
grp0a = J(4,1,0);
grp0b = J(4,1,0);
grp1  = J(4,1,0);
grp1a = J(4,1,0);
grp1b = J(4,1,0);

DO q = 1 to nq;

za0 = w0*alphaa + mua + sda*bq[q];
zb0 = w0*alphab + mub + sdb*bq[q];
za1 = w1*alphaa + mua + sda*bq[q];
zb1 = w1*alphab + mub + sdb*bq[q];

grp0 =grp0 + (1/(1+(exp(za0)+exp(zb0))))*aq[q];
grp0a=grp0a+ (exp(za0)/(1+(exp(za0)+exp(zb0))))*aq[q];
grp0b=grp0b+ (exp(zb0)/(1+(exp(za0)+exp(zb0))))*aq[q];

grp1 =grp1 + (1/(1+(exp(za1)+exp(zb1))))*aq[q];
grp1a=grp1a+ (exp(za1)/(1+(exp(za1)+exp(zb1))))*aq[q];
grp1b=grp1b+ (exp(zb1)/(1+(exp(za1)+exp(zb1))))*aq[q];

END;

print 'Quadrature method - 10 points';

print 'marginal probability for group 0 - category 0' grp0  [FORMAT=8.4];
print 'marginal probability for group 0 - category 1' grp0a [FORMAT=8.4];
print 'marginal probability for group 0 - category 2' grp0b [FORMAT=8.4];

print 'marginal probability for group 1 - category 0' grp1  [FORMAT=8.4];
print 'marginal probability for group 1 - category 1' grp1a [FORMAT=8.4];
print 'marginal probability for group 1 - category 2' grp1b [FORMAT=8.4];