MoM Gamma python3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e0ef60be",
   "metadata": {},
   "source": [
    "# Standard Gamma model: introduction\n",
    "\n",
    "We consider the following Gamma model, with the two parameters $\\alpha$ and $\\beta$ unknown\n",
    "\n",
    "\n",
    "$$\n",
    "\\bigg\\{  Ga\\big( \\alpha, \\beta   \\big) ; \\alpha > 0 , \\ \\beta > 0 \\bigg\\} \n",
    "$$\n",
    "\n",
    "The Probability Density Function (PDF) of the Gamma distribution is given by\n",
    "\n",
    "$$\n",
    "f_{\\alpha, \\beta}(x) = \\frac{  \\beta^{\\alpha}   }{ \\Gamma(\\alpha)  } x^{\\alpha - 1}   e^{- \\beta x} \n",
    "$$\n",
    "\n",
    "\n",
    "The gamma function $\\Gamma(\\cdot)$, also called the Euler integral of the second kind is defined, for a nonnegative integer $\\alpha$, as $\\Gamma(\\alpha) = \\int\\limits_{0}^{\\infty} \\ x^{\\alpha -1}   e^{-x} dx.\n",
    "\n",
    "# Method of Moments (MoM) estimation\n",
    "\n",
    "\n",
    "\\textbf{Method of Moments (MoM)}: If $ E[\\varphi(X) ] = h(\\theta)$, then we have that $\\hat{\\theta}_{MOM} = h^{-1}\\bigg(  \\frac{1}{n} \\sum\\limits_{i=1}^{n} \\varphi(x_{i})  \\bigg)$ is a method of moments estimator for $\\theta$.\n",
    "\n",
    "To obtain a \\textbf{Method of Moments estimator} for the parameter vector $\\boldsymbol{\\theta} = ( \\alpha \\ , \\  \\beta)$ from a sample $x_{1},..., x_{n}$ of $n$ $i.i.d.$ realizations of a Gamma  distributed r.v $X$, we need the moments of $X$, given by \n",
    "\n",
    "$$ E[X] = \\frac{\\alpha}{\\beta}  \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\  E[X^2] = \\frac{\\alpha(\\alpha+1)}{\\beta^{2}}$$\n",
    "\n",
    "\n",
    "# Derivation of the first moment\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{align*} \n",
    "& E[X] = \\int\\limits_{-\\infty}^{+\\infty} x \\  f_{\\alpha, \\beta}(x)     \\ dx =   \\int\\limits_{0}^{\\infty} x \\ \\frac{  \\beta^{\\alpha}   }{ \\Gamma(\\alpha)  } x^{\\alpha - 1}   e^{- \\beta x}      \\ dx  \\\\\n",
    "& \\ \\ \\ \\ \\ \\ \\  = \\frac{  \\beta^{\\alpha}   }{ \\Gamma(\\alpha)  }           \\int\\limits_{0}^{\\infty} x \\  x^{\\alpha - 1}   e^{- \\beta x}      \\ dx    \\ \\  \\ (u= \\beta x  \t\\Leftrightarrow  x= u/ \\beta, \\ \\ \\ du = \\beta dx  \\Leftrightarrow  dx= du / \\beta ) \\\\\n",
    "& \\ \\ \\ \\ \\ \\ \\  = \\frac{  \\beta^{\\alpha}   }{ \\Gamma(\\alpha)  }           \\int\\limits_{0}^{\\infty}      \\bigg(   \\frac{u}{\\beta}   \\bigg)^{\\alpha}   e^{- u}      \\ \\frac{du}{\\beta}    \\\\          \n",
    "& \\ \\ \\ \\ \\ \\ \\  = \\frac{  \\beta^{\\alpha}   }{ \\Gamma(\\alpha)  }         \\frac{1}{\\beta^{\\alpha}}         \\frac{1}{\\beta}         \\int\\limits_{0}^{\\infty}     u^{\\alpha}   e^{- u}      \\  du      \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\     \\int\\limits_{0}^{\\infty}     x^{\\alpha}   e^{- x}      \\  dx  =  \\Gamma(\\alpha +1)  = \\alpha  \\Gamma(\\alpha)    \\\\ \n",
    "& \\ \\ \\ \\  \\ \\ \\  = \\frac{  \\beta^{\\alpha}   }{ \\Gamma(\\alpha)  }         \\frac{1}{\\beta^{\\alpha +1}}     \\alpha  \\Gamma(\\alpha)   \\\\\n",
    "& \\ \\ \\ \\  \\ \\ \\  = \\frac{\\alpha}{\\beta}  \\ \\ \\ \\ \\ \\ \\ \\  \\  \\\\  \n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "# Derivation of the second moment\n",
    "\n",
    "$$\n",
    "\\begin{align*} \n",
    "& E[X^2] = \\int\\limits_{-\\infty}^{+\\infty} x^2 \\  f_{\\alpha, \\beta}(x)     \\ dx =   \\int\\limits_{0}^{\\infty} x^2 \\ \\frac{  \\beta^{\\alpha}   }{ \\Gamma(\\alpha)  } x^{\\alpha - 1}   e^{- \\beta x}      \\ dx  \\\\\n",
    "& \\ \\ \\ \\ \\ \\ \\  = \\frac{  \\beta^{\\alpha}   }{ \\Gamma(\\alpha)  }           \\int\\limits_{0}^{\\infty} x^2 \\  x^{\\alpha - 1}   e^{- \\beta x}      \\ dx    \\ \\ \\ \\ \\ \\ \\ \\ (u= \\beta x  \t\\Leftrightarrow  x= u/ \\beta, \\ \\ \\ du = \\beta dx  \\Leftrightarrow  dx= du / \\beta ) \\\\\n",
    "& \\ \\ \\ \\ \\ \\ \\  = \\frac{  \\beta^{\\alpha}   }{ \\Gamma(\\alpha)  }           \\int\\limits_{0}^{\\infty}      \\bigg(   \\frac{u}{\\beta}   \\bigg)^{\\alpha + 1}   e^{- u}      \\ \\frac{du}{\\beta}    \\\\       \n",
    "& \\ \\ \\ \\ \\ \\ \\  = \\frac{  \\beta^{\\alpha}   }{ \\Gamma(\\alpha)  }         \\frac{1}{\\beta^{\\alpha + 1}}         \\frac{1}{\\beta}         \\int\\limits_{0}^{\\infty}     u^{\\alpha}   e^{- u}      \\  du      \n",
    "\\ \\ \\     \\int\\limits_{0}^{\\infty}     x^{\\alpha + 1}   e^{- x}      \\  dx  =  \\Gamma(\\alpha + 2)   = (\\alpha+1) \\alpha \\Gamma(\\alpha)   \\\\ \n",
    "& \\ \\ \\ \\  \\ \\ \\  = \\frac{  \\beta^{\\alpha}   }{ \\Gamma(\\alpha)  }         \\frac{1}{\\beta^{\\alpha +2}}     (\\alpha+1) \\alpha \\Gamma(\\alpha) \\\\\n",
    " & \\ \\ \\ \\  \\ \\ \\ = \\frac{(\\alpha + 1) \\alpha}{\\beta^{2}}  \\ \\ \\ \\ \\ \\  \\\\  \n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "# Gamma MoM estimators\n",
    "\n",
    "\n",
    "So we set $\\overline{x} = u$ and $\\overline{x^{2}} = v$ and solve for $\\alpha$ and $\\beta$ in the following equation:\n",
    "\n",
    "$$   h \\bigg( \\begin{matrix} \\alpha\\\\ \\beta \\end{matrix}  \\bigg) = \\Bigg( \\ \\begin{matrix}   \\alpha / \\beta \\\\   \\alpha(\\alpha+1) / \\beta^{2}  \\end{matrix}  \\Bigg)  = \\Bigg( \\ \\begin{matrix}   u \\\\   v  \\end{matrix}  \\Bigg) \n",
    "$$\n",
    "\n",
    "$$\n",
    " \\begin{array}{lll}\n",
    "  u = \\alpha / \\beta    \\\\\n",
    "  v = \\alpha(\\alpha+1) / \\beta^{2}  \\ \\ \\  \\ \\ \\Leftrightarrow \\ \\ \\ \\ \\ \\ \\ \\ v  =  u^{2}+ \\frac{u}{\\beta}   \\ \\ \\  \\ \\ \\Leftrightarrow \\ \\ \\ \\ \\ \\ \\ \\ v -  u^{2}= \\frac{u}{\\beta}   \\ \\ \\  \\ \\ \\Leftrightarrow \\ \\ \\ \\ \\ \\ \\ \\ \\beta  = \\frac{u}{v-u^{2}}\\\\\n",
    " \\end{array}\n",
    "$$\n",
    "\n",
    "\n",
    "Then it is easy to show that  $\\alpha = u \\beta  \\ \\ \\  \\ \\ \\Leftrightarrow \\ \\ \\ \\ \\ \\ \\ \\ \\alpha = \\frac{u^{2}}{v-u^{2}}$\n",
    "\n",
    "So we get as our Method of Moments estimators:\n",
    "\n",
    "$$\n",
    "\\bigg( \\frac{\\overline{x}^{2} }{\\overline{x^{2}} - \\overline{x}^{2}},    \\frac{\\overline{x} }{\\overline{x^{2}} - \\overline{x}^{2}}     \\bigg) =  \\bigg(  \\hat{\\alpha}_{MOM}, \\hat{\\beta}_{MOM}   \\bigg)\n",
    "$$\n",
    "\n",
    "# Example 1 - Lifetime of a product\n",
    "\n",
    "\n",
    "The Gamma distribution is often used in reliability analysis to model time until failure kind of events, in particular Poisson distributed events. Suppose that we test the lifetime of a specific product on a sample of size $100$ and obtain the following data, in weeks (20 first entries of the dataset)\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "x_{i} & = 17.49, 14.93, 16.01, 17.85,  7.37,  5.96, 12.27, 16.44, 28.34, 30.28, \\\\\n",
    "& \\ \\ \\ \\ 10.66, 25.77, 17.21, 34.99, 28.53, 12.98,  9.32, 22.67, 12.55, 16.81\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "\n",
    "(i) \\ \\ From which Gamma distribution have the data been generated if we use the Method of Moment estimators for the parameters?\n",
    "\n",
    "\n",
    "(ii) \\ \\ What are the cumulative hazard rates after 12, 24 and 36 weeks (or the total number of times we expect a failure for one product)?\n",
    "\n",
    "# Example 1 - Estimation\n",
    "\n",
    "\n",
    "The full dataset is available here:\n",
    "\n",
    "\\urlf{https://github.com/JRigh/Gamma-distribution-estimation/blob/ \\\\\n",
    "main/dataset.csv}\n",
    "\n",
    "\n",
    "The data have been generated from a $Ga\\big( \\alpha = 2, \\beta = 0.08   \\big)$ and computing the MoM estimators in R and Python, we can obtain a $Ga\\big( \\alpha = 1.900579, \\beta = 0.07399452    \\big)$ with R and a $Ga\\big( \\alpha = 1.900579, \\beta = 0.073995   \\big)$ with Python, so equivalent results up to rounding.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "d6af0204",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'alpha_mom': x    1.900579\n",
      "dtype: float64, 'beta_mom': x    0.073995\n",
      "dtype: float64}\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# function that implements the estimators above\n",
    "def MoMgamma(x):\n",
    "    n = len(x)\n",
    "    sample_moment_1 = np.sum(x) / n\n",
    "    sample_moment_2 = np.sum(x**2) / n\n",
    "\n",
    "    alpha_mom = sample_moment_1**2 / (sample_moment_2 - sample_moment_1**2)\n",
    "    beta_mom = sample_moment_1 / (sample_moment_2 - sample_moment_1**2)\n",
    "\n",
    "    output = {\n",
    "        'alpha_mom': alpha_mom,\n",
    "        'beta_mom': beta_mom\n",
    "    }\n",
    "\n",
    "    return output\n",
    "\n",
    "# read the dataset\n",
    "x = pd.read_csv(\"C:/Users/julia/OneDrive/Desktop/github/37. momgamma/dataset.csv\")\n",
    "\n",
    "# perform MoM estimation\n",
    "result = MoMgamma(x)\n",
    "print(result)"
   ]
  }
 ],
 "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}