{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 23"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm, ncx2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ks = [5, 10, 100]\n",
    "sim = 100_000\n",
    "\n",
    "mplot = 200\n",
    "thetanorm2 = np.linspace(.0001, 10, mplot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exact formula for expectation from [1]:\n",
    "$$\\mathbb{E}[1/\\lVert Y\\rVert^2] = \\frac12\\exp\\{-\\lVert\\theta\\rVert^2/2\\}\\frac{\\Gamma(k/2-1)}{\\Gamma(k/2)}{}_1F_1(k/2-1, k/2, \\lVert\\theta\\rVert^2/2)$$\n",
    "\n",
    "1. A. K. Gupta and D. K. Nagar. *Matrix Variate Distributions*. Chapman and Hall, 1999."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.special import hyp1f1, gamma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simulations and plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, sharey=True, figsize=(12, 6))\n",
    "\n",
    "for i in range(3):\n",
    "    k = ks[i]\n",
    "    chi2mean = .5 * np.exp(-thetanorm2/2) * gamma(k/2 - 1)/gamma(k/2) * hyp1f1(k/2 - 1, k/2, thetanorm2/2)\n",
    "    \n",
    "    R = 1 - (k-2)**2/k*chi2mean\n",
    "    \n",
    "    Y2sample = ncx2(k, thetanorm2).rvs((sim, mplot))\n",
    "    D = (k**(-1)*Y2sample - 2)*(Y2sample <= k-2) + (-(k-2)**2/k*1/Y2sample)*(Y2sample > k-2)\n",
    "    Dmean = np.mean(D, axis=0)\n",
    "    Rstein = 1 + Dmean\n",
    "    \n",
    "    axs[i].set_title(f'k = {ks[i]}')\n",
    "    axs[i].set_xlabel(r'$||\\theta||^2$')\n",
    "    axs[i].plot(thetanorm2, np.ones(mplot), label='Y', color='black', linestyle='dotted')\n",
    "    axs[i].plot(thetanorm2, R, label=r'$\\hat\\theta$', color='black', linestyle='dashed')\n",
    "    axs[i].plot(thetanorm2, Rstein, label='Stein', color='black')\n",
    "\n",
    "axs[0].set_ylabel('risk')\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}