/* --- information for a qtl with small effect in the presence of a second 
unlinked qtl of varying effect --- */

/* define normal density function */

normpdf(x,m,s) := exp( -(1/2)*( (x-m)/s )^2 ) / (s*sqrt(2*%PI));

/* define binomial density function */

binomialpdf(x,p) := if (x=0) then (1-p) else if (x=1) then p else 0;

/* conditional distribution of phenotype given genotype */
/* note the 2-qtl model */

ypdf(y,b1,b2,g1,g2) := normpdf(y,b1*(2*g1-1)+b2*(2*g2-1),1);

/* prior distribution of genotype */
/* the two qtls are unlinked and hence independent */

gpdf(g1,g2,q1,q2) := binomialpdf(g1,q1)*binomialpdf(g2,q2);

/* joint distribution of phenotype and genotype */
/* we get this by just multiplying the prior with the likelihood */

ygpdf(y,g1,g2,b1,b2,q1,q2) := ypdf(y,b1,b2,g1,g2) * gpdf(g1,g2,q1,q2);

/* marginal distribution of phenotype */
/* obtained by integrating (summing) over the missing data (g1 and g2) */

ymarpdf(y,b1,b2,q1,q2) := sum( sum( ygpdf(y,g1,g2,b1,b2,q1,q2),
    g1,0,1 ), g2,0,1);

/* posterior distribution of genotype */
/* joint dist / marginal dist */

gpostpdf(g1,g2,y,b1,b2,q1,q2) := ygpdf(y,g1,g2,b1,b2,q1,q2)/
                                 ymarpdf(y,b1,b2,q1,q2);

/* establish constraints */
/* q's lie between 0 and 1 */

assume(q1>0);
assume(q1<1);
assume(q2>0);
assume(q2<1);

/* for simplicity assume that we are at an ungenotyped location */
q1:1/2;
q2:1/2;

/* information of the missing data likelihood */

/* we get each entry in the missing information matrix by
differentiating the missing data likelihood twice and then summing
over the posterior distribution of the qtl genotypes */

im11 : ratsimp(sum(sum(diff(diff(
            -log(gpostpdf(g1,g2,y,b1,b2,q1,q2)),b1),b1)
            *gpostpdf(g1,g2,y,b1,b2,q1,q2),g1,0,1),g2,0,1));
im12 : ratsimp(sum(sum(ratsimp(diff(diff(
            -log(gpostpdf(g1,g2,y,b1,b2,q1,q2)),b1),b2))
            *gpostpdf(g1,g2,y,b1,b2,q1,q2),g1,0,1),g2,0,1));
im21 : ratsimp(sum(sum(ratsimp(diff(diff(
            -log(gpostpdf(g1,g2,y,b1,b2,q1,q2)),b2),b1))
            *gpostpdf(g1,g2,y,b1,b2,q1,q2),g1,0,1),g2,0,1));
im22 : ratsimp(sum(sum(ratsimp(diff(diff(
            -log(gpostpdf(g1,g2,y,b1,b2,q1,q2)),b2),b2))
            *gpostpdf(g1,g2,y,b1,b2,q1,q2),g1,0,1),g2,0,1));

/* make the missing information matrix */

Im : matrix( [im11,im12], [im21,im22] );

/* simplify the expression by considering the case when the first qtl
has negligibly small effect; subsitute 0 for b1 */

Im : subst(0,b1,Im);
Im : ratsimp(Im);

/* complete information matrix */

/* obtained by taking the second derivative of the complete data
likeihood and then summing over the posterior distribution of the
missing data (qtl genotypes) */

ic11 : ratsimp(sum(sum(diff(
            -log(ygpdf(y,g1,g2,b1,b2,q1,q2)),b1,2)
            *gpostpdf(g1,g2,y,b1,b2,q1,q2),g1,0,1),g2,0,1));
ic12 : ratsimp(sum(sum(ratsimp(diff(diff(
            -log(ygpdf(y,g1,g2,b1,b2,q1,q2)),b1),b2))
            *gpostpdf(g1,g2,y,b1,b2,q1,q2),g1,0,1),g2,0,1));
ic21 : ratsimp(sum(sum(ratsimp(diff(diff(
            -log(ygpdf(y,g1,g2,b1,b2,q1,q2)),b2),b1))
            *gpostpdf(g1,g2,y,b1,b2,q1,q2),g1,0,1),g2,0,1));
ic22 : ratsimp(sum(sum(ratsimp(diff(diff(
            -log(ygpdf(y,g1,g2,b1,b2,q1,q2)),b2),b2))
            *gpostpdf(g1,g2,y,b1,b2,q1,q2),g1,0,1),g2,0,1));

a : gpostpdf(0,0,y,b1,b2,q1,q2) - gpostpdf(0,1,y,b1,b2,q1,q2) -
	gpostpdf(1,0,y,b1,b2,q1,q2) + gpostpdf(1,1,y,b1,b2,q1,q2);

Ic : matrix( [ic11,ic12], [ic21,ic22] );
Ic : subst(0,b1,Im);