/* --- 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,b3,g1,g2) := 
     normpdf( y, b1*(2*g1-1) + b2*(2*g2-1) + b3*(2*g1-1)*(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,b3,q1,q2) := ypdf(y,b1,b2,b3,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,b3,q1,q2) := sum( sum( ygpdf(y,g1,g2,b1,b2,b3,q1,q2),
    g1,0,1 ), g2,0,1);

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

gpostpdf(g1,g2,y,b1,b2,b3,q1,q2) := ygpdf(y,g1,g2,b1,b2,b3,q1,q2)/
                                 ymarpdf(y,b1,b2,b3,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 : sum(sum(diff(diff(
            -log(gpostpdf(g1,g2,y,b1,b2,b3,q1,q2)),b1),b1)
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1);
im12 : sum(sum(diff(diff(
            -log(gpostpdf(g1,g2,y,b1,b2,b3,q1,q2)),b1),b2)
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1);
im13 : sum(sum(diff(diff(
            -log(gpostpdf(g1,g2,y,b1,b2,b3,q1,q2)),b1),b3)
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1);
im21 : sum(sum(diff(diff(
            -log(gpostpdf(g1,g2,y,b1,b2,b3,q1,q2)),b2),b1)
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1);
im22 : sum(sum(diff(diff(
            -log(gpostpdf(g1,g2,y,b1,b2,b3,q1,q2)),b2),b2)
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1);
im23 : sum(sum(diff(diff(
            -log(gpostpdf(g1,g2,y,b1,b2,b3,q1,q2)),b2),b3)
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1);
im31 : sum(sum(diff(diff(
            -log(gpostpdf(g1,g2,y,b1,b2,b3,q1,q2)),b3),b1)
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1);
im32 : sum(sum(diff(diff(
            -log(gpostpdf(g1,g2,y,b1,b2,b3,q1,q2)),b3),b2)
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1);
im33 : sum(sum(diff(diff(
            -log(gpostpdf(g1,g2,y,b1,b2,b3,q1,q2)),b3),b3)
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1);

/* make the missing information matrix */

Im : matrix( [im11,im12,im13], [im21,im22,im23], [im31,im32,im33] );

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(diff(
            -log(ygpdf(y,g1,g2,b1,b2,b3,q1,q2)),b1),b1)
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1));
ic12 : ratsimp(sum(sum(ratsimp(diff(diff(
            -log(ygpdf(y,g1,g2,b1,b2,b3,q1,q2)),b1),b2))
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1));
ic13 : ratsimp(sum(sum(ratsimp(diff(diff(
            -log(ygpdf(y,g1,g2,b1,b2,b3,q1,q2)),b1),b3))
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1));
ic21 : ratsimp(sum(sum(ratsimp(diff(diff(
            -log(ygpdf(y,g1,g2,b1,b2,b3,q1,q2)),b2),b1))
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1));
ic22 : ratsimp(sum(sum(ratsimp(diff(diff(
            -log(ygpdf(y,g1,g2,b1,b2,b3,q1,q2)),b2),b2))
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1));
ic23 : ratsimp(sum(sum(ratsimp(diff(diff(
            -log(ygpdf(y,g1,g2,b1,b2,b3,q1,q2)),b2),b3))
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1));
ic31 : ratsimp(sum(sum(ratsimp(diff(diff(
            -log(ygpdf(y,g1,g2,b1,b2,b3,q1,q2)),b3),b1))
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1));
ic32 : ratsimp(sum(sum(ratsimp(diff(diff(
            -log(ygpdf(y,g1,g2,b1,b2,b3,q1,q2)),b3),b2))
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1));
ic33 : ratsimp(sum(sum(ratsimp(diff(diff(
            -log(ygpdf(y,g1,g2,b1,b2,b3,q1,q2)),b3),b3))
            *gpostpdf(g1,g2,y,b1,b2,b3,q1,q2),g1,0,1),g2,0,1));

Ic : matrix( [ic11,ic12,ic13], [ic21,ic22,ic23], [ic31,ic32,ic33] );

/* the observed information matrix */

Io : Ic-Im;
Io : ratsimp(Io);


/* Now we simplify the observed information matrix for the case when the
   epistatic effect is small (b3=0), one QTL has small effect (b2=0), as
   a function of the effect of one qtl (b1=b). */

A : subst(0,b3,Io);
A : subst(0,b2,A);
A : subst(b,b1,A);
A : ratsimp(A);

ratsimp(expand(num(A[1,1])));
ratsimp(expand(denom(A[1,1])));
ratsimp(expand(num(A[2,3])));
ratsimp(expand(num(A[3,3])));
ratsimp(expand(denom(A[3,3])));


/* Now we simplify the observed information matrix for the case when
   the epistatic effect is small (b3=0), as a function of the additive
   QTL of equal effect (b1=b, b2=b). */

B : subst(0,b3,Io);
B : subst(b,b2,B);
B : subst(b,b1,B);
B : ratsimp(B);


/* Both these expressions look quite complicated.  We will now try to 
   extract the important elements so what we can write down a formula
   for numerical computation. */

/* Although the expressions look daunting, they have a common denominator.
   So, to simplyfy matters, we multiply the matrix by that denominator. */

C : denom(B[1,1])*B;

/* The denominator is a polynomial of exp(2*b*y) and exp(2*b^2). */

D : expand(denom(B[1,1]));

/* we get each of the parts as follows */

a1 : part(D,1,2);
a2 : part(D,2,2);
a3 : part(D,3,2);
a4 : part(D,4);
a5 : part(D,5,2);

/* this is a function to get the coefficients of an expression (expr)
   corresponding to a1, a2, a3, a4, a5 after simplification. */

getcoef(expr,a1,a2,a3,a4,a5) :=	block(
c1 : coeff(expr,a1),
c2 : coeff(expr,a2),
c3 : coeff(expr,a3),
c4 : coeff(expr,a4),
c5 : coeff(expr,a5),
cc : [c1,c2,c3,c4,c5],
aa : [a1,a2,a3,a4,a5],
c6 : ratsimp(expr - aa.cc),
cc : [c1,c2,c3,c4,c5,c6],
return(cc));

aa : [a1,a2,a3,a4,a5,1];

/* for each of the non-trivial elements of the information matrix, 
   we extract the coefficients corresponding to the polynomial in 
   the denominator */

z11 : getcoef(expand(C[1,1]),a1,a2,a3,a4,a5);
z12 : getcoef(expand(C[1,2]),a1,a2,a3,a4,a5);
z13 : getcoef(expand(C[1,3]),a1,a2,a3,a4,a5);
z33 : getcoef(expand(C[3,3]),a1,a2,a3,a4,a5);

/* print them in transpose as it is easier to see them */

transpose(z11);
transpose(z12);
transpose(z13);
transpose(z33);