dan_cos.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "from scipy.stats import rv_continuous\n",
    "from scipy import stats\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using scipy's slightly horrible but largely wonderful rv_continuous\n",
    "The syntax is whacky, but it's otherwise a miracle..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class _cos_pdf(stats.rv_continuous):\n",
    "\n",
    "    def _pdf(self, X): return np.cos(X) / 2\n",
    "\n",
    "# This is for the \"base\" model, which can be scaled and offset \n",
    "#  by loc and scale.  The bounds refer to scale = 1 and loc = 0.\n",
    "# This creates an instance... it's very convoluted.\n",
    "cos_rv = _cos_pdf(a = -np.pi/2, b = np.pi/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiUAAACtCAYAAABvEQ93AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEaJJREFUeJzt3XHMXXV5wPHvIzWWUN5lSCXR5bXBzHQrWhfeRJdt0YgL6NxG7P5gzIUmQomGJRsk2j8oqyJBcP1L2aABxxBdGEnBSckSyDQL2aJpYgo26UhA69BhWmJeaVWc+OyPe265Xt63vbfvPff8zjnfT3LSe8/v3Pf93dPnnPe5z+937onMRJIkqWmvaboDkiRJYFIiSZIKYVIiSZKKYFIiSZKKYFIiSZKKYFIiSZKKYFIiSZKKYFIiSZKKYFIiSZKKUGRSctlllyXg4jKrZe6MYZcZL40wjl1muEykyKTk2LFjTXdBWhNjWF1gHGveikxKJElS/5iUSJKkIpiUSJKkIpiUSJKkIqxrugN9tmnn/pOPv/uZP2qwJ1L9hvFurKsPjPczY6VEkiQVwaREkiQVweEbSZIm4JB7/ayUqBci4v6I+N+I+HFEPB0RV4+0XRIRhyPiJxHxtYh480hbRMRtEfFCtdweEdHMu+iPTTv3/8ofAEn9YFKivrgV2JSZC8CfAJ+OiIsj4nxgH7ALOA84ADww8rodwOXAVuDtwAeBa+fZ8VKYKKhPjPdmOHzTMpYPz0xmHhp9Wi1vAS4GDmXmgwARsRs4FhGbM/MwcBWwJzOfq9r3ANcAd86x+60yTYzWta2kdjIpUW9ExN8D24GzgW8BjwK3AAeH22TmiYh4BtgCHK7+PTjyYw5W61b6+TsYVFZYXFyc/Rso3Fo/VU7zei+3VBuZWJ+ewzfqjcz8GHAu8AcMhmxeAjYAy2ObLlfbsUL7MrBhpXklmbk3M5cyc2njxo2z7r4kdZ6VEvVKZr4MPBERHwY+ChwHFsY2WwBerB6Pty8AxzNz4ltx61fNcpzeT55qykrVOuegrJ1JSYd5wj6ldQzmlBxiMG8EgIg4Z2Q91b9bgW9Wz7eOtKkgfYj3iLgfuAQ4B3geuD0z767aLgHuABaBbwDbM/NI1RbAZ4DhVWf3AJ8wuf5VJhXNMykpWB9OsvMQEW8A3gs8AvwUeB/w58CVwH8Cn42IbcB+4CbgyWqSK8B9wPUR8SiDybE3AJ+b7zvoL/9IvMqtwEcy86WI2Ax8PSK+BRxhMCR5NfBV4GYGV5G9q3rd6FVkCTwGPIsTto2xwpiUtIQHzpokg6GaOxnMozoC/HVmfgWgSkg+D9zP4BPmFSOvvQu4EHiqen53tU4TMG5ny6vI1HUmJeq8zDwKvPsU7Y8Dm1dpS+Dj1SI1zqvI1GUmJYXxk6WkU8nMj0XEXwG/C7yHV64iOzq26URXkY3PK8nMvcBegKWlJeeczIGXuL/CpETSGSs1ie76Sd6ryNRVEyUlzvgukxNhpd7zKrKGlZqYt9WkX57mfUMK530aNCvDWDKeyhIRb4iIKyJiQ0ScFRGXMriK7N+Bh4CLImJbRKxn9avI3hQRb2RwFdm9DbyNuTGG22miSokzviUNeaJvjFeRrZHV5fJNPKfEGd+S1ByvIlMfTJyUOOO7HfwUK0lqq6muvnHGd/tZvtQ0THLVVcZ2mc70kmBnfEuSNEN+aJwgKfG+IZKktlprRaTp1/fNJJUSZ3x3gAeGJJXDc/LKTpuUOOO7uywVSpJK4tfMS+osE291VVdj26SkEJbyJEl9Z1IyZ6UnH13NviWpTfp6Lp703jeSJEm1slJSo9KrIqPa1FdJUjeZlEg91tcSsdQmw+O0D8eowzeSJKkIVkok9YJVIbVdH4bZTUpmrA9BI0lSHRy+kSRJRTApkSRJRXD4RlJvOc9EKotJidQTp7us0PlQkppmUiJJWDWRSmBSIknqBKt97WdSolX5yVF91adv0FR3dOGcbVKyBp64pHbyE7W6pEvx7CXBkiSpCFZKNJGVqkKrZedWjiRJZ8JKiTovIl4XEfdExJGIeDEivhUR7x9pvyQiDkfETyLiaxHx5pG2iIjbIuKFark9IqKZdyJp3Kad+08uaj8rJTPgwVC8dcD/AO8Gvgd8APiXiHgbcBzYB1wNfBW4GXgAeFf12h3A5cBWIIHHgGeBO+fY/5kyXtUFxnE3mZSo8zLzBLB7ZNUjEfEd4GLg9cChzHwQICJ2A8ciYnNmHgauAvZk5nNV+x7gGlqclEhSqUxKNJUufDqJiAuAtwKHgI8CB4dtmXkiIp4BtgCHq38Pjrz8YLVupZ+7g0FlhcXFxVr6Lkld5pwS9UpEvBb4EvBPVSVkA7A8ttkycG71eLx9Gdiw0rySzNybmUuZubRx48bZd1695two9YFJiXojIl4DfBH4OXBdtfo4sDC26QLw4irtC8DxzMwauyqtZHRu1K8BuxjMjdoUEeczmBu1CzgPOMBgbtTQ6NyotwMfBK6dX9elyZiUqBeqT4X3ABcA2zLz/6qmQwxO1MPtzgHeUq1/VXv1+BDSnGXmiczcnZnfzcxfZuYjwHBu1Ieo5kZl5s8YzKHaGhGbq5efnBuVmd8H9gDb5/8upFNzTsmUujCnoqf+Afgt4H2Z+dOR9Q8Bn42IbcB+4CbgyWpoB+A+4PqIeJTB1Tc3AJ+bX7ellfVlblQXvjpdkzttpcRxTLVdFZPXAu8Ano+I49XyF5l5FNgG3AL8CHgncMXIy+9icKnwU8C3GSQud82z/9I450apqyaplPT+Ox6sjrRbZh4BVk2GM/NxYPMqbQl8vFqkxjk3Sl122qTE73iQ2suEulvG5kZ9YGxu1FUj2602N+qb1fOi50atFrfGc/dNPaekL+OYUlt54u4050ap06ZKSsbHMSNiA3B0bLOJxjHHy4aZuRfYC7C0tGRJUVLjSppkOTI36iUGc6OGTddm5peqhOTzwP3AN3j13KgLGcyNArgb50Z1WkmxO42JkxLHMSWpOc6NUh9M9D0lfseDJEmq26RfnjYcx/zjFcYxL4qIbRGxntXHMd8UEW9kMI5572y6LkmSuuS0wzd9Hsd0wqAkSfMzySXBjmNKkqTaee8bzdymnfutMkmSpua9b1Sbtl6SJq3EeFbbtSGGTUqkjrA6JWklbTo3OHwjSZKKYKVkTJsySkmSusRKiSRJKoKVEklSEaxUq9eVEi9dlSSpHFZKMDuXJKkEJiWSNKU2fN+D1Ea9Hr6RJEnlMCmRJElFcPhGkqSeKXUI0qREczc8GEo6ENrKSdqSusThG0mSVAQrJZoLP9FLWonnBo2yUiJJa+CXMEqz07tKiScPSZLKZKVEkiQVwaREkiQVoXfDN5KkZjmMXpaSvqahF0mJB4DartQvOpKkWepFUiJ1iUm2pK5yTol6ISKui4gDEfFSRNw71nZJRByOiJ9ExNci4s0jbRERt0XEC9Vye0TE3N+AJPWAlRL1xQ+ATwOXAmcPV0bE+cA+4Grgq8DNwAPAu6pNdgCXA1uBBB4DngXunFfH1Q4OsantSohhkxI1ZrVhiDoOhszcBxARS8BvjDR9CDiUmQ9W7buBYxGxOTMPA1cBezLzuap9D3ANJiWas4i4DtgOvA3458zcPtJ2CXAHsAh8A9iemUeqtgA+wyDxBrgH+ERm5tw6X3HoUadjUqK+2wIcHD7JzBMR8Uy1/vB4e/V4y0o/KCJ2MKissLi4WFd/1QI1feK02qfOm2hOiePx6rANwPLYumXg3FXal4ENK8VxZu7NzKXMXNq4ceNMOudXmGsoM/dl5sPAC2NNJ6t9mfkzYDewNSI2V+0nq32Z+X1gD4OKi1ScSSe6DjP0L4yuHMnQdwHnAQcYZOhDoxn624EPAteurcvqujn/IT4OLIytWwBeXKV9ATjeROlbWsWrqn3AsNr3qnZOUe2DQcWv+hB64OjRozV0V1rdRMM3jserww4xiFMAIuIc4C3V+mH7VuCb1fOtI21SCTYA49nDRNW+lZLrzNwL7AVYWlpac/JtpU/TWOslwTPL0M3OVaeIWBcR64GzgLMiYn1ErAMeAi6KiG1V+03Ak1VSDXAfcH1EvCki3gjcANzbwFuQVmO1T52x1qSk6PF4acSNwE+BncCHq8c3ZuZRYBtwC/Aj4J3AFSOvu4vB5MGngG8D+6t1UimG1TzglNW+Iat9KtZar74xQ1crZOZuBhMAV2p7HNi8SlsCH68WqTFVZW8dI9U+4BcMqn2fjYhtDJLm1ap9jzK4+uYG4HPz7r80ibUmJY7HS9J83Aj87cjzDwOfzMzdVULyeeB+Bt9TMl7tu5BBtQ/gbqz26TSa+iK1iZKStmboTrCS1BVW+9QHk84pcTxekiTVatJLgndjhi5JkmrkXYIlSVIRvPeNVBjnQknqK5MSSdJMmVjrTDl8I0k18qaK0uQ6Vynx4JdUoqa+90Fqk84lJeoOT+KS1Lx5nosdvpEkSUUwKZEkSUVw+EaSNBPO6dNamZRIBfBkLkkmJWoJJ71KUvc5p0SSJBXBSokkzZmVP2llVkokSVIROlEpcZKgJDXD869myUqJJEmaSN33cmptpcTsXJKkZtQ1L8pKiSRJKkJrKyVSF1jxk6RXWCmRJElFMCmRJElFMClR69Q9+1uS1AznlEhzZkIlqUtmeSWOlRJJapCVP+kVrauUePBKktRNrUtKpCFvaiY1ww+HqovDN5IkqQi1JyURcV5EPBQRJyLiSERcWffvlGbJGFYXGMdqg3kM39wB/By4AHgHsD8iDmbmoTn8bmkWjGF1gXGs4tVaKYmIc4BtwK7MPJ6ZTwD/Cvxlnb9XmhVjWF1gHKst6q6UvBV4OTOfHll3EHj3+IYRsQPYUT09HhH/XXPfVnI+cKyB33um2tTfWvsat52y+d8y87Iz/NFti2EwLupUW39rjGGoN47b9n84a31//zCyD04RxxPFcN1JyQZgeWzdMnDu+IaZuRfYW3N/TikiDmTmUpN9mEab+tumvo5pVQxDu/Z1m/oK7evviNriuMX7ZCb6/v5htvug7omux4GFsXULwIs1/15pVoxhdYFxrFaoOyl5GlgXEb85sm4r4MQqtYUxrC4wjtUKtSYlmXkC2Ad8KiLOiYjfA/4U+GKdv3cNGi+9T6lN/W1TX09qYQxDu/Z1m/oK7esvUHsct3KfzFDf3z/McB9EZs7qZ638CyLOA74A/CHwArAzM79c6y+VZsgYVhcYx2qD2pMSSZKkSfg185IkqQgmJZIkqQgmJSuIiOsi4kBEvBQR9zbdn3FtuodF6fuyC6aJh4jYHhEvR8TxkeU9BfXvbyLi+YhYjogvRMTr6uzbmfa1if1YgtJjrS5tiuE6zPO4mMe9b9roB8CngUuBsxvuy0radA+L0vdlF0wbD/+Vmb8/t95N2L+IuBTYCbyXQdw8BHyyWldUXyvz3o8lKD3W6tKmGK7D3I4LKyUryMx9mfkwgxnqRWnbPSxK3pddUHo8TNm/q4B7MvNQZv4IuBnYXmhfe6ev+6dNMVyHef+/m5S0z2r3sNjSUH/UrDOJh9+JiGMR8XRE7IqIOium0/RvS9U2ut0FEfH6Gvs3atp9Oc/9WILSY60ubYrhOsz1uGhjgPTdxPewUC9MGw//AVwEHGFwUnkA+AVwawH9G992+Phc5lNpm6av896PJSg91urSphiuw1yPi95VSiLi6xGRqyxPNN2/CXgPix6ZIF6niofMfDYzv5OZv8zMp4BPAX9W41uYpn/j2w4fzyu2J+5rA/uxdh2Itbq0KYbrMNfjondJSWa+JzNjlaUNE7K8h0WPTBCva42HBGLW/R4xTf8OVW2j2/0wM+f1CXMt+7Lu/Vi7DsRaXdoUw3WY63HRu6RkEhGxLiLWA2cBZ0XE+lLGQtt2L5aS92UXTBsPEfH+iLigerwZ2AV8pZD+3Qd8JCJ+OyJ+HbgRuLeuvq2lr/PejyUoPdbq0qYYrsPcj4vMdBlbgN0MMrzRZXfT/Rrp33nAw8AJ4HvAlU33qa37sgvLqeIBWGRQfl2snv8d8MNq22cZlFdf20T/xvtWrbu+6t+PgX8EXlfCvixhP5awlB5rTcdFta7RGG7y/c/i/9x730iSpCI4fCNJkopgUiJJkopgUiJJkopgUiJJkopgUiJJkopgUiJJkopgUiJJkopgUiJJkorw/7Y1VWCfNx1PAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x180 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 3, figsize = (9, 2.5))\n",
    "\n",
    "for xi, axi in enumerate(ax):\n",
    "    axi.hist(cos_rv(loc = 0, scale = 1.0 / (xi+1)).rvs(10000), bins = 50);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = cos_rv(loc = 0, scale = 1.2345).rvs(10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.fit.html#scipy.stats.rv_continuous.fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0010984080928232165, 1.2311065892727473)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cos_rv.fit(data) # MLE fit!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To do it by hand..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trunc_cos(p, X):\n",
    "    \n",
    "    return np.where(np.abs(X) < np.pi / 2 / p[0], p[0] * np.cos(p[0] * X) / 2, 1e-12)\n",
    "\n",
    "def cos_nll(p, data): \n",
    "    \n",
    "    return - np.sum(np.log(trunc_cos(p, data)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.5, 2.488769531250003)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "θ = 2.5\n",
    "\n",
    "cos_inst = cos_rv(loc = 0, scale = 1 / θ)\n",
    "data = cos_inst.rvs(10000)\n",
    "\n",
    "result = minimize(cos_nll, x0 = [1], args = (data, ), \n",
    "                  method = 'Nelder-Mead', options={'maxiter' : 5e3, 'disp': False})\n",
    "\n",
    "θ, result.x[0]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}