mlFun <- function(par) {
  xbeta <- X %*% par[1:K]
  Sig <- par[(K+1)]
  ll <- -(1/2) * log(2*pi) - (1/2) * log(Sig^2) - (1/(2*Sig^2))*(Y-xbeta)^2
  return(sum(ll))
}
n <- 100
K <- 1                                  #解释变量个数
paras <- c(1,1)                         #初值（斜率系数，扰动项方差)
X <- matrix(rnorm(n),ncol=K)
Y <- matrix(2*X+rnorm(n),ncol=1)
mlfit <- optim(par=paras,fn=mlFun)
mlfit$par
[1] 7.906692e-01 1.884238e-52