{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytic maximum likelihood scaled transit fit\n",
    "\n",
    "We can use a factor $s$ to modifiy the transit depth in a nominal model ${\\bf m_0}$ calculated with approximately the correct depth that is scaled as follows ...\n",
    "\n",
    "$ {\\bf m}(s) = 1 + s\\times({\\bf m_0}- 1)$\n",
    "\n",
    "The data are normalised fluxes ${\\bf f} = f_1, \\dots, f_N$ with nominal errors ${\\bf \\sigma} = \\sigma_1,\\dots, \\sigma_N$. \n",
    "\n",
    "Assume that the actual standard errors are underestimated by some factor $b$, and that these are normally distributed and uncorrelated, so that the log-likelihood is\n",
    "\n",
    "$\\ln {\\cal L}  = -\\frac{1}{2b^2}\\chi^2 - \\frac{1}{2}\\sum_{i=1}^N \\ln \\sigma_i^2  - N\\ln b - \\frac{N}{2}\\ln(2\\pi)$\n",
    "\n",
    "where\n",
    "\n",
    "$\\chi^2 = \\sum_i^N \\frac{(f_i - 1 - s(m_{0,i}-1)^2}{\\sigma_i^2} $\n",
    "\n",
    "\n",
    "The maximum likelihood occurs for parameter values $\\hat{s}$, and $\\hat{b}$ such that\n",
    "\n",
    "$\\left. \\frac{\\partial  \\ln {\\cal L}}{\\partial s}\\right|_{\\hat{s},\\hat{b}} = 0$\n",
    "and \n",
    "$\\left. \\frac{\\partial  \\ln {\\cal L}}{\\partial b}\\right|_{\\hat{s},\\hat{b}} = 0$,\n",
    "\n",
    "from which we obtain \n",
    "\n",
    " $\\hat{s} = \\sum_{i=1}^N \\frac{(f_i - 1)(m_{0,i}-1) }{\\sigma_i^2} \\left[ \\sum_{i=1}^N \\frac{(m_{0,i}-1)^2}{\\sigma_i^2}\\right]^{-1}$\n",
    " \n",
    " and \n",
    " \n",
    " $\\hat{b} = \\sqrt{\\chi^2/N}$\n",
    " \n",
    " For the standard errors on these parameters we use \n",
    " \n",
    " $\\sigma_s^{-2}  = -\\frac{\\partial^2\\ln{\\cal L}}{\\partial^2 s^2} \\left.\\right|_{\\hat{s},\\hat{b}}$ \n",
    " \n",
    " and\n",
    " \n",
    " $\\sigma_b^{-2}  = -\\frac{\\partial^2\\ln{\\cal L}}{\\partial^2 b^2} \\left.\\right|_{\\hat{s},\\hat{b}}$ \n",
    "\n",
    "to derive\n",
    "\n",
    "$\\sigma_s = b\\left[\\sum_{i=1}^N \\frac{(m_i-1)^2}{\\sigma_i^2}\\right]^{-1/2}$\n",
    " \n",
    " and\n",
    " \n",
    " $\\sigma_b  = \\left[3\\chi^2/b^4 -N/b^2\\right]^{-1/2} .$\n",
    " \n",
    "**N.B.** Whether or how much of the data outside transit to include depends on whether these data can be assumed to have the same noise characteristics as the data in transit. Note that including these data has no effect on $s$ or $\\sigma_s$, because of the factors $(m_{0,i}-1)$ in their calculation, but will affect the estimates of $b$ and $\\sigma_b$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "from pycheops.models import TransitModel, scaled_transit_fit\n",
    "from pycheops.instrument import transit_noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simulate some data\n",
    "\n",
    "The reference model ${\\bf m_0}$ wil be a transit depth of $D = 250\\,$ppm and a transit width $W=3\\,$hours.  \n",
    "Use limb darkening appropriate for a K-type dwarf.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 20936 characters]\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1295c0a90>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "WIDTH_HOURS = 3\n",
    "DEPTH_PPM = 250\n",
    "PERIOD = 13\n",
    "T_0 = 0.5\n",
    "SIGMA_TRUE_PPM = 150\n",
    "SIGMA_NOMINAL_PPM = 100\n",
    "\n",
    "h_1 = 0.680\n",
    "h_2 = 0.432\n",
    "\n",
    "time = np.linspace(0,1,1440)\n",
    "\n",
    "D = DEPTH_PPM*1e-6  # Set transit depth\n",
    "W = WIDTH_HOURS/24/PERIOD  # Transit Width in phase units\n",
    "b = 0 # Impact parameter S is the ratio (W_23 / W_14)**2 = ((1-k**2)-b**2)/((1+k)**2-b**2) \n",
    "k = sqrt(D)  # Ratio of the radii\n",
    "S = ((1-k)**2-b**2)/((1+k)**2-b**2)  # Transit Shape\n",
    "\n",
    "TM = TransitModel()\n",
    "pars = TM.make_params(T_0=T_0, P=PERIOD, D=D, W=W, S=S, h_1=h_1, h_2=h_2)\n",
    "model_0 = TM.eval(params=pars, t=time) \n",
    "\n",
    "flux = model_0 + np.random.normal(0, SIGMA_TRUE_PPM*1e-6, size=len(time))\n",
    "\n",
    "fig,ax=subplots(1,1,figsize=(5,4))\n",
    "ax.plot(time, flux-1,'b.',ms=1)\n",
    "ax.plot(time, model_0-1)\n",
    "ax.set_xlim(min(time),max(time))\n",
    "ax.set_xlabel(\"Time\")\n",
    "ax.set_ylabel(\"Flux - 1\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit the transit depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Transit scale factor s = 0.98 +/- 0.04\n",
      "Error scale factor  b = 1.49 +/- 0.03\n",
      "Depth = 245 +/- 10 ppm\n"
     ]
    }
   ],
   "source": [
    "sigma = SIGMA_NOMINAL_PPM*1e-6\n",
    "s, b, sigma_s, sigma_b = scaled_transit_fit(flux, sigma ,model_0)\n",
    "print('Transit scale factor s = {:0.2f} +/- {:0.2f}'.format(s, sigma_s))\n",
    "print('Error scale factor  b = {:0.2f} +/- {:0.2f}'.format(b, sigma_b))\n",
    "depth = s*DEPTH_PPM\n",
    "e_depth = sigma_s*DEPTH_PPM\n",
    "print('Depth = {:0.0f} +/- {:0.0f} ppm'.format(depth, e_depth))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&copy; Pierre Maxted, Keele University, 2019"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
