---
title: "Fit MLEs"
author: "Prof. Eric A. Suess"
format: 
  html:
    embed-resources: true
toc: true
---

# Exercises

## Exercise 1

### a. 

$$f(x|\theta) = \theta x^{\theta−1}$$

where $0 < x < 1$, $0 < \theta$, and 0 otherwise.

The model is $Beta(\theta, 1)$, where $\theta$ is the parameter of interest.

```{r}
library(fitdistrplus)

n <- 100
theta <- 0.5

x <- rbeta(n, theta, 1)

plotdist(x)

fit <- fitdist(x, "beta", start = list(shape1 = 1), fix.arg = list(shape2 = 1))

summary(fit)

llplot(fit, pal.col = heat.colors(100), fit.show = TRUE)

fit$estimate

fit$sd

confint(fit)
```

Compare with the math.

Estimate

```{r}
theta.hat <- -mean(log(x))^{-1}
theta.hat
```

Estimate AV

```{r}
av.hat <- theta.hat^2/n
av.hat

se.hat <- sqrt(av.hat)
se.hat

confint <- c(theta.hat - 1.96*se.hat, theta.hat + 1.96*se.hat)
confint
```

### b. 

$$f(x|\theta) = (\theta + 1) x^{−\theta−2}$$

where $1 < x$, $0 < \theta$, and 0 otherwise.

The model is $Pareto( (\theta+1), 1)$, where $\theta$ is the parameter of interest.

```{r}
library(fitdistrplus)
library(Pareto)

n <- 100
theta <- 5

x <- rPareto(n, t = 1, alpha = (theta + 1))

plotdist(x)

fit <- fitdist(x, dPareto, start = list(alpha = 10), fix.arg = list(t = 1))

summary(fit)

llplot(fit, pal.col = heat.colors(100), fit.show = TRUE)

fit$estimate 

fit$sd

confint(fit)

```

Compare with the math.

Note: $\theta = \alpha - 1$.  So subtract 1 from the MLE and from the Confidence Interval to the the MLE and CI for $\theta$.

Estimate

```{r}
theta.hat <- (n/sum(log(x))) - 1

theta.hat

```

Estimate AV

```{r}
av.hat <- ((theta.hat+1)^2)/n
av.hat

se.hat <- sqrt(av.hat)
se.hat

confint <- c(theta.hat - 1.96*se.hat, theta.hat + 1.96*se.hat) 
confint
```

### c.

$$f(x|\theta) = \theta^{2} x e^{−\theta x}$$

where $0 < x$, $0 < \theta$, and 0 otherwise.

The model is $Gamma(2, \theta)$, where $\theta$ is the parameter of interest.

```{r}
library(fitdistrplus)

n <- 100
theta <- 5

x <- rgamma(n, shape = 2, rate = theta)

plotdist(x)

fit <- fitdist(x, "gamma", start = list(rate = 2), fix.arg = list(shape = 2))

summary(fit)

llplot(fit, pal.col = heat.colors(100), fit.show = TRUE)

fit$estimate

fit$sd

confint(fit)


```

Compare with the math.

Estimate


```{r}

theta.hat <- 2/mean(x)

theta.hat

```

Estimate AV

```{r}
av.hat <- (theta.hat^2)/(2*n)
av.hat

se.hat <- sqrt(av.hat)
se.hat

confint <- c(theta.hat - 1.96*se.hat, theta.hat + 1.96*se.hat)
confint
```






 


### d.

$$f(x|\theta) = \theta(1 − \theta)^{x−1}$$

for $x = 1, 2,...$, $0 < \theta < 1$, and 0 otherwise

Note that the model given is for 

$X =$ number of trials until the first success

$$
f(x|\theta) = \theta (1-\theta)^{x-1}
$$

where $x = 1, 2, 3, \ldots$.

But the model in R is for

$Y = X - 1 =$ number of failures before the first success


$$
f(x|\theta) = \theta (1-\theta)^{x}
$$

where $x = 0, 1, 2, 3, \ldots$.

The model is $Geometric(\theta)$, where $\theta$ is the parameter of interest.  Use the model in R.

```{r}
library(fitdistrplus)

n <- 100
theta <- 0.5

x <- rgeom(n, theta)

plotdist(x)

fit <- fitdist(x, "geom", start = list(prob = 0.5))

summary(fit)

llplot(fit, pal.col = heat.colors(100), fit.show = TRUE)

fit$estimate

fit$sd

confint(fit)
```

Compare with the math.

Since the model is different, we need to adjust the math.  

The MLE is $\hat{\theta } = \frac{1}{\bar{x} + 1}$.

Estimate

```{r}
theta.hat <- (mean(x) +  1)^{-1}
theta.hat
```

Estimate AV

```{r}
av.hat <- (theta.hat^2)*(1-theta.hat)/n
av.hat

se.hat <- sqrt(av.hat)
se.hat

confint <- c(theta.hat - 1.96*se.hat, theta.hat + 1.96*se.hat)
confint

```