More Distributions for Modeling Claim

In this section, we will discuss two more distributions that are commonly used to model the insurance claims, which are the Pareto distribution and the Weibull distribution.

Pareto distribution

Characteristics of the Pareto distribution

If a random variable $X$ follows a Pareto distribution, we denote it as

$$X \sim \text{Pareto}(\alpha, \lambda),$$

where $\alpha$ and $\lambda$ are the shape and scale parameters, respectively.

• PDF:

  $$f(x)=\frac{\alpha\lambda^{\alpha}}{(\lambda+x)^{\alpha+1}},\mathrm{for}\ x>0$$

  The PDF can be derived by the following mixture distribution:

  Suppose $X|\theta \sim \text{exp}(\theta)$ and $\theta \sim \text{gamma}(\alpha, \lambda)$, and the density function of $\lambda$ is given by:

  $$f(\theta|\alpha, \lambda) = \frac{\lambda^\alpha}{\Gamma(\alpha)} \theta^{\alpha - 1} e^{-\lambda \theta}, \quad \theta > 0$$

  Then, the density function of $X$ is given by:

  $$\begin{align*} f_X(x|\alpha, \lambda) &= \int_0^\infty f(x|\theta) f(\theta|\alpha, \lambda) d\theta \\ &= \int_0^\infty (\theta e^{-\theta x}) \left[ \frac{\lambda^\alpha}{\Gamma(\alpha)} \theta^{\alpha - 1} e^{-\lambda \theta} \right] d\theta \\ &= \int_0^\infty \frac{\lambda^\alpha}{\Gamma(\alpha)} \theta^\alpha e^{-(\lambda + x)\theta} d\theta \\ &= \frac{\lambda^\alpha}{\Gamma(\alpha)} \int_0^\infty \theta^\alpha e^{-(\lambda + x)\theta} d\theta \\ &= \frac{\lambda^\alpha \Gamma(\alpha+1)}{\Gamma(\alpha)} \int_0^\infty \frac{1}{\Gamma(\alpha+1)} \theta^\alpha e^{-(\lambda + x)\theta} d\theta \\ &= \frac{\lambda^\alpha \Gamma(\alpha+1)}{\Gamma(\alpha)(\lambda + x)^{\alpha+1}} \underbrace{\left( \int_0^\infty \frac{(\lambda + x)^{\alpha+1}}{\Gamma(\alpha+1)} \theta^\alpha e^{-(\lambda + x)\theta} d\theta \right)}_{1} \\ &= \frac{\lambda^\alpha \Gamma(\alpha+1)}{\Gamma(\alpha)(\lambda + x)^{\alpha+1}} = \frac{\alpha \lambda^\alpha}{(\lambda + x)^{\alpha+1}} \end{align*}$$

• CDF:

  $$F(x)=1-\left(\frac{\lambda}{\lambda+x}\right)^{\alpha},\mathrm{for}\ x>0$$

• Survival function:

  $$\bar{F}(x)=\left(\frac{\lambda}{\lambda+x}\right)^{\alpha},\mathrm{for}\ x>0$$

  The survival function can be derived by the following:

  $$\begin{align*} \bar{F}(x|\alpha, \lambda) &= \int_0^\infty F(x|\theta) f(\theta|\alpha, \lambda) d\theta \\ &= \int_0^\infty e^{-\theta x} \left[ \frac{\lambda^\alpha}{\Gamma(\alpha)} \theta^{\alpha - 1} e^{-\lambda \theta} \right] d\theta \\ &= \frac{\lambda^\alpha}{\Gamma(\alpha)} \int_0^\infty \theta^{\alpha - 1} e^{-(\lambda + x)\theta} d\theta \\ &= \frac{\lambda^\alpha \Gamma(\alpha)}{\Gamma(\alpha)} \int_0^\infty \frac{1}{\Gamma(\alpha)} \theta^{\alpha-1} e^{-(\lambda + x)\theta} d\theta \\ &= \frac{\lambda^\alpha \Gamma(\alpha)}{\Gamma(\alpha)(\lambda + x)^\alpha} \underbrace{\left( \int_0^\infty \frac{(\lambda + x)^\alpha}{\Gamma(\alpha)} \theta^{\alpha - 1} e^{-(\lambda + x)\theta} d\theta \right)}_{1} \\ &= \frac{\lambda^\alpha \Gamma(\alpha)}{\Gamma(\alpha)(\lambda + x)^\alpha} = \frac{\lambda^\alpha}{(\lambda + x)^\alpha} = \left( \frac{\lambda}{\lambda + x} \right)^\alpha \end{align*}$$

Propositions of the Pareto distribution

• Mean:

  $$E(X) = \frac{\lambda}{\alpha - 1}, \quad \alpha > 1$$

  Following the previous setting, we can have:

  $$\begin{align*} E(X) &= \int_0^\infty E(X|\theta) f(\theta|\alpha, \lambda) d\theta \\ &= \int_0^\infty \frac{1}{\theta} f(\theta|\alpha, \lambda) d\theta \\ &= \int_0^\infty \frac{1}{\theta} \frac{\lambda^\alpha}{\Gamma(\alpha)} \theta^{\alpha - 1} e^{-\lambda \theta} d\theta \\ &= \frac{\lambda^\alpha}{\Gamma(\alpha)} \int_0^\infty \theta^{\alpha - 2} e^{-\lambda \theta} d\theta \\ &= \frac{\lambda^\alpha \Gamma(\alpha - 1)}{\Gamma(\alpha)} \int_0^\infty \frac{1}{\Gamma(\alpha - 1)} \theta^{\alpha - 2} e^{-\lambda \theta} d\theta \\ &= \frac{\lambda^\alpha \Gamma(\alpha - 1)}{\Gamma(\alpha) \lambda^{\alpha - 1}} \underbrace{\left( \int_0^\infty \frac{\lambda^{\alpha - 1}}{\Gamma(\alpha - 1)} \theta^{\alpha - 2} e^{-\lambda \theta} d\theta \right)}_{1} \\ &= \frac{\lambda^\alpha \Gamma(\alpha - 1)}{\Gamma(\alpha) \lambda^{\alpha - 1}} = \frac{\lambda}{\alpha - 1} \end{align*}$$

• K-th origin moment:

  $$E(X^k) = \frac{\lambda^k \Gamma(\alpha + k)}{\Gamma(\alpha) (\alpha - 1)^k}, \quad \alpha > 1$$

  Derive the $k$-th origin moment

  $$\begin{align*} E(X^k) &= \int_0^\infty E(X^k|\theta) f(\theta|\alpha, \lambda) d\theta \\ &= \int_0^\infty \frac{\Gamma(k+1)}{\theta^k} f(\theta|\alpha, \lambda) d\theta \\ &= \int_0^\infty \frac{\Gamma(k+1)}{\theta^k} \frac{\lambda^\alpha}{\Gamma(\alpha)} \theta^{\alpha - 1} e^{-\lambda \theta} d\theta \\ &= \frac{\lambda^\alpha \Gamma(k+1)}{\Gamma(\alpha)} \int_0^\infty \theta^{\alpha - k - 1} e^{-\lambda \theta} d\theta \\ &= \frac{\lambda^\alpha \Gamma(k+1) \Gamma(\alpha - k)}{\Gamma(\alpha)} \int_0^\infty \frac{1}{\Gamma(\alpha - k)} \theta^{\alpha - k - 1} e^{-\lambda \theta} d\theta \\ &= \frac{\lambda^\alpha \Gamma(k+1) \Gamma(\alpha - k)}{\Gamma(\alpha) \lambda^{\alpha - k}} \underbrace{\left( \int_0^\infty \frac{\lambda^{\alpha - k}}{\Gamma(\alpha - k)} \theta^{\alpha - k - 1} e^{-\lambda \theta} d\theta \right)}_{1} \\ &= \frac{\lambda^\alpha \Gamma(k+1) \Gamma(\alpha - k)}{\Gamma(\alpha) \lambda^{\alpha - k}} = \frac{\lambda^k \Gamma(k+1) \Gamma(\alpha - k)}{\Gamma(\alpha)} = \frac{\lambda^k k! (\alpha - k - 1)!}{(\alpha - 1)!} \\ &= \frac{\lambda^k k!}{(a-1) \cdots (a-k)}, \quad -1 < k < \alpha \end{align*}$$

• Variance:

  $$\text{Var}(X) = \frac{\lambda^2}{(\alpha - 1)^2 (\alpha - 2)}, \quad \alpha > 2$$

  Derive the variance of Pareto distribution by standard formula

  $$\begin{align*} \text{Var}(X) &= E(X^2) - E^2(X) = \frac{\lambda^2 \Gamma(\alpha + 2)}{\Gamma(\alpha) (\alpha - 1)^2} - \left( \frac{\lambda}{\alpha - 1} \right)^2 \\ &= \frac{\lambda^2 \Gamma(\alpha + 2)}{\Gamma(\alpha) (\alpha - 1)^2} - \frac{\lambda^2}{(\alpha - 1)^2} \\ &= \frac{\lambda^2 \Gamma(\alpha + 2)}{\Gamma(\alpha) (\alpha - 1)^2} - \frac{\lambda^2 \Gamma(\alpha + 1)}{\Gamma(\alpha) (\alpha - 1)^2} \\ &= \frac{\lambda^2 \Gamma(\alpha + 2) - \lambda^2 \Gamma(\alpha + 1)}{\Gamma(\alpha) (\alpha - 1)^2} \\ &= \frac{\lambda^2 \Gamma(\alpha + 2 - 1) (\alpha + 1 - 1) - \lambda^2 \Gamma(\alpha + 1)}{\Gamma(\alpha) (\alpha - 1)^2} \\ &= \frac{\lambda^2 \Gamma(\alpha + 2)}{\Gamma(\alpha) (\alpha - 1)^2} - \frac{\lambda^2 \Gamma(\alpha + 1)}{\Gamma(\alpha) (\alpha - 1)^2} \\ &= \frac{\lambda^2 \Gamma(\alpha + 2) - \lambda^2 \Gamma(\alpha + 1)}{\Gamma(\alpha) (\alpha - 1)^2} \\ &= \frac{\lambda^2 \Gamma(\alpha + 2 - 1) (\alpha + 1 - 1) - \lambda^2 \Gamma(\alpha + 1)}{\Gamma(\alpha) (\alpha - 1)^2} \\ &= \frac{\lambda^2 \Gamma(\alpha + 1) (\alpha + 1 - 1) - \lambda^2 \Gamma(\alpha + 1)}{\Gamma(\alpha) (\alpha - 1)^2} \\ &= \frac{\lambda^2 \Gamma(\alpha + 1)}{\Gamma(\alpha) (\alpha - 1)^2} = \frac{\lambda^2}{(\alpha - 1)^2 (\alpha - 2)} \end{align*}$$

Estimation of parameters

Method of moments (MM)

In the Pareto distribution, we have two parameters to estimate, so we need to use the method of moments to estimate the parameters $\alpha$ and $\lambda$. The first moment of the Pareto distribution is $\frac{\lambda}{\alpha - 1}$, and the second moment is $\frac{\lambda^2}{(\alpha - 1)^2 (\alpha - 2)}$. Both of them should be equal to the sample moments.

We can solve the following criteria:

$$E(X) = \frac{\lambda}{\alpha - 1} = \frac{1}{n} \sum_{i=1}^n X_i = \bar{X}$$

$$Var(X) = \frac{\lambda^2}{(\alpha - 1)^2 (\alpha - 2)} = \frac{1}{n} \sum_{i=1}^n (X_i - \bar{X})^2 = S^2$$

Solve the first moment,

$$\hat{\lambda} = \bar{X}(\hat{\alpha} - 1)$$

Then, replace $\lambda$ with $\bar{X}(\hat{\alpha} - 1)$ in the second moment, we can get the estimator of $\alpha$:

$$\begin{align*} \frac{\bar{X}(\hat{\alpha} - 1)}{(\hat{\alpha} - 1)^2 (\hat{\alpha} - 2)} & = S^2 \\ \bar{X}(\hat{\alpha} - 1) & = S^2 (\hat{\alpha} - 1)^2 (\hat{\alpha} - 2) \\ \bar{X} & = S^2 (\hat{\alpha} - 1) (\hat{\alpha} - 2) \\ \bar{X} & = S^2 (\hat{\alpha}^2 - 3\hat{\alpha} + 2) \\ \hat{\alpha}^2 - 3\hat{\alpha} + 2 & = \frac{\bar{X}}{S^2} \\ \hat{\alpha}^2 - 3\hat{\alpha} + 2 - \frac{\bar{X}}{S^2} & = 0 \\ \end{align*}$$

According to the standard solution, we can solve the equation:

$$\hat{\alpha} = \frac{3 \pm \sqrt{5 + 4\frac{\bar{X}}{S^2}}}{2}$$

Then, we can solve the estimator of $\lambda$:

$$\begin{align*} \hat{\lambda} & = \bar{X}(\hat{\alpha} - 1) \\ & = \bar{X}(\frac{3 \pm \sqrt{5 + 4\frac{\bar{X}}{S^2}}}{2} - 1) \\ & = \bar{X}(\frac{1 \pm \sqrt{5 + 4\frac{\bar{X}}{S^2}}}{2}) \end{align*}$$

The estimator of $\alpha$ and $\lambda$ are:

$$\hat{\alpha} = \frac{3 \pm \sqrt{5 + 4\frac{\bar{X}}{S^2}}}{2}, \quad \hat{\lambda} = \bar{X}(\frac{1 \pm \sqrt{5 + 4\frac{\bar{X}}{S^2}}}{2})$$

Maximum likelihood estimation (MLE)

The likelihood function of the Pareto distribution is given by:

$$L(\alpha, \lambda) = \prod_{i=1}^n \frac{\alpha\lambda^{\alpha}}{(\lambda+x_i)^{\alpha+1}}$$

Then the log-likelihood function is:

$$\begin{array}{rl} l & = \sum_{i=1}^n \left( \log(\alpha) + \alpha \log(\lambda) - (\alpha + 1) \log(\lambda + x_i) \right) \\ \\ & = n \log(\alpha) + n \alpha \log(\lambda) - (\alpha + 1) \sum_{i=1}^n \log(\lambda + x_i) \end{array}$$

Taking the derivative of the log-likelihood function with respect to $\alpha$ and $\lambda$,

we have the following equations:

$$\frac{dl}{d\alpha} = \frac{n}{\alpha} + n \log(\lambda) - \sum_{i=1}^n \log(\lambda + x_i) = 0$$

Solve the equation, we have the MLE of $\alpha$:

$$\begin{align*} \frac{n}{\alpha} & = \sum_{i=1}^n \log(\lambda + x_i) - n \log(\lambda) \\ \hat{\alpha} & = \frac{n}{\sum_{i=1}^n \log(\lambda + x_i) - n \log(\lambda)} \\ & = \frac{n}{\sum_{i=1}^n \log(1+ \frac{x_i}{\lambda})} \\ \end{align*}$$

Then, taking the derivative of the log-likelihood function with respect to $\lambda$,

$$\frac{dl}{d\lambda} = \frac{n\alpha}{\lambda} -(\alpha + 1)\sum_{i=1}^n \frac{1}{\lambda + x_i} = 0$$

Solve the equation, we have the MLE of $\lambda$:

$$\begin{align*} \frac{n\alpha}{\lambda} & = (\alpha + 1)\sum_{i=1}^n \frac{1}{\lambda + x_i} \\ (\frac{n}{\lambda} - \sum_{i=1}^n \frac{1}{\lambda + x_i})\alpha & = \sum_{i=1}^n \frac{1}{\lambda + x_i}\\ (\sum_{i=1}^n (\frac{1}{\lambda} - \frac{1}{\lambda + x_i})) \alpha & = \sum_{i=1}^n \frac{1}{\lambda + x_i}\\ (\sum_{i=1}^{n}\frac{x_i}{\lambda(\lambda + x_i)}) \alpha & = \sum_{i=1}^{n} \frac{1}{\lambda + x_i} \\ \hat{\alpha} & = \frac{\sum_{i=1}^{n} \frac{1}{\lambda + x_i}}{(\sum_{i=1}^{n}\frac{x_i}{\lambda(\lambda + x_i)})} \end{align*}$$

Then, we can have the equation:

$$\frac{n}{\sum_{i=1}^n \log(1+ \frac{x_i}{\lambda})} -\frac{\sum_{i=1}^{n} \frac{1}{\lambda + x_i}}{(\sum_{i=1}^{n}\frac{x_i}{\lambda(\lambda + x_i)})} = 0$$

x <- rpareto(1000, shape = 2, scale = 1)

f <- function(lambda) {
  n / sum(log(1 + x / lambda)) - sum(1 / (lambda + x)) / sum(x / (lambda *
(lambda + x)))
}

lambda <- uniroot(f, c(1e-8, 1e8))$root

alpha <- length(x) / sum(log(1 + x / lambda))

paste0("alpha:", alpha)
paste0("lambda:", lambda)

Weibull distribution

Characteristics of the Weibull distribution

If a random variable $X$ follows a Weibull distribution, we denote it as

$$X \sim \text{Weibull}(c, \gamma),$$

where $c$ and $\gamma$ are the shape and scale parameters, respectively.

• PDF:

  $$f_X(x) = c \gamma x^{\gamma - 1} e ^{-c x ^{\gamma}}$$

  Weibull distribution can be derived from the exponential distribution.

  Suppose $Y \sim exp(c)$, then $X = Y^{1/\gamma} \sim Weibull(c, \gamma)$

  First, recall the change of the variable formula:

  $$f_X(x) = \{ \begin{align*} f_Y(g^{-1}(x)), & \ \text{if} \ x > 0 \\ 0 &, \text{otherwise} \end{align*}$$

  Then, recall the density function of the exponential distribution:

  $$f_Y(y) = ce^{-cy}, \quad \text{for} \quad y > 0$$

  Then, we have

  $$X = g(Y) = Y^{1 / \gamma}$$

  The inverse function of G is given by

  $$Y = g^{-1}(x) = x^{\gamma}$$

  Therefore, the density function of X is given by:

  $$\begin{align*} f_X(x) & = F_Y(g^{-1}(x)) | \frac{dg^{-1}}{dx} \\ & = ce^{-cx^{\gamma}}\gamma x^{\gamma -1} \\ & = c \gamma x^{\gamma -1} e^{-cx^{\gamma}} \end{align*}$$

• CDF:

  $$F_X(x) = P[X\leq x] = P[Y^{1/\gamma} \leq x] = P[Y \leq x^{\gamma}] = 1 - e^{-cx^{\gamma}}$$

• Survival function:

  $$\bar{F}_X(x) = P[X > x] = P[Y^{1/\gamma} > x] = P[Y > x^{\gamma}] = e^{-cx^{\gamma}}$$

Propositions of the Weibull distribution

• Mean:

  $$E(X) = E(Y^{1/\gamma}) = \frac{\Gamma(1 + 1/\gamma)}{c^{1/\gamma}} = \Gamma(1 + 1/\gamma) c^{\gamma}$$

  Recall the expectation of the exponential distribution:

  $$E(Y^{k}) = \frac{k!}{c^k} = \frac{\Gamma(1 + k)}{c^k}$$

  Then, extend to the $k$-th origin moment:

  $$E(X^k) = E(Y^{k/\gamma}) = \frac{\Gamma(1+k/\gamma)}{c^{k/\gamma}} = \Gamma(1 + k/\gamma) c^{-k/\gamma}$$

• Variance:

  $$\begin{align*} Var(X) & = E(X^2) - E^2(X)\\ & = \Gamma(1 + 2/\gamma) c^{-2/\gamma} - (\Gamma(1 + 1/\gamma) c^{-1/\gamma})^2 \\ & = \Gamma(1 + 2/\gamma) c^{-2/\gamma} - (\Gamma(1 + 1/\gamma)^2 c^{-2/\gamma}) \\ & = \Gamma(1 + 2/\gamma) c^{-2/\gamma} - \Gamma^2(1 + 1/\gamma) c^{-2/\gamma} \\ & = c^{-2/\gamma} (\Gamma(1 + 2/\gamma) - \Gamma^2(1 + 1/\gamma)) \end{align*}$$

Estimation of parameters

Method of moments (MM)

The first two moments of the Weibull distribution are:

$$\begin{align*} \mu_1 & = E(X) = \Gamma(1 + 1/\gamma) c^{\gamma} \\ \mu_2 & = Var(X) = c^{-2/\gamma} (\Gamma(1 + 2/\gamma) - \Gamma^2(1 + 1/\gamma)) \end{align*}$$

Solve the first equation with taking logarithem:

$$\begin{align*} \log(\mu_1) & = \log(\Gamma(1 + 1/\gamma)) + \gamma \log(c) \\ \rightarrow \gamma & = \frac{\log(\mu_1) - \log(\Gamma(1 + 1/\gamma))}{\log(c)} \end{align*}$$

Then, we can get the estimator of $c$ in the second equation:

$$\begin{align*} \log(\mu_2) & = -\frac{2}{\gamma} \log(c) + \log(\Gamma(1 + 2/\gamma)) - \log(\Gamma^2(1 + 1/\gamma)) \\ \frac{2}{\gamma} \log(c) & = \log(\Gamma(1 + 2/\gamma)) - \log(\Gamma^2(1 + 1/\gamma)) - log(\mu_2) \\ \log(c) & = \frac{\log(\Gamma(1 + 2/\gamma)) - \log(\Gamma^2(1 + 1/\gamma)) - \log(\mu_2)}{2/\gamma} \\ c & = \exp\left(\frac{\log(\Gamma(1 + 2/\gamma)) - \log(\Gamma^2(1 + 1/\gamma)) - \log(\mu_2)}{2/\gamma}\right) \\ \end{align*}$$

Then, we replace $c = \exp\left(\frac{\log(\Gamma(1 + 2/\gamma)) - \log(\Gamma^2(1 + 1/\gamma)) - \log(\mu_2)}{2/\gamma}\right)$ in the first equation, we can get the estimator of $\gamma$.

$$\gamma = \frac{\log(\mu_1) - \log(\Gamma(1 + 1/\gamma))}{\log\left(\exp\left(\frac{\log(\Gamma(1 + 2/\gamma)) - \log(\Gamma^2(1 + 1/\gamma)) - \log(\mu_2)}{2/\gamma}\right)\right)}$$

Next, we replace the first and second moments with the sample moments:

$$\begin{align*} \mu_1 & = \frac{1}{n} \sum_{i=1}^n x_i = \bar{X} \\ \mu_2 & = Var(X) = \frac{\sum_{i=1}^n (x_i - \bar{x})^2}{n} = S^2 \end{align*}$$

Then, we can estimate the parameters $\gamma$ and $c$ by the following equations:

$$\begin{align*} \hat{\gamma} & = \frac{\log(\bar{X}) - \log(\Gamma(1 + 1/\hat{\gamma}))}{\log\left(\exp\left(\frac{\log(\Gamma(1 + 2/\hat{\gamma})) - \log(\Gamma^2(1 + 1/\hat{\gamma})) - \log(S^2)}{2/\hat{\gamma}}\right)\right)} \\ \hat{c} & = \exp\left(\frac{\log(\Gamma(1 + 2/\hat{\gamma})) - \log(\Gamma^2(1 + 1/\hat{\gamma})) - \log(S^2)}{2/\hat{\gamma}}\right) \end{align*}$$

Finally, we can solve the equation by R or Python.

# generate 1000 Weibull random variables
x <- rweibull(1000, shape = 2, scale = 1)

aux <- log(mean(x))

f <- function(gamma) {
  (log(mean(x)) - log(gamma(1 + 1 / gamma))) / log(exp((log(gamma(1 + 2 / gamma)) - log(gamma(1 + 1 / gamma))^2 - log(var(x))) / (2 / gamma))) - gamma
}

gamma <- uniroot(f, c(1e-8, 1e8))$root

c <- exp((log(gamma(1 + 2 / gamma)) - log(gamma(1 + 1 / gamma))^2 - log(var(x))) / (2 / gamma))

cat(
"Estimated shape parameter: ", gamma, "\n",
"Estimated scale parameter: ", c, "\n"
)

Maximum likelihood estimation (MLE)

The likelihood functon of Weibull distribution is given by:

$$L = \prod_{i = 1}^{n} c \gamma x_i^{y-1} e^{-cx_i^{\gamma}}$$

The log-likelihood function will be:

$$\begin{align*} l & = \sum_{i = 1}^n\log(c) +\sum_{i = 1}^n\log(\gamma) + (\gamma - 1) \sum_{i = 1}^n - c\sum_{i = 1}^n x_i^{\gamma} \\ & =n \log(c) + n \log(\gamma) + (\gamma - 1) \sum_{i = 1}^n \log(x+i) - c\sum_{i = 1}^n x_i^{\gamma} \end{align*}$$

Taking the derivative of the log-likelihood function with respect to $c$ and $\gamma$, we have the following equations:

$$\frac{dl}{dc} = \frac{n}{c} - \sum_{i = 1}^n x_i^{\gamma} = 0$$

Solve the equation, we have the MLE of $c$:

$$\frac{n}{c} = \sum_{i = 1}^n x_i^{\gamma} \Rightarrow \hat{c} = \frac{n}{\sum_{i = 1}^n x_i^{\gamma}}$$

Then, taking the derivative of the log-likelihood function with respect to $\gamma$,

$$\frac{dl}{d\gamma} = \frac{n}{\gamma} + \sum_{i = 1}^n \log(x_i) - c\sum_{i = 1}^n x_i^{\gamma} \log(x_i) = 0$$

Replace $c$ with $\frac{n}{\sum_{i = 1}^n x_i^{\gamma}}$, we can get the MLE of $\gamma$:

$$\begin{align*} \frac{n}{\gamma} + \sum_{i = 1}^n \log(x_i) - \frac{n}{\sum_{i = 1}^n x_i^{\gamma}}\sum_{i = 1}^n x_i^{\gamma} \log(x_i) & = 0 \\ \frac{1}{\gamma} + \overline{\log(x_i)} - \frac{1}{\sum_{i = 1}^n x_i^{\gamma}}\sum_{i = 1}^n x_i^{\gamma} \log(x_i) & = 0 \\ \end{align*}$$

This equation can be solved by R or Python.

# generate 1000 Weibull random variables
x <- rweibull(1000, shape = 2, scale = 1)

aux <- mean(log(x))

f <- function(gamma) {
1 / gamma + aux - sum(x^gamma * log(x)) / sum(x^gamma)
}

gamma <- uniroot(f, c(1e-8, 1e8))$root
c <- length(x) / sum(x^gamma)

cat(
"Estimated shape parameter: ", gamma, "\n",
"Estimated scale parameter: ", c, "\n"
)