{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ade1a501-2021-4f19-9a6b-46e83ceaa38d",
   "metadata": {},
   "source": [
    "# 🏆  Homework"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2043ee76",
   "metadata": {},
   "source": [
    "## Maximum Likelihood Estimation of Parameters of Gamma Distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd5a11be",
   "metadata": {},
   "source": [
    "## Answer by Muhammed Talha Kaya"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f5fcd5b-93df-406c-9097-4c2cad5dab69",
   "metadata": {},
   "source": [
    "Let's consider the maximum likelihood estimation for the parameters of Gamma distribution with unknown scale parameter $\\theta$ and shape parameter $\\kappa$ based on a random sample of size n."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8f90b2c-7efd-4ec1-9a3d-bdad6f56c4c5",
   "metadata": {},
   "source": [
    "The probability density function of Gamma distribution with unknown parameters ($\\theta, \\kappa$) is:\n",
    "\n",
    "$f(x;\\theta, \\kappa) = \\frac{1}{\\theta^\\kappa \\Gamma(\\kappa)}x^{\\kappa-1}e^{-\\frac{x}{\\theta}} $ for $x>0$ and $\\theta,\\kappa>0$.\n",
    "\n",
    "Then the likelihood and log-likelihood functions based on a random sample of size n are defined, respectively, as follows:\n",
    "\n",
    "$L(\\theta, \\kappa) = \\frac{1}{\\theta^{n\\kappa} \\Gamma(\\kappa)^n} \\big(\\prod_{i}^{n} x_i\\big)^{\\kappa-1}e^{-\\sum_{i}^{n} \\frac{x_i}{\\theta}} $ for $x>0$ and $\\theta,\\kappa>0$, and\n",
    "\n",
    "$\\log(L(\\theta, \\kappa)) = -n \\kappa\\log(\\theta)-n \\log(\\Gamma(\\kappa)) + (\\kappa-1)\\log(\\prod_{i}^{n} x_i)-\\frac{\\sum_{i}^{n} x_i}{\\theta}$.\n",
    "\n",
    "The partial derivatives are:\n",
    "\n",
    "$\\frac{\\partial\\log(L(\\theta, \\kappa))}{\\partial\\theta }= -\\frac{n \\kappa}{\\theta} + \\frac{\\sum_{i}^{n} x_i}{\\theta^2}$ and \n",
    "\n",
    "$\\frac{\\partial\\log(L(\\theta, \\kappa))}{\\partial\\kappa }= -n \\log(\\theta) - n \\frac{\\Gamma'(\\kappa)}{\\Gamma(\\kappa)} +\n",
    "\\log(\\prod_{i}^{n} x_i)$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "515b41d5-ab10-4760-9585-1006d4c656d3",
   "metadata": {},
   "source": [
    "Let $\\psi(\\kappa) = \\frac{\\Gamma'(\\kappa)}{\\Gamma(\\kappa)}$ denote the psi function and $\\tilde X = \\big(\\prod_{i}^{n} x_i\\big)^{\\frac{1}{n}}$, where $n \\log(\\tilde X) = \\log(\\prod_{i}^{n} x_i)$, denote the geometric mean of the sample, then set the derivatives equal to zero:\n",
    "\n",
    "$\\frac{\\partial\\log(L(\\theta, \\kappa))}{\\partial\\theta }= -\\frac{n \\kappa}{\\theta} + \\frac{\\sum_{i}^{n} x_i}{\\theta^2}=0 \\Rightarrow $ \n",
    "\n",
    "$\\hat{\\theta}=  \\frac{\\sum_{i}^{n} x_i}{n \\hat{\\kappa}} = \\frac{\\bar x}{\\hat{\\kappa}}$.\n",
    "\n",
    "So $\\hat{\\theta}_{MLE}$ depends on $\\hat{\\kappa}_{MLE}$. \n",
    "\n",
    "$\\frac{\\partial\\log(L(\\theta, \\kappa))}{\\partial\\kappa }= -n \\log(\\frac{\\bar x}{\\hat{\\kappa}}) - n \\psi(\\hat{\\kappa}) +\n",
    "n \\log(\\tilde X)\n",
    "=-n \\log(\\bar x) + n \\log(\\hat{\\kappa})-n \\psi(\\hat{\\kappa}) + n  \\log(\\tilde X) = \\log(\\hat{\\kappa})-\\psi(\\hat{\\kappa}) + \\log(\\frac{\\tilde X}{\\bar x})=0$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c402767-a89c-4b56-a9b5-b793b940fa4d",
   "metadata": {},
   "source": [
    "The last ML equation cannot be solved in closed form due to non-linear form of $\\kappa$. We need an iterative algorithm to find the root."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7629eb47-c153-4148-9fdc-fc7828017cf2",
   "metadata": {},
   "source": [
    "İbrahim Talha suggests to use Newton-Raphson algoritm to find the root of the last ML equation iteratively as follows:\n",
    "\n",
    "$\\hat{\\kappa}_{n+1} = \\hat{\\kappa}_{n}-\\frac{f(\\hat{\\kappa}_{n})}{f'(\\hat{\\kappa}_{n})}$,\n",
    "\n",
    "where $f(\\hat{\\kappa}_{n})= \\log(\\hat{\\kappa}_{n})-\\psi(\\hat{\\kappa}_{n}) + \\log(\\frac{\\tilde X}{\\bar x})$, $f'(\\kappa_{n})= \\frac{1}{\\hat{\\kappa}_{n}}-\\psi'(\\hat{\\kappa}_{n})$ with $\\psi(\\hat{\\kappa}_{n}) = \\frac{\\Gamma'(\\hat{\\kappa}_{n})}{\\Gamma(\\hat{\\kappa}_{n})}$ and $\\tilde X = \\big(\\prod_{i}^{n} x_i\\big)^{\\frac{1}{n}}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3e255c8f-e8b4-4c97-bc2a-1bf788eca664",
   "metadata": {},
   "outputs": [],
   "source": [
    "#import required libraries\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import scipy.stats as stats\n",
    "from scipy import special\n",
    "import math\n",
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b5c489b1-41dd-43f2-a1be-931b1aab574a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Initial Kappa = 1\n",
      "Function Value in Initial Kappa = 0.40741223354311584\n",
      "\n",
      "New Kappa = 1.6317114484805046\n",
      "Function Value in New Kappa = 0.16690439321216854\n",
      "\n",
      "New Kappa = 2.377444934540144\n",
      "Function Value in New Kappa = 0.05500775249943013\n",
      "\n",
      "New Kappa = 2.924951073208216\n",
      "Function Value in New Kappa = 0.01077193924582498\n",
      "\n",
      "New Kappa = 3.0907780282741752\n",
      "Function Value in New Kappa = 0.0006042910603489826\n",
      "\n",
      "New Kappa = 3.10121944491736\n",
      "Function Value in New Kappa = 2.1279091184656096e-06\n",
      "\n",
      "Theta Hat = 1.9624247419732048\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'\\nsort_data = np.sort(data)\\n\\nplt.plot(sort_data, stats.gamma.pdf(sort_data, a = 2, scale = 2))\\n\\nplt.show() \\n'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "#special.polygamma(1, x)\n",
    "\n",
    "#Generate 10000 data with gamma distribution\n",
    "#np.random.gamma(Kappa,Theta,Size of Data)\n",
    "data = np.random.gamma(shape=2,scale=2,size=500)\n",
    "#we are estimating shape parameter\n",
    "\n",
    "#define x_tilda based on the data\n",
    "def x_tilda(data):\n",
    "\n",
    "\tx_tilda = 0\n",
    "\n",
    "\tfor item in data:\n",
    "\n",
    "\t\tx_tilda = x_tilda+ math.log(item)\n",
    "\n",
    "\treturn math.exp(x_tilda*(1/np.size(data)))\n",
    "\n",
    "x_tilda = x_tilda(data)\n",
    "\n",
    "x_bar = np.mean(data)\n",
    "\n",
    "#scipy.special.polygamma(n, x) = nth derivative of the psi(x) function\n",
    "def psi(k_hat): \n",
    "\t\n",
    "\tpsi_value = special.polygamma(0, k_hat)\n",
    "\n",
    "\treturn psi_value\n",
    "\n",
    "def fd_psi(k_hat): #first derivative of psi\n",
    "\t\n",
    "\tfd_psi_value = special.polygamma(1, k_hat)\n",
    "\n",
    "\treturn fd_psi_value\n",
    "\n",
    "\n",
    "def kappa_function(k_hat,x_tilda,x_bar): #F(Kappa Hat) =? 0 \n",
    "\n",
    "\tkappa_function_result = math.log(k_hat) - psi(k_hat) + math.log(x_tilda) - math.log(x_bar)\n",
    "\t#print(kappa_function_result)\n",
    "\treturn kappa_function_result\n",
    "\n",
    "def fd_kappa_function(k_hat): #first derivative of kappa_function\n",
    "\n",
    "\tfd_kappa_function_result = (1/k_hat) - fd_psi(k_hat)\n",
    "\t#print(fd_kappa_function_result)\n",
    "\treturn fd_kappa_function_result\n",
    "\n",
    "#Newthon-Raphson Method\n",
    "def New_Rap(k_n,x_tilda,x_bar): #Newthon-Raphson Method will return a value closer to the root for the given initial value.\n",
    "\n",
    "\tnew_k_n = k_n - (kappa_function(k_n,x_tilda,x_bar)/fd_kappa_function(k_n))\n",
    "\t#print(new_k_n)\n",
    "\n",
    "\treturn new_k_n\n",
    "\n",
    "k_n = 1 # take initial kappa value 1 (k>0) \n",
    "#Initial value can be method of moment estimate:np.mean(data)**2/np.var(data) \n",
    "##https://ocw.mit.edu/courses/mathematics/18-443-statistics-for-applications-spring-2015/lecture-notes/MIT18_443S15_LEC3.pdf\n",
    "tolerance = 10**-5\n",
    "\n",
    "print(\"\\nInitial Kappa = {}\".format(k_n))\n",
    "print(\"Function Value in Initial Kappa = {}\\n\".format(kappa_function(k_n,x_tilda,x_bar)))\n",
    "\n",
    "while(abs(kappa_function(k_n,x_tilda,x_bar))>tolerance):\n",
    "\n",
    "\tk_n = New_Rap(k_n,x_tilda,x_bar)\n",
    "\tprint(\"New Kappa = {}\".format(k_n))\n",
    "\tprint(\"Function Value in New Kappa = {}\\n\".format(kappa_function(k_n,x_tilda,x_bar)))\n",
    "\n",
    "#Find theta_hat\n",
    "theta_hat = np.mean(data)/k_n\n",
    "print(\"Theta Hat = {}\".format(theta_hat))\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "sort_data = np.sort(data)\n",
    "\n",
    "plt.plot(sort_data, stats.gamma.pdf(sort_data, a = 2, scale = 2))\n",
    "\n",
    "plt.show() \n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46b8c15b",
   "metadata": {},
   "source": [
    "🏆: ☕ + 🍰 at Espresso Lab, MED."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f6dced3-bfa3-43e3-ac43-96a957f085de",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}