options ls=80;
libname sample '.';
proc iml;
    create lrdat var {loc l1 cm a sig2e l0 la lrt };
    use sample.qtlcart;
    read ALL;
    use sample.qtlcart;
	read all var ('c1m1':'c1m16') into c1;
	map={0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150};
    nobs=nrow(t1);
    y=t1;
    chrom=1;
    base=0;
    do marker=1 to 15;
    M=c1[,marker];
    N=c1[,marker+1]; 
    i=0;
	*print marker;
    do j=1 to nobs;
        if y[j,1] ^= . & M[j,1] ^= . & N[j,1] ^= .  then do;
            i=i+1;
            y[i,1]=y[j,1];
            M[i,1]=M[j,1];
            N[i,1]=N[j,1];
            end;
    end;
    nobs=i;
    y=y[1:nobs,];
    M=M[1:nobs,];
    N=N[1:nobs,];
	*print m n y;
    morgan=(map[1,marker+1]-map[1,marker])/100;
    R=.5*(1-exp(-2*morgan));
    morgan=-.5*log(1-2*r);

    X=j(nobs,1)||design(c1m9)||design(c2m4)||design(c2m14)||design(c3m9)||design(c3m14);
	if marker>6 & marker<11 then X=j(nobs,1)||design(c2m4)||design(c2m14)||design(c3m9)||design(c3m14);
	*X=j(nobs,1,1);
    bhat=ginv(x`*x)*x`*y;
    ehat=y-x*bhat;
    sse=ehat`*ehat;
    sig2e=sse/nobs;
    L0=-(nobs/2)*log(sig2e)+sum(log(exp(-.5*ehat`*ehat/sig2e)));
	l00=-(nobs/2)*log(sig2e)-nobs/2;
	*print l0 l00 bhat sig2e;
    bhat=bhat//{0};
    X1=(j(nobs,1,0))||X;
    X2=(j(nobs,1,-1))||X;
    lrt=0;

    do l1= 0.0001 to morgan by 2/100;  /* Walk the interval in 2 cM units */
        loc=base+l1;
        r1=.5*(1-exp(-2*l1));
        l2=morgan-l1;
        r2=.5*(1-exp(-2*l2));
        p2=  (M=4)# (N=4)*r1*r2/(1-R)
            +(M=2)# (N=2)*(1-r1*r2/(1-R))
            +(M=4)# (N=2)*r1*(1-r2)/R
            +(M=2)# (N=4)*r2*(1-r1)/R;
        p1=1-p2;
        w1=p1;
        w2=p2;
        
        do until(-1.e-20 < delt & delt < 1.e-20);
            lrtold=lrt;
            alpha1=w1/(w1+w2);
            alpha2=w2/(w1+w2);
            xpx=x1`*diag(alpha1)*x1+x2`*Diag(alpha2)*x2;
            xpy=x1`*diag(alpha1)*y+x2`*diag(alpha2)*y;
            bhat=ginv(xpx)*xpy;
            e1=y-x1*bhat;e2=y-x2*bhat;
            sse=e1`*diag(alpha1)*e1+e2`*diag(alpha2)*e2;
            sig2e=sse/nobs;
            w1=p1#exp(-.5*e1#e1/sig2e);
            w2=p2#exp(-.5*e2#e2/sig2e);
            La=-(nobs/2)*log(sig2e)+sum(log(w1+w2));
            lrt=2*(la-l0);
            delt=lrt-lrtold;
            end;
        a=bhat[1,1];
        cm=l1*100+map[1,marker];
        print cm sig2e a l0 la lrt delt;
        append;
        end;

	base=base+morgan;
    end;
 quit;
run;
proc print data=lrdat;
var cm sig2e a l0 la lrt;
run;
symbol1 i=join line=1 value=none color=black ;
symbol2 i=join line=3 value=none color=black;
axis1 label=(angle=90 'LRT' ) order=(0 to 40 by 10);
axis2 label=('Location (cM)');
legend1 label=("Design") Value=("Backcross" "F2");
goption device=psepsf  hsize=7in vsize=4in gsfmode=replace;

filename grafout 'lrt.eps';
goption gsfname=grafout;
proc gplot data=lrdat;
    plot lrt*cm/ haxis=axis2 vaxis=axis1;
run;