极大似然估计（Fitting a Model by Maximum Likelihood）

#拟合正态分布
产生一个正态分布的数据
```{r}
set.seed(1001)
N <- 100
x <- rnorm(N, mean = 3, sd = 2)
mean(x)
sd(x)
```

定义似然函数
```{r}
LL <- function(mu, sigma) {
  R = dnorm(x, mu, sigma)
  -sum(log(R))
}
```

```{r}
library(stats4)
#由于计算时产生了负值，所以出现警告信息
mle(LL, start = list(mu = 1, sigma=1))
#比如,mean=1,sd=-1,标准差不能为负数
dnorm(x, 1, -1)
```


mle()调用optim(),mle()默认方法BFGS
```{r}
library(stats4)
mle(LL, start = list(mu = 1, sigma=1), method = "L-BFGS-B", lower = c(-Inf, 0),
    upper = c(Inf, Inf))
```

修改似然函数
```{r}
LL <- function(mu, sigma) {
  R = suppressWarnings(dnorm(x, mu, sigma))
  -sum(log(R))
}
mle(LL, start = list(mu = 1, sigma=1))
#初始值选择错误，将对结果又较大影响
mle(LL, start = list(mu = 0, sigma=1))
```

#拟合线性模型
产生数据
```{r}
x <- runif(N)
y <- 5 * x + 3 + rnorm(N)
fit <- lm(y ~ x)
summary(fit)
plot(x, y)
abline(fit, col = "red")
```

由于模型不是概率密度函数（pdf），无法像正态分布那样进行精确的似然函数的定义。由于线性模型要求残差符合正态分布，所以定义残差的似然函数
```{r}
LL <- function(beta0, beta1, mu, sigma) {
    # Find residuals
    R = y - x * beta1 - beta0
    # Calculate the likelihood for the residuals (with mu and sigma as parameters)
    R = suppressWarnings(dnorm(R, mu, sigma))
    # Sum the log likelihoods for all of the data points
    -sum(log(R))
}

LL <- function(beta0, beta1, mu, sigma) {
    R = y - x * beta1 - beta0
    #
    R = suppressWarnings(dnorm(R, mu, sigma, log = TRUE))
    #
    -sum(R)
}
```

```{r}
#同样，初始值很重要，错误的选择将有问题
#fit <- mle(LL, start = list(beta0 = 3, beta1 = 1, mu = 0, sigma=1))
#fit <- mle(LL, start = list(beta0 = 4, beta1 = 2, mu = 0, sigma=1))
fit <- mle(LL, start = list(beta0 = 5, beta1 = 3, mu = 0, sigma=1))
summary(fit)
fit <- mle(LL, start = list(beta0 = 2, beta1 = 1.5, sigma=1), fixed = list(mu = 0),
             nobs = length(y))
summary(fit)
```

模型比较
```{r}
AIC(fit)
BIC(fit)
logLik(fit)
```

mle()函数调用失败的问题在于对求Hessian 矩阵的逆矩阵。stats4包中的mle()方法除了有较好的初始值外，没有更好的解决办法。bbmle包中的mle2()提供了更好的替代方法。
```{r}
library(bbmle)
fit <- mle2(LL, start = list(beta0 = 3, beta1 = 1, mu = 0, sigma = 1))
summary(fit)
```

参考：http://www.r-bloggers.com/fitting-a-model-by-maximum-likelihood/