Lab_Lec3_sherpa_fit.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "from sherpa.astro import ui"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ui.load_pha('3c273.pi')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ui.plot_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ui.notice(0.1, 6.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ui.subtract()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ui.plot_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ui.set_xlog()\n",
    "ui.set_ylog()\n",
    "ui.plot_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ui.list_models()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sherpa.astro import *\n",
    "from sherpa.models import *\n",
    "#from sherpa.models.astro.xspec import get_xsversion\n",
    "#from sherpa.models import xspec\n",
    "#ui.set_source(PowLaw1d.p1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "truth = PowLaw1D()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ui.guess(truth.p1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's try the image\n",
    "Fitting an image using a PSF image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ui.load_image(\"center_box_0.25pix.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ui.image_getregion()\n",
    "#ui.notice2d(\"rotbox(88.16875,75.8625,70.416667,68.508334,0)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ui.load_psf(\"psf0\",\"psf_f1_norm_0.25pix.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ui.set_psf(\"psf0\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psf0.center = [128,129]\n",
    "psf0.size = [72,72]\n",
    "ui.set_psf(psf0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ui.get_psf())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We now set a Gaussian source model with a constant background and thaw ellipticity \n",
    "# and rotation angle. We also set the limits for background\n",
    "ui.set_source(ui.gauss2d.g1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g1.fwhm = 6\n",
    "g1.xpos = 90\n",
    "g1.ypos = 80\n",
    "g1.ampl = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we set the statistic to the Cash statisitc, and the optimization method to Simplex.\n",
    "ui.set_stat('cstat')\n",
    "ui.set_method('simplex')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let us fit, by finding a MAP estimate. First, we have sherpa guessing initial values\n",
    "ui.guess(g1)\n",
    "ui.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let us fit again, just in case the optimizer got stuck\n",
    "ui.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We now print the model\n",
    "print(ui.get_source())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let us add background\n",
    "ui.set_source(g1 + ui.const2d.c1)\n",
    "c1.c0 = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ui.freeze(c1.c0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ui.get_source())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit\n",
    "ui.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do it again\n",
    "ui.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fit is run again with the background component thawed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ui.thaw(c1.c0)\n",
    "ui.fit() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we will get the confidence level of the fit, and in doing so, recalculate the best-fit parameters using the confidence method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ui.conf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let us get a covariance matrix\n",
    "ui.covar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we run MCMCM\n",
    "draws = ui.get_draws(niter=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svals, accept, pvals = draws"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Acceptance rate\n",
    "accept.sum() * 1.0 / 5000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "# Now let us plot the traces\n",
    "plt.subplots(3,2,figsize=(10,10))\n",
    "plt.subplot(3,2,1)\n",
    "plt.plot(pvals[0, :])\n",
    "plt.xlabel('Iteration')\n",
    "plt.ylabel('fwhm')\n",
    "plt.subplot(3,2,2)\n",
    "plt.plot(pvals[1, :])\n",
    "plt.xlabel('Iteration')\n",
    "plt.ylabel('xpos')\n",
    "plt.subplot(3,2,3)\n",
    "plt.plot(pvals[2, :])\n",
    "plt.xlabel('Iteration')\n",
    "plt.ylabel('ypos')\n",
    "plt.subplot(3,2,4)\n",
    "plt.plot(pvals[3, :])\n",
    "plt.xlabel('Iteration')\n",
    "plt.ylabel('ampl')\n",
    "plt.subplot(3,2,5)\n",
    "plt.plot(pvals[4, :])\n",
    "plt.xlabel('Iteration')\n",
    "plt.ylabel('bkg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sherpa import plot\n",
    "pdf = plot.PDFPlot()\n",
    "pdf.prepare(pvals[0, :], 20, False, 'fwhm', name='fwhm')\n",
    "pdf.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdf.prepare(pvals[1, :], 20, False, 'xpos', name='xpos')\n",
    "pdf.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdf.prepare(pvals[2, :], 20, False, 'ypos', name='ypos')\n",
    "pdf.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdf.prepare(pvals[3, :], 20, False, 'ampl', name='ampl')\n",
    "pdf.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdf.prepare(pvals[4, :], 20, False, 'bkg', name='bkg')\n",
    "pdf.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}