options ls=80;
libname barley 'barley_data';
proc iml;
    create lrdat var {l1 cm a sig2e l0 la lrt };
    use barley.barley1;
    read ALL;
    nobs=nrow(env9);
    y=env9;
    M=MWG036B;
    N=RISIC10A;
    i=0;
	print m n y;
    do j=1 to nobs;
        if y[j,1] > 0 & 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=2.7439/100;
    R=.5*(1-exp(-2*morgan));
    morgan=-.5*log(1-2*r);
    X=j(nobs,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=X||(j(nobs,1,0));
    x2=X||(j(nobs,1,-1));
    lrt=0;
    do l1= 0.0001 to morgan by .1/100;  /* Walk the interval in 1 cM units */
        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-10 < delt & delt < 1.e-10);
            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[2,1];
        cm=l1*100;
        print l1 bhat sig2e l0 la lrt delt nobs;
        append;
        end;
    quit;
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 5 by .5);
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;