mu<- 4 sd<- 2 n<- 50 x<- rnorm(n,mu,sd) # maximum likelihood estimated mean round(mean(x),3) # maximum likelihood estimated standard deviation round(sqrt(((n-1)/n)*var(x)),3) norm.like<- function(n,x,m,s){ -n*log(s)-0.5*sum( ( (x-m)/s)^2) } dx<- 0.001 # initial values m<- 3.0 s<- 1.0 n<- 50 #x<- rnorm(n,mu,sd) etemp<- NULL cat("Iteration","Mean","Standard Deviation","Gradient","H11","H12","H22","\n") for (i in 1:50) { f0<- norm.like(n,x,m,s) # approximate first partial derivatives f1<- norm.like(n,x, m+dx,s) f2<- norm.like(n,x, m,s+dx) b1<-(f1-f0)/dx b2<- (f2-f0)/dx # approximate second partial derivatives f11<- norm.like(n,x, m+2*dx,s) f22<- norm.like(n,x, m,s+2*dx) f12<- norm.like(n,x, m+dx,s+dx) g11<- (f11-2*f1+f0)/dx^2 g22<- (f22-2*f2+f0)/dx^2 g12<- (f12-f1-f2+f0)/dx^2 est<- rbind( m,s) b<- cbind(b1,b2) g<- cbind(c(g11,g12),c(g12,g22)) etemp<- est-solve(g)%*%t(b) m<- etemp[1] s<- etemp[2] cat(" ",i," ",round( m,3)," ",round(s,3)," ",round(b1,3)," ",round(b2,3), " ",round(g11,3)," ",round(g12,3)," ",round(g22,3),"\n") if(sum(abs(b))<0.000001) break } loglik<-norm.like(n,x,m,s) cat("m=",round( m,3),"s=",round(s,3),"loglik=",round(loglik,3),"\n") gtemp<- -solve(g) "variance/covariance" round(gtemp,6)