#twin.r - - estimation of p and d from twin data
 mm<- 32
 mf<- 41
 ff<- 36
 x<-c(mm,mf,ff)
# log-likelihood function
discrete.like<- function(x,p,d){	
    m<- 1-d
    q<- 1-p
    x[1]*log(p*m+p^2*d)+x[2]*log(2*p*q*d)+x[3]*log(q*m+q^2*d)
}

    dx<- 0.001
# initial values
    d<- 0.6
    p<- 0.4
    etemp<- NULL
    for (i in 1:10) {
        f0<- discrete.like(x,p,d)	
# approximate first partial derivatives	
	f1<- discrete.like(x,p+dx,d)
	f2<- discrete.like(x,p,d+dx)
	b1<-(f1-f0)/dx
	b2<- (f2-f0)/dx
# approximate second partial derivatives
	f11<- discrete.like(x,p+2*dx,d)
	f22<- discrete.like(x,p,d+2*dx)
	f12<- discrete.like(x,p+dx,d+dx)
	g11<- (f11-2*f1+f0)/dx^2
	g22<- (f22-2*f2+f0)/dx^2
	g12<- (f12-f1-f2+f0)/dx^2
	est<- rbind(p,d)
	b<- cbind(b1,b2)
	g<- cbind(c(g11,g12),c(g12,g22))
	etemp<- est-solve(g)%*%t(b)
	p<- etemp[1]
	d<- etemp[2]
	if(sum(abs(b))<0.000001) break
    }
cat("p=",round(p,3),"d=",round(d,3),"\n")
gtemp<- -solve(g)
"variance/covariance"
round(gtemp,6)