{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 273,
   "id": "98d26b7b-0ddb-4ef8-93af-106818267229",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import scipy.special as sp\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.optimize as opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bf2aad2a-b9ad-4d38-9763-4182752b5d54",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def f(L, z):\n",
    "    x, y = L[0], L[1]\n",
    "    return np.power(z*x, 2) + np.power(y/z, 2) - np.cos(2 * np.pi * x) - np.cos(2 * np.pi * y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bc00dbe1-2edf-4a90-839e-1aa8d53fc7dc",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  message: Optimization terminated successfully.\n",
      "  success: True\n",
      "   status: 0\n",
      "      fun: -1.0125323577209524\n",
      "        x: [ 1.409e-07  1.975e+00]\n",
      "      nit: 6\n",
      "      jac: [ 7.048e-06  3.651e-06]\n",
      " hess_inv: [[ 2.137e-02  2.198e-04]\n",
      "            [ 2.198e-04  2.546e-02]]\n",
      "     nfev: 30\n",
      "     njev: 10\n"
     ]
    }
   ],
   "source": [
    "res2 = opt.minimize(f, x0=(1, 2), args=2, method=\"BFGS\")\n",
    "print(res2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "9a50016b-d8c1-489d-affe-81bae33463dc",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                    message: ['requested number of basinhopping iterations completed successfully']\n",
      "                    success: True\n",
      "                        fun: -2.0\n",
      "                          x: [ 1.358e-09  6.023e-10]\n",
      "                        nit: 100\n",
      "      minimization_failures: 0\n",
      "                       nfev: 3369\n",
      "                       njev: 1123\n",
      " lowest_optimization_result:  message: Optimization terminated successfully.\n",
      "                              success: True\n",
      "                               status: 0\n",
      "                                  fun: -2.0\n",
      "                                    x: [ 1.358e-09  6.023e-10]\n",
      "                                  nit: 7\n",
      "                                  jac: [ 4.172e-07  3.278e-07]\n",
      "                             hess_inv: [[ 2.140e-02  3.244e-04]\n",
      "                                        [ 3.244e-04  2.532e-02]]\n",
      "                                 nfev: 30\n",
      "                                 njev: 10\n"
     ]
    }
   ],
   "source": [
    "mk = {\"method\": \"BFGS\", \"args\": 2}\n",
    "res = opt.basinhopping(f, x0=(1, 2), minimizer_kwargs=mk)\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "id": "4bbd6193-9658-421a-a746-81b9fd530f2c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 80956 characters]",
      "text/plain": [
       "<Figure size 1800x500 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(18, 5))\n",
    "for i, nb_repl in enumerate([100, 1000, 10000]):\n",
    "    plt.subplot(1, 3, i+1)\n",
    "    sample_X1 = np.random.normal(0, 1, nb_repl)\n",
    "    sample_X2 = np.random.normal(3, np.sqrt(5), nb_repl)\n",
    "    sample_e = np.random.normal(0, np.sqrt(1/4), nb_repl)\n",
    "    Y = 5 * sample_X1 - 4 * sample_X2 + 2 + sample_e\n",
    "\n",
    "    intervalle = np.linspace(np.min(Y), np.max(Y), 100)\n",
    "    plt.hist(Y, bins=intervalle, density=True, label='Echantillon de Y')\n",
    "\n",
    "    densite = stats.norm(-10, np.sqrt(105.25)).pdf\n",
    "    plt.plot(intervalle, densite(intervalle), label='Fonction densité')\n",
    "\n",
    "    plt.title(f\"Graphique de la somme de gaussiennes pour N={nb_repl}\")\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "id": "508aa866-d4ed-42a6-8f01-09c9b8721566",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def theta_hat(Y):\n",
    "    return np.mean(Y)\n",
    "\n",
    "def sigma_hat(Y):\n",
    "    return 1/nb_repl * np.sum(np.power(Y - theta_hat(Y), 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "id": "eb62aac4-2a4f-4ae8-84e9-702bc212315e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def log_likehood_gauss(X, Y):\n",
    "    theta, sigma_2 = X[0], X[1]\n",
    "    return 1/2*np.log(2*np.pi) + 1/2*np.log(sigma_2) + 1/(2*nb_repl*sigma_2) * np.sum(np.power(Y - theta, 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "id": "ae68062d-1c5d-4881-b750-afaefd62b8af",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-10.15254483 105.22161157]\n",
      "[-10.152550623018831, 105.22156995434973]\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "nb_repl = 5000\n",
    "sample_X1 = np.random.normal(0, 1, nb_repl)\n",
    "sample_X2 = np.random.normal(3, np.sqrt(5), nb_repl)\n",
    "sample_e = np.random.normal(0, np.sqrt(1/4), nb_repl)\n",
    "Y = 5 * sample_X1 - 4 * sample_X2 + 2 + sample_e\n",
    "    \n",
    "mk = {\"method\": \"BFGS\", \"args\": Y}\n",
    "res = opt.basinhopping(log_likehood_gauss, x0=(-1, 98.75), minimizer_kwargs=mk)\n",
    "print(res.x)\n",
    "print([theta_hat(Y), sigma_hat(Y)])\n",
    "print(np.allclose([theta_hat(Y), sigma_hat(Y)], res.x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "id": "2ce462b1-00e8-4386-bb8c-163f2ba69e68",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def simule(a, b, n):\n",
    "    X = np.random.gamma(a, 1/b, n)\n",
    "    intervalle = np.linspace(0, np.max(X), 100)\n",
    "    plt.hist(X, bins=intervalle, density=True, label='Echantillon de X')\n",
    "\n",
    "    densite = stats.gamma.pdf(intervalle, a, 0, 1/b)\n",
    "    plt.plot(intervalle, densite, label='Fonction densité Gamma(2, 1)')\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 259,
   "id": "b8861952-f676-4b38-8f24-ecd71b137bfd",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 34308 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "simule(2, 1, 10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 290,
   "id": "dcb07fde-0af2-46e0-b009-bf4d799249c1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def log_likehood_gamma(X, sample):\n",
    "    a, b = X[0], X[1]\n",
    "    n = len(sample) \n",
    "    return -n*a*np.log(b) + n * np.log(sp.gamma(a)) - (a-1) * np.sum(np.log(sample)) + b * np.sum(sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 296,
   "id": "ea041002-cc5c-47e1-aa52-9a607def05d2",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.99590773 0.96694005]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_360/437431687.py:4: RuntimeWarning: invalid value encountered in log\n",
      "  return -n*a*np.log(b) + n * np.log(sp.gamma(a)) - (a-1) * np.sum(np.log(sample)) + b * np.sum(sample)\n",
      "/tmp/ipykernel_360/437431687.py:4: RuntimeWarning: divide by zero encountered in log\n",
      "  return -n*a*np.log(b) + n * np.log(sp.gamma(a)) - (a-1) * np.sum(np.log(sample)) + b * np.sum(sample)\n",
      "/tmp/ipykernel_360/437431687.py:4: RuntimeWarning: invalid value encountered in log\n",
      "  return -n*a*np.log(b) + n * np.log(sp.gamma(a)) - (a-1) * np.sum(np.log(sample)) + b * np.sum(sample)\n",
      "/tmp/ipykernel_360/437431687.py:4: RuntimeWarning: divide by zero encountered in log\n",
      "  return -n*a*np.log(b) + n * np.log(sp.gamma(a)) - (a-1) * np.sum(np.log(sample)) + b * np.sum(sample)\n",
      "/tmp/ipykernel_360/437431687.py:4: RuntimeWarning: invalid value encountered in log\n",
      "  return -n*a*np.log(b) + n * np.log(sp.gamma(a)) - (a-1) * np.sum(np.log(sample)) + b * np.sum(sample)\n",
      "/tmp/ipykernel_360/437431687.py:4: RuntimeWarning: divide by zero encountered in log\n",
      "  return -n*a*np.log(b) + n * np.log(sp.gamma(a)) - (a-1) * np.sum(np.log(sample)) + b * np.sum(sample)\n"
     ]
    }
   ],
   "source": [
    "nb_repl = 1000\n",
    "a, b = 2, 1\n",
    "\n",
    "sample = np.random.gamma(a, 1/b, nb_repl)\n",
    "mk = {\"method\": \"BFGS\", \"args\": sample}\n",
    "res = opt.basinhopping(log_likehood_gamma, x0=(1, 1), minimizer_kwargs=mk)\n",
    "print(res.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abfb9307-2531-4e48-818d-169ebe17072b",
   "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
