{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Making Custom Distributions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By using the InterpolatedUnivariateDistribution class, you can easily create a single-variable distribution by specifying its PDF as a callable function. Here, we'll demonstrate this functionality by implementing the asymmetric Lorentz distribution of [Stancik and Brauns](http://www.sciencedirect.com/science/article/pii/S0924203108000453)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Preamble" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As always, we start by setting up the Python environment for inline plotting and true division." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/cgranade/anaconda/envs/qinfer-binder/lib/python3.5/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.\n", " warnings.warn(self.msg_depr % (key, alt_key))\n" ] } ], "source": [ "from __future__ import division\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "try: plt.style.use('ggplot')\n", "except: pass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we import the InterpolatedUnivariateDistribution class." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/cgranade/anaconda/envs/qinfer-binder/lib/python3.5/site-packages/qinfer/metrics.py:51: UserWarning: Could not import scikit-learn. Some features may not work.\n", " warnings.warn(\"Could not import scikit-learn. Some features may not work.\")\n", "/home/cgranade/anaconda/envs/qinfer-binder/lib/python3.5/site-packages/IPython/parallel.py:13: ShimWarning: The IPython.parallel package has been deprecated. You should import from ipyparallel instead.\n", " \"You should import from ipyparallel instead.\", ShimWarning)\n", "/home/cgranade/anaconda/envs/qinfer-binder/lib/python3.5/site-packages/qinfer/parallel.py:53: UserWarning: Could not import IPython parallel. Parallelization support will be disabled.\n", " \"Could not import IPython parallel. \"\n" ] } ], "source": [ "from qinfer.distributions import InterpolatedUnivariateDistribution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Defining Distributions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The asymmetric Lorentz distribution is defined by letting the scale parameter $\\gamma$ of a Lorentz distribution be a function of the random variable $x$,\n", "$$\n", " \\gamma(x) = \\frac{2\\gamma_0}{1 + \\exp(a [x - x_0])}.\n", "$$\n", "It is straightforward to implement this in a vectorized way by defining this function and then substituting it into the PDF of a Lorentz distribution." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def asym_lorentz_scale(x, x_0, gamma_0, a):\n", " return 2 * gamma_0 / (1 + np.exp(a * (x - x_0)))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def asym_lorentz_pdf(x, x_0, gamma_0, a):\n", " gamma = asym_lorentz_scale(x, x_0, gamma_0, a)\n", " return 2 * gamma / (np.pi * gamma_0 * (1 + 4 * ((x - x_0) / (gamma_0))**2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Once we have this, we can pass the PDF as a lambda function to InterpolatedUnivariateDistribution in order to specify\n", "the values of the location $x_0$, the nominal scale $\\gamma_0$ and the asymmetry parameter $a$." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dist = InterpolatedUnivariateDistribution(lambda x: asym_lorentz_pdf(x, 0, 1, 2), 2, 1200)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The resulting distribution can be sampled like any other, such that we can quickly check that it produces something of the desired shape." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "hist(dist.sample(n=10000), bins=100);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We note that making this distribution object is fast enough that it can conceivably be embedded within a likelihood function itself, so as to enable using the method of hyperparameters to estimate the parameters of the asymmetric Lorentz distribution." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1000 loops, best of 3: 764 µs per loop\n" ] } ], "source": [ "%timeit dist = InterpolatedUnivariateDistribution(lambda x: asym_lorentz_pdf(x, 0, 1, 2), 2, 120)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1000 loops, best of 3: 957 µs per loop\n" ] } ], "source": [ "%timeit dist = InterpolatedUnivariateDistribution(lambda x: asym_lorentz_pdf(x, 0, 1, 2), 2, 1200)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "100 loops, best of 3: 2.65 ms per loop\n" ] } ], "source": [ "%timeit dist = InterpolatedUnivariateDistribution(lambda x: asym_lorentz_pdf(x, 0, 1, 2), 2, 12000)" ] }, { "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.5.1" }, "widgets": { "state": {}, "version": "1.1.1" } }, "nbformat": 4, "nbformat_minor": 0 }