Testing bounding ellipsoid in nestle

estimated_ellipsoid_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import nestle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def random_ellipsoid(n):\n",
    "    \"\"\"Return a random `n`-d ellipsoid centered at the origin\"\"\"\n",
    "\n",
    "    # `a` in the ellipsoid must be positive definite, so we have to construct\n",
    "    # a positive definite matrix. For any real, non-singular matrix A,\n",
    "    # `A^T A` will be positive definite.\n",
    "    det = 0.\n",
    "    while abs(det) < 1.e-10:  # ensure a non-singular matrix\n",
    "        A = np.random.rand(n, n)\n",
    "        det = np.linalg.det(A)\n",
    "\n",
    "    return nestle.Ellipsoid(np.zeros(n), np.dot(A.T, A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "SQRTEPS = math.sqrt(float(np.finfo(np.float64).eps))\n",
    "\n",
    "def est_ellipsoid(x):\n",
    "    \"\"\"Like bounding ellipsoid, but don't scale to include the outermost point.\n",
    "    This is useful for comparing to the result of bounding_ellipsoid to see if\n",
    "    the outermost point is usually outside the 'estimated' ellipsoid.\"\"\"\n",
    "    \n",
    "    npoints, ndim = x.shape\n",
    "    \n",
    "    ctr = np.mean(x, axis=0)\n",
    "    delta = x - ctr\n",
    "    cov = np.atleast_2d(np.cov(delta, rowvar=0)) * (ndim + 2)\n",
    "    \n",
    "    a = np.linalg.inv(cov)\n",
    "\n",
    "    return nestle.Ellipsoid(ctr, a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ndim    true_vol       est_vol  bounding_vol\n",
      "----    --------       -------  ------------\n",
      " 1        7.0890        7.1841        7.1841\n",
      " 2       15.8982       16.7563       16.7563\n",
      " 3      163.6792      152.9843      186.6362\n",
      " 4       26.5588       27.3430       33.5129\n",
      " 5      114.1826      109.4350      228.6318\n",
      " 6      253.1310      235.7950      416.4938\n",
      " 7      145.5416      129.0424      355.7259\n",
      " 8       58.3859       45.3929      142.8048\n",
      " 9      635.2109      519.4966     1026.8658\n",
      "10       47.5843       43.4041      275.6352\n",
      "11       57.1798       39.7135      153.4214\n",
      "12       80.2361       75.3141      219.1016\n",
      "13        5.9599        4.5178       20.2549\n",
      "14      386.6646      225.3474     1055.4631\n",
      "15       49.9046       28.9416      127.8040\n",
      "16      165.1167       81.5517      587.5828\n",
      "17        3.7018        1.3528       25.6952\n",
      "18        1.4007        0.5769       10.0011\n",
      "19        0.3621        0.1116        1.3829\n",
      "20        0.7569        0.2576        2.4605\n"
     ]
    }
   ],
   "source": [
    "# For dimensions 1-20, draw points uniformly within a random ellipsoid,\n",
    "# calculate the bounding ellipsoid and the 'estimated ellipsoid', which\n",
    "# is the bounding ellipsoid before scaling to include the outermost point.\n",
    "# Compare the volumes of the true ellipsoid, and the two bounding ellipsoids.\n",
    "\n",
    "npoints = 100\n",
    "\n",
    "print(\"ndim    true_vol       est_vol  bounding_vol\")\n",
    "print(\"----    --------       -------  ------------\")\n",
    "for ndim in range(1, 21):\n",
    "    ell_gen = random_ellipsoid(ndim)\n",
    "    x = ell_gen.samples(npoints)\n",
    "    ell1 = est_ellipsoid(x)\n",
    "    ell2 = nestle.bounding_ellipsoid(x)\n",
    "\n",
    "    print(\"{:2d}  {:12.4f}  {:12.4f}  {:12.4f}\"\n",
    "          .format(ndim, ell_gen.vol, ell1.vol, ell2.vol))\n"
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}