Sample Notebook

gistfile1.json

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from scipy.linalg import inv"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def additive_log_ratio(x):\n",
      "    y = np.log(x / x[-1])\n",
      "    return y[:-1]\n",
      "\n",
      "def make_fcns(alpha, mu, sigma, sigma_zero, mu_zero):\n",
      "    # alpha, mu, vectors\n",
      "    # sigma, matrix\n",
      "    # sigma_zero, matrix\n",
      "    # mu_zero, vector\n",
      "    # returns H(x), f1(x), f2(x)\n",
      "    \n",
      "    \n",
      "    def H(x):\n",
      "        y = additive_log_ratio(x)\n",
      "        val = np.exp((1/4 * alpha.T * sigma * alpha * (y - mu).T * inv(sigma) * (y - mu)\n",
      "                      + 1/2 * ((y - mu_zero).T * inv(sigma_zero) * (y - mu_zero))))\n",
      "        return val\n",
      "\n",
      "    def f1(x):\n",
      "        y = additive_log_ratio(x)\n",
      "        a = - 1/2 * ((y - mu_zero).T * inv(sigma_zero) * (y - mu_zero))\n",
      "        val = np.exp(a)\n",
      "        \n",
      "        return val\n",
      "    \n",
      "    def f2(x):\n",
      "        inv_prod = 1.0\n",
      "        for i in range(len(x)):\n",
      "            inv_prod *= 1/x[i]\n",
      "        return inv_prod\n",
      "    \n",
      "    \n",
      "    def h(x):\n",
      "        D = len(x)\n",
      "        val = np.zeros((1/2 * D * (D+1)))\n",
      "        val[:D] = np.log(x)\n",
      "        \n",
      "        s = D\n",
      "        for i in range(D-1):\n",
      "            for j in range(i+1, D):\n",
      "                val[s] =-1/2*(log(x[i]) - log(x[j]))**2\n",
      "                s += 1\n",
      "        val = np.matrix(val)\n",
      "        return val\n",
      "    \n",
      "    \n",
      "    def h_prime(x):\n",
      "        val = h(x).T\n",
      "        return val\n",
      "                                    \n",
      "    \n",
      "    def integrand_log_partition(x):\n",
      "        _H = H(x)\n",
      "        _f1 = f1(x)\n",
      "        _f2 = f2(x)\n",
      "        print \"H: %.3f\" %_H\n",
      "        print \"f1: %.3f\" %_f1\n",
      "        print \"f2: %.3f\" %_f2\n",
      "        val = _H * _f1 * _f2\n",
      "        return val\n",
      "    \n",
      "    \n",
      "    def integrand_hess_log_partion(x):\n",
      "        return integrand_log_partition(x) * h(x) * h_prime(x)\n",
      "\n",
      "    return integrand_hess_log_partion"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Stupid Example:\n",
      "alpha = np.matrix([1, 0, 0]).T\n",
      "\n",
      "mu = np.zeros_like(alpha)\n",
      "mu_zero = np.zeros_like(alpha)\n",
      "\n",
      "sigma_zero = np.eye(3)\n",
      "sigma = 4*sigma_zero\n",
      "\n",
      "x = np.array([0.25, 0.25, 0.25, 0.25])\n",
      "x /= np.sum(x)\n",
      "x = np.matrix(x).T"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "test_fun = make_fcns(alpha, mu, sigma, sigma_zero, mu_zero)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#print x\n",
      "test_fun(x)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "H: 1.000\n",
        "f1: 1.000\n",
        "f2: 256.000\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "matrix([[ 256.]])"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}