R code:
### The following R code does just this for a normal linear regression (eq1 above): # Generate some data a <- 5 b <- 2 sigma <- 0.5 x <- runif(1000,0,20) y <- a+(b*x)+rnorm(1000)*sigma dat <- data.frame(x=x,y=y) plot(dat$x,dat$y) # Write the likelihood function for the statistical model (in this case, normal linear regression) loglike<- function(par,dat) { a <- par['a'] b <- par['b'] sigma <- par['sigma'] return(sum(dnorm(dat$y-a-(b*dat$x),mean=0,sd=sigma,log=T))) } # Maximise the log likelihood res <- optim(par=c(a=1.5,b=1.5,sigma=0.6), fn=loglike, dat=dat, method='L-BFGS-B', lower=rep(0.2,3), upper=c(a=5,b=5,sigma=2), control=list(fnscale=-1)) res$value res$par # Plot the fitted result over the data to see if it looks like a good fit (which we know it will be). x <- seq(from=min(dat$x),to=max(dat$x),by=0.1) y <- res$par['a']+(res$par['b']*x) plot(dat$x, dat$y, type='p') points(x, y, type='l', col='red') ### And compare to LM lmfit <- lm(y ~ x, dat) summary(lmfit)$coeff res$par
