Log-likelihood:

```{r}
ll <- function(beta, delta, t, x, u) {
	n <- length(t)
	z <- (log(t) - beta[1] - beta[2] * x - beta[3] * u) / delta
	-n * log(delta) + sum(z - exp(z))
}
```

Data:

```{r}
data(leuk, package = "MASS")
leuk <- transform(leuk, u = ifelse(ag == "present", 1, 0), x = log(wbc / 10000))
```

Using $\delta$:

```{r}
f <- function(par)
    -ll(par[1:3], par[4], leuk$time, leuk$x, leuk$u)
par0 <- c(log(mean(leuk$time)), 0, 0, 1)
v <- optim(par0, f, method = "BFGS", hessian = TRUE)
v$convergence
v$counts
(par <- v$par)
(se <- sqrt(diag(solve(v$hessian))))
par0 <- c(coef(lm(log(time) ~ x + u, data = leuk)), 1)
v <- optim(par0, f, method = "BFGS", hessian = TRUE)
v$convergence
v$counts
(par <- v$par)
(se <- sqrt(diag(solve(v$hessian))))
```

Using $\log\delta$:

```{r}
f_ld <- function(par)
    -ll(par[1:3], exp(par[4]), leuk$time, leuk$x, leuk$u)
par0_ld <- c(log(mean(leuk$time)), 0, 0, 0)
v <- optim(par0_ld, f_ld, method = "BFGS", hessian = TRUE)
v$convergence
v$counts
par <- v$par
se <- sqrt(diag(solve(v$hessian)))
c(par[1:3], exp(par[4]))
c(se[1:3], se[4] * exp(par[4]))
```

```{r}
llg <- function(beta, delta, t, x, u) {
	n <- length(t)
	z <- (log(t) - beta[1] - beta[2] * x - beta[3] * u) / delta
    exp_z <- exp(z)
    exp_z_m_1 <- exp_z - 1
    c(sum(exp_z_m_1) / delta,
      sum(exp_z_m_1 * x) / delta,
      sum(exp_z_m_1 * u) / delta,
      (sum(exp_z * z - z) - n) / delta)
}
fg <- function(par) - llg(par[1:3], par[4], leuk$t, leuk$x, leuk$u)
par0 <- c(log(mean(leuk$time)), 0, 0, 1)
v <- optim(par0, f, fg, method = "BFGS", hessian = TRUE)
v$convergence
v$counts
(par <- v$par)
(se <- sqrt(diag(solve(v$hessian))))
```

```{r}
fg_ld <- function(par) {
    g <- - llg(par[1:3], exp(par[4]), leuk$t, leuk$x, leuk$u)
    g[4] <- exp(par[4]) * g[4]
    g
}
v <- optim(par0_ld, f_ld, fg_ld, method = "BFGS", hessian = TRUE)
v$convergence
par <- v$par
se <- sqrt(diag(solve(v$hessian)))
c(par[1:3], exp(par[4]))
c(se[1:3], se[4] * exp(par[4]))
```
