{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# MultiNorm\n",
"\n",
"This is a tutorial notebook for the ``MultiNorm`` class from the ``multinorm`` Python package ([Github](https://github.com/cdeil/multinorm)).\n",
"\n",
"- A quick tutorial introduction (~ 10 min read)\n",
"- Explain and illustrate some of the theory / statistical aspects\n",
"- Hands-on introduction how to work with ``MultiNorm`` class.\n",
"\n",
"Suggestions:\n",
"\n",
"- Read this tutorial and judge if `multinorm` is useful for you.\n",
"- If yes, install it with `pip install multinorm`, download this notebook, and execute / play with these examples on your machine.\n",
"- Learn more at https://multinorm.readthedocs.io (complete description of features and theory).\n",
"\n",
"We will use [pandas](http://pandas.pydata.org), so after reading this you will be a data scientist!\n",
"In the future we will probably change the example model to a deep neural network and use [tensorflow](https://www.tensorflow.org/), because likelihood fitting is so 1990s."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n",
"\n",
"- TODO: What is MVN (theory)?\n",
"- https://en.wikipedia.org/wiki/Multivariate_normal_distribution\n",
"- TODO: What is `MultiNorm` (code)?\n",
"\n",
"\n",
"We will be using [scipy.optimize.curve_fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) to fit a model to data and obtain an example covariance matrix. Note that most analysis problems are not a simple least squares curve fit, but involve custom data and likelihood. Information how to obtain parameter error and covariance estimates with other popular Python packages is described in the `multinorm` documentation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Imports\n",
"\n",
"Let's start by importing everything we'll use in this notebook."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"import scipy.optimize\n",
"from multinorm import MultiNorm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example Setup\n",
"\n",
"As mentioned in the introduction, the kind of analysis we are considering here is that there is a **model** with **parameters**, as well as some **data** and a **likelihood function**. \n",
"\n",
"We will use the [Rat42](https://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml) dataset from the [NIST Standard Reference Database](https://www.itl.nist.gov/div898/strd/). To quote from [here](https://www.itl.nist.gov/div898/strd/nls/data/LINKS/i-ratkowsky2.shtml):\n",
"\n",
"> This model and data are an example of fitting sigmoidal growth curves taken from Ratkowsky (1983).\n",
"> The response variable is pasture yield, and the predictor variable is growing time.\n",
"\n",
"Did I hear you say\n",
"\n",
"> Why? pasture yield? WTF?\n",
"\n",
"Well, the data is small (9 observations x, y) and , and the model is simple (one line in Python) and has 3 parameters, which is nice to illustrate properties of the correlation matrix and multivariate normal distribution. And the fit is using a least square likelihood, which is built-in in `scipy` and available for pretty much any fitting package, or can be implemented also with one line in Python."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
0
\n",
"
1
\n",
"
2
\n",
"
3
\n",
"
4
\n",
"
5
\n",
"
6
\n",
"
7
\n",
"
8
\n",
"
\n",
" \n",
" \n",
"
\n",
"
x
\n",
"
9.00
\n",
"
14.0
\n",
"
21.00
\n",
"
28.00
\n",
"
42.00
\n",
"
57.00
\n",
"
63.00
\n",
"
70.00
\n",
"
79.00
\n",
"
\n",
"
\n",
"
y
\n",
"
8.93
\n",
"
10.8
\n",
"
18.59
\n",
"
22.33
\n",
"
39.35
\n",
"
56.11
\n",
"
61.73
\n",
"
64.62
\n",
"
67.08
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" 0 1 2 3 4 5 6 7 8\n",
"x 9.00 14.0 21.00 28.00 42.00 57.00 63.00 70.00 79.00\n",
"y 8.93 10.8 18.59 22.33 39.35 56.11 61.73 64.62 67.08"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data = pd.DataFrame({\n",
" 'x': [9, 14, 21, 28, 42, 57, 63, 70, 79],\n",
" 'y': [8.93, 10.80, 18.59, 22.33, 39.35, 56.11, 61.73, 64.62, 67.08],\n",
"})\n",
"data.T"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And having a DOI (http://dx.doi.org/10.18434/T43G6C) and expected reference results for the analysis is also nice. This is the result we should find:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
value
\n",
"
error
\n",
"
\n",
"
\n",
"
par
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
b1
\n",
"
72.462238
\n",
"
1.734028
\n",
"
\n",
"
\n",
"
b2
\n",
"
2.618077
\n",
"
0.088295
\n",
"
\n",
"
\n",
"
b3
\n",
"
0.067359
\n",
"
0.003447
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" value error\n",
"par \n",
"b1 72.462238 1.734028\n",
"b2 2.618077 0.088295\n",
"b3 0.067359 0.003447"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"expected = pd.DataFrame.from_records(\n",
" [\n",
" (\"b1\", 7.2462237576e01, 1.7340283401e00),\n",
" (\"b2\", 2.6180768402e00, 8.8295217536e-02),\n",
" (\"b3\", 6.7359200066e-02, 3.4465663377e-03),\n",
" ],\n",
" columns=[\"par\", \"value\", \"error\"],\n",
" index='par',\n",
")\n",
"expected"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Model\n",
"\n",
"Our model is just a simple line."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def model(x, b1, b2, b3):\n",
" return b1 / (1 + np.exp(b2 - b3 *x))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Fit\n",
"\n",
"If we fit the model to the data, we obtain a parameter estimate (`popt`, the \"optimial\" or \"best\" parameter values) and a covariance matrix (`pcov`)."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"p0 = (100, 1, 0.1)\n",
"popt, pcov = scipy.optimize.curve_fit(\n",
" f=model,\n",
" xdata=data[\"x\"],\n",
" ydata=data[\"y\"],\n",
" p0=p0,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's make a plot to see the best-fit model and the data:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xl8FPX9x/HXJyFAOMNtCGBQKYoiohFRUVG08UCJB5coIFhaq63WFpX+2moPq/6wtvZXa4sXFLlvRAEpxVptPUBuMAUFgQQICEGOQJbk+/tjFg0STAK7md3Z9/PxyGN3JrPkvY+BN5PvfmfGnHOIiEj8S/I7gIiIRIYKXUQkIFToIiIBoUIXEQkIFbqISECo0EVEAkKFLiISECp0EZGAUKGLiAREjer8YU2bNnWZmZnV+SNFROLekiVLdjrnmlW0XbUWemZmJosXL67OHykiEvfM7LPKbKchFxGRgFChi4gEhApdRCQgVOgiIgGhQhcRCQgVuohIQKjQRUQColrnoYuIBFqoCPbvhAOfw4GdcGBXeHknXPIDSG0U1R9fYaGbWXtgUplVpwG/AP4WXp8JbAT6OOd2Rz6iiIiPnPOKeW8+fLHVe9y7HfYXwL7tsG9H+HkBFO8r/8+wZOjY2/9Cd87lAucBmFkykAfMAB4BFjrnnjSzR8LLD0cxq4hI5BXvh8LNsGcL7NlU5vlm+CIP9m6DkuJjX1c7Deq1gHrNoWVnqNsc6jWDOk2hblOo0yT8vAnUaghJ0R/hruqQSw/gE+fcZ2bWC+geXj8GeAsVuojEov2fw+frYfcG2LWhzONG7+i6rKQa0KAlNGwNrS+C+unecv308PN0r8hr1PLlrXyTqhZ6P2BC+HkL59xWAOfcVjNrHtFkIiJVUVrqFfWOXPh8Hez8L+xc7z0W7SqzoUGDDGjcFr6VDY0yIe1USGsNDVt5pZ2U7Ne7OCmVLnQzqwncBIyoyg8ws2HAMIA2bdpUKZyISLkO7IJtK6FgDWxf7X3t+BhCB77apm5zaNoOOtwETdp5zxu1hbQ2kFLbv+xRVJUj9OuAj5xz28PL280sPXx0ng4UlPci59woYBRAVlaWO6m0IpJ49n8OW5dC/jLYugzyl3tj3UfUaQotOsD5g7zHZmd55Z2a5l9mn1Sl0Pvz1XALwGxgEPBk+HFWBHOJSCIqCcH2VbD5g6++ypZ349OgVRZ0uRtO6QgtzvE+lBSgkoVuZnWAa4Dvlln9JDDZzIYCm4DekY8nIoFWvB82vQefveuVd96Sr4ZN6reE1l2gy3eg5XlwyrkJedRdFZUqdOfcAaDJ19Z9jjfrRUSkckIHYcuHsPFfsOFt2LIYSkPePO30c+H8gV6Jt77I+4BSqkRniopIdH3+CaxbAOve9I7EDx8ES4L0TnDx9yHzcmjTFWrV8ztp3FOhi0hkhQ7CZ+/AugXsW/UG9fZ7Y+CfWQaH2/bm9C7Xw6mXaPgkClToInLyivd7R+FrZ8N/50PxPkqSarGkpAMLD1/JW6Xnscm1IHVdMk907EiOyjwqVOgicmIOfgH/nQdrZsH6hXC4yDvd/Zxb4cyeXD2thA0HSo96SVGohJHzc8npnOFT6GBToYtI5ZWEvPJeMRFy53rj4fVOgc53eCfwtLkEkr1a2bjn9XL/iPzCoupMnFBU6CLyzZyD/KWwYhKsnOpdCja1sTcj5ZzboNWF5V54qmVaKnnllHfLtNTqSJ2QVOgiUr6i3bB8IiwZ7Z1Wn1wT2l8H5/aDM66GGjW/8eXDs9szYvpKikIlX65LTUlmeHb7KAdPXCp0EfmKc5D3ESx+GVZN88bFMy6Ann+As3OqdD3vI+PkI+fnkl9YRMu0VIZnt9f4eRSp0EUSzMyleceW7NmNYOVkr8i3LoeUutCpL2QN8eaLn6Cczhkq8GqkQhdJIDOX5h01DFJUuJ0tM8ZxaO5CahXvhuYd4Pqn4dy+ULuBz2mlqlToIglk5PxcikIlnGrbuDv5DXon/5PaFuKdwxfSbfBjcOqlYOZ3TDlBKnSRBNJ4z2p+ljKL7KTFhEhmRkk3Xii5gU9dBhsyu/kdT06SCl0kEWxdDm89yWu13mCPq8OfS25izOFsduCdsZmhqYSBoEIXCbJtq+CtJ+DjOVC7IWvO/CEDV3dm5+Gv7oepqYTBoUIXCaKCj+Gt33qn5ddqAN1HwEXfo0NqGj8rb5aLZqIEggpdJEj27fCKfMlob+rh5Q95l6gtM39cUwmDS4UuEgShg/Den+Ffz3h3/LnwO3DFw1C3ScWvlcBQoYvEsHJPAip7dO2cd0bn33/p3Xuz/fVwza+8myRLwlGhi8Sor58ElFdYxIjpK4HwafVbV8DrD3q3dDulI+S8Bm0v9zOy+EyFLhKjjpwEVFZRqIQ/zVtGzvbn4P3nvase9noOOvWHpGSfkkqsUKGLxKhjrxvu+HbSYh47OAbe2wUX3AVXP1qlC2ZJsFWq0M0sDXgROAdwwBAgF5gEZAIbgT7Oud1RSSmSgMpeTzyDHTyWMoZrkj9ivZ0KQyZB6y4+J5RYc+xV6cv3LDDPOXcm0AlYCzwCLHTOtQMWhpdFJEKGZ7enTopxZ/KbLKj1EJcmreap0jtY3XO2ylzKVeERupk1AC4HBgM454qBYjPrBXQPbzYGeAt4OBohRRJRTttSLmnxLM13vsc/S87l2Tr3MfDabvTSHHI5jsoMuZwG7ABeMbNOwBLgfqCFc24rgHNuq5k1j15MkQTiHCx9FeaNoDkObnyWK84fxBW6CqJUoDJDLjWA84HnnXOdgf1UYXjFzIaZ2WIzW7xjx44TjCmSIPZug/F9YfZ93o0l7nkXLhisS9pKpVSm0LcAW5xz74eXp+IV/HYzSwcIPxaU92Ln3CjnXJZzLqtZs2aRyCwSTGtfg+cugg3/hGufhEGvQaNMv1NJHKmw0J1z24DNZnbkcmw9gDXAbGBQeN0gYFZUEooEXeggvDEcJt0BjdvC996BrvdAUmXnLIh4KjsP/QfAODOrCXwK3IX3n8FkMxsKbAJ6RyeiSIB9/glMGQzbVsDF90GPR6FGTb9TSZyqVKE755YBWeV8q0dk44gkkBVTYM4DkJwC/SdC++v8TiRxTmeKilS34gMw9yFYOhZad4XbXoKGrfxOJQGgQhepTrs3wsQBsH01XPYT78YTyfpnKJGhv0ki1eXTt7zxclcKA6ZCu6v9TiQBo0IXiYKjrmPesDZ/OeM9Oq55Gpq2h37joMnpfkeUAFKhi0RY2euY16KYnxx4jo6r3yU//RpaDn4FatX3O6IElApdJMKOXMe8JTv5a81nONs+Y2SoDzN39eNdlblEkQpdJMLyC4voZOt5sebT1CLE0NBPWFTaGdtz0O9oEnA6FU0kwvrWX8HEmr+hyNXi5uJfsai0M+Bd31wkmnSELhJJ7z3PE6GnWMHpDCn+MZ/TEIDUlGSGZ7ev4MUiJ0dH6CKRUFoCcx+GeY9gZ97App6TqJ12CgZkpKXyxC0dvRs7i0SRjtBFTlbxfph2N+S+4V2P5ZpfcWNSMjdeeIbfySTBqNBFTsb+nTCuN2xdBtc/DV2+43ciSWAqdJETtWcL/C0H9myGfuN1cS3xnQpd5ETsXOeV+aEv4M4ZcOolficSUaGLVFn+Mnj1FrAkGDzHu1WcSAzQLBeRqtj4DozuCSl14K55KnOJKSp0kcrKnQuv3goN0mHIfGiqWSwSW1ToIpWxapp3HfPmZ3lH5g01p1xij8bQRSqyYgrMGAatL4IBU3S1RIlZOkIX+SbLJ3ll3uYS76YUKnOJYSp0keNZNh5mfBcyu8GAyVCrnt+JRL5RpYZczGwjsBcoAQ4757LMrDEwCcgENgJ9nHO7oxNTpJp9NBZm/wBO6+6dNFSzjt+JRCpUlSP0K51z5znnssLLjwALnXPtgIXhZZH4t2Q0zL4PTr8K+k9QmUvcOJkhl17AmPDzMUDOyccR8dmS0fDa/dDu296ReYquYS7xo7KF7oA3zWyJmQ0Lr2vhnNsKEH5sHo2AItVm+UR47QGvzPu+Cim1/U4kUiWVnbZ4qXMu38yaAwvM7OPK/oDwfwDDANq0aXMCEUWqweoZMPMeaHs59BkLNWr5nUikyip1hO6cyw8/FgAzgC7AdjNLBwg/FhzntaOcc1nOuaxmzZpFJrVIJOXO865n3voib8xcR+YSpyosdDOra2b1jzwHvg2sAmYDg8KbDQJmRSukSNR88g+YfCecci7cPhlq1vU7kcgJq8yQSwtghpkd2X68c26emX0ITDazocAmoHf0YopEwcZ3YcLt0PRbcMc0qN3A70QiJ6XCQnfOfQocc0k559znQI9ohBKJui1LYHwfSGsDd86EOo39TiRy0nSmqCSego9h3K1QtykMnAX19NmOBIMKXRJL4SYYezMk1/LKvEG634lEIkZXW5TEsX+nV+ah/XDXXGiU6XcikYhSoUtiOLTXuznFnjwYOBNanO13IpGIU6FL8IUOwsTbYdtKb555m65+JxKJChW6BFtpCUy/Gza8DTePgm9l+51IJGr0oagEl3Pw+oOw9jXIfgI69fU7kUhUqdAluN4e6V09sduDcPH3/U4jEnUqdAmmpa/CosehU3/o8Qu/04hUCxW6BM+6v8PsH8JpV8KNfwTvshUigadCl2DJXwqTB0KLDtB3LNSo6XcikWqjQpfg2L0RxvWBOk1gwFSoVd/vRCLVStMWJRgO7IJXb4OSYhg8B+qf4ncikWqnQpf4FyqCCf2867QMnAXN2vudSMQXKnSJb6WlMON7sPkD6D0aTr3Y70QivlGhS3xb+EtYMxOu+TWcneN3GhFf6UNRiV9LRsO7f4CsIXDJD/xOI+I7FbrEp/ULYc6DcMbVcN1IzTUXQYUu8Wj7apg8CJqfBbe9AskaORQBFbrEm73bvLnmterB7ZN1Y2eRMnRoI/GjeD+M7wtFu2HIXGiY4XcikZhS6SN0M0s2s6VmNie83NbM3jezdWY2ycx0jrVET2kpTB8G21bAbS9Deie/E4nEnKoMudwPrC2z/BTwe+dcO2A3MDSSwUSOsvAx+HiOd13z9tf6nUYkJlWq0M2sFXAD8GJ42YCrgKnhTcYAmgQsETdzaR5P/Oan8O6zTEu+lpk1e/odSSRmVXYM/Q/AQ8CRqx01AQqdc4fDy1uAcgc0zWwYMAygTZs2J55UEs7MpXnMmD6BF5P+wtulHXno4ABqzlgFZuR01vi5yNdVeIRuZj2BAufckrKry9nUlfd659wo51yWcy6rWbNmJxhTEtHEuf/g2aRn2OBO4d7Q/ZSQTFGohJHzc/2OJhKTKnOEfilwk5ldD9QGGuAdsaeZWY3wUXorID96MSXhHNjFEwd/Q4klMSQ0nL3U+fJb+YVFPgYTiV0VHqE750Y451o55zKBfsA/nHMDgEXAbeHNBgGzopZSEsvhYph0JxlJOxlW/CBbXPOjvt0yLdWnYCKx7WROLHoYeNDM1uONqb8UmUiS0JyD138En73DivMfZ02NDkd9OzUlmeHZujyuSHmqdGKRc+4t4K3w80+BLpGPJAntP3/ybvB8+XCyrvoeT7TOY+T8XPILi2iZlsrw7Pb6QFTkOHSmqMSO3Lnw5s+hQy/o/lMAcjpnqMBFKknXcpHYsG0VTLsbWp4HOX+BJP3VFKkq/asR/+0r8G4hV6s+9JsANetU/BoROYaGXMRfoYMwcQDs3+ldcKtBut+JROKWCl384xzMvg+2fAB9xkLLzn4nEolrGnIR/7z9NKycAlf9HDrc5HcakbinQhd/rJ4Ji34D5/aFy37sdxqRQFChS/XL+whmfA9aXwQ3/lH3AxWJEBW6VK8v8mHi7VC3GfQdBym1/U4kEhj6UFSqT/F+b3riob0w9E2op6tvikSSCl2qR2mpN8yybSX0nwgtzvY7kUjgqNCleix6HNbOhuzfwrey/U4jEkgaQ5foWz4J/vU0nD8Iun7f7zQigaVCl+ja9J538lDmZXDD7zSjRSSKVOgSPbs2eDNa0tpAn79BcorfiUQCTYUu0VFUCOP7giuF2ydDncZ+JxIJPH0oKpFXEoIpg2HXpzBwJjQ53e9EIglBhS6R5RzMfQg+XQS9noPMbn4nEkkYGnKRyHrveVj8Mlz6AHS+w+80IglFhS6RkzsP5v8UzroRejzqdxqRhKNCl8jYugKmDYX0TnDzKN1CTsQHFf6rM7PaZvaBmS03s9Vm9svw+rZm9r6ZrTOzSWZWM/pxJSbtyYPxfaB2Q++0ft1CTsQXlTmMOgRc5ZzrBJwHXGtmXYGngN8759oBu4Gh0YspMevQXm964qF9MGCKbiEn4qMKC9159oUXU8JfDrgKmBpePwbIiUpCiV0lh73piQVroM8YXXBLxGeVGug0s2QzWwYUAAuAT4BC59zh8CZbgIzjvHaYmS02s8U7duyIRGaJBc7B3OGw/u/Q8xk4o4ffiUQSXqUK3TlX4pw7D2gFdAHOKm+z47x2lHMuyzmX1ayZrn8dGP/+v6+mJ14w2O80IkIVZ7k45wqBt4CuQJqZHTkxqRWQH9loErPWzIIFP4ezb9b0RJEYUplZLs3MLC38PBW4GlgLLAJuC282CJgVrZASQzZ/ANOHQasukPO8pieKxJDKnPqfDowxs2S8/wAmO+fmmNkaYKKZ/QZYCrwUxZwSC3au92a01E+H/hMgJdXvRCJSRoWF7pxbAXQuZ/2neOPpkgj2bodXbwFLgjunQ92mficSka/RxbmkYof2wvjesH8HDJ4DjU/zO5GIlEOFLt+sJASTB8G2Vd4wS8YFficSkeNQocvxOQezfwifLISb/k83dxaJcZqiIMe36HFYPh66j4DzB/qdRkQqoEKX8n34Irw90ivyKx72O42IVIIKXY61ahq8/hNolw03/B7M/E4kIpWgQpejrf87TP8utOkKvUdDsj5mEYkXKnT5yuYPYdKd0OxMXddcJA6p0MVTsBbG3Qb1WsAd0yA1ze9EIlJFKnSB3Z/B2JuhRm0YOBPqt/A7kYicAA2QJrp9BTA2B0IH4K650CjT70QicoJU6InswC4Yewt8sRUGztIdh0TinAo9UR3cA6/eCjtzvVP621zkdyIROUkq9ER0aB+M6wPbVkCfsXDG1X4nEpEI0IeiiSZUxI4XbqFk8wfce/D7XDqzNjOX5vmdSkQiQEfoieTwIba/cBvNdnzAg6F7eL20KxQWMWL6SgByOpd7n28RiRM6Qk8UJSGYchctCt5hxOG7mVna7ctvFYVKGDk/18dwIhIJKvREUHLYuw9o7uv8IjSYSSVXHrNJfmGRD8FEJJJU6EFXEoJpQ2D1dLjm1yys36vczVqm6f6gIvFOhR5kh4thymBYMwuyfwuX/pDh2e1JTUk+arPUlGSGZ7f3J6OIREyFhW5mrc1skZmtNbPVZnZ/eH1jM1tgZuvCj42iH1cq7fAhmDwQPp4D1/0vXHwv4H3w+cQtHclIS8WAjLRUnriloz4QFQkAc8598wZm6UC6c+4jM6sPLAFygMHALufck2b2CNDIOfeNd0LIyspyixcvjkxyOb7QQZh8J6x7E274HVx4t9+JROQkmNkS51xWRdtVeITunNvqnPso/HwvsBbIAHoBY8KbjcErefFbqAgm9od1C+DGZ1XmIgmkSvPQzSwT6Ay8D7Rwzm0Fr/TNrHnE00nVHNrnlfmGf0Gv56DzAL8TiUg1qvSHomZWD5gGPOCc+6IKrxtmZovNbPGOHTtOJKNUxoFd8LdesPEduPkvKnORBFSpI3QzS8Er83HOuenh1dvNLD18dJ4OFJT3WufcKGAUeGPoEcicUGYuzWPk/FzyC4tomZbK8Oz2x36A+UW+dz3zXRu8a7Oc1dOfsCLiq8rMcjHgJWCtc+6ZMt+aDQwKPx8EzIp8vMQ2c2keI6avJK+wCAfkhU/TP+raKzvXw0vZsCfPu9OQylwkYVVmyOVS4E7gKjNbFv66HngSuMbM1gHXhJclgkbOz6UoVHLUuqNO09+6HF7O9m5OMXgOtL3Mh5QiEisqHHJxzr0D2HG+3SOycaSs452On19Y5H3wOaG/d+/PO2dC0zOqOZ2IxBqdKRrDjnc6fv/6y7ybUzRoCUPmq8xFBFChx7RjT9N3fL/mGzweGgnp58KQedBQZ3iKiEfXQ49hR2azjJyfy/bCfTxddyw5JW9ChxxvamKKLqglIl9Roce4nM4Z5JxVH6beBev/Dt1+BFf9ApL0y5WIHE2FHuv2bIHxfaFgrXcq/wWD/U4kIjFKhR7L8pfBhH5QvB/umAqnX+V3IhGJYfq9PVYtn+TNMU+q4c1kUZmLSAV0hB5rSkKw4Bfw3p/h1G7QezTUa+Z3KhGJAyr0WLJ/p3eHoY3/govugW//GpJT/E4lInFChR4r8pfCxDvgwE64+a/QqZ/fiUQkzqjQY8Gy8fDaA1CvuTde3vI8vxOJSBxSofvp0D544yewfAJkXuaNl9dt6ncqEYlTKnS/bF0OU+6C3Rug+wi4fDgkJVf8OhGR41ChVzfn4INR8ObPoE4TGPQaZHbzO5WIBIAKvTod2AWz7oPc16FdNuQ8D3Wb+J1KRAJChV5d1i+E2T+AfQWQ/QR0vQfseJeZFxGpOhV6tB3a6w2vLBkNTb8FQ9+EjPP9TiUiAaRCj6ZP/+kNsezZDJf8AK78H13yVkSiRoUeDcX7YcGj8OEL0Ph0b255m4v8TiUiAadCj6CZS/NY9MZkfnToeTKTtvPJaXdyer//hZp1/I4mIglAV1uMkLn/WUbKjLt5NvQYDuh76Of0XNeTmat3+x1NRBJEhYVuZi+bWYGZrSqzrrGZLTCzdeHHRtGNGcNKS+D9v3LZ/Ou42j7k96Fbubb4Kd53Z1EUKmHk/Fy/E4pIgqjMEfpo4NqvrXsEWOicawcsDC8nnrwl8MKVMPchPio5neziJ3m25FYOUfPLTfILi3wMKCKJpMJCd869Dez62upewJjw8zFAToRzxbYvtnqzV17oAXu3w22vMKLOL9no0o/ZtGWaZrWISPU40Q9FWzjntgI457aaWfMIZopdxfvh3T/Cv//o3Yji4nvhioehdgOGh/IYMX0lRaGSLzdPTUlmeHZ7HwOLSCKJ+iwXMxsGDANo06ZNtH/cSZm5NI+R83PJLyyiZVoqw7Pbk9M5wxsnXzYO/vE47NsGZ98MPR6Fxm2/fG1O5wyA8l8vIlINzDlX8UZmmcAc59w54eVcoHv46DwdeMs5V+GhaFZWllu8ePHJJY6SmUvLO8JO4uVLd3Pxp3+CgtXQqgtkPw6tu/iYVEQSjZktcc5lVbTdiR6hzwYGAU+GH2ed4J8TM0bOzy1T5o7uSct5wKZx3nufQKNM71rlHXJ0/RURiVkVFrqZTQC6A03NbAvwKF6RTzazocAmoHc0Q1YHbzaK44qkFTxQYxqdk9azxTXlkdB3ePLe30KNmhX+GSIifqqw0J1z/Y/zrR4RzuIf57il/hruODSxTJHfzbSSy2meVl9lLiJxIbFP/T98CFZOgX//id+F1pJvTRkRGsrUkisIUUOzVEQkriRmoR/YBYtf9u4ctG87tDgHcv7ChyVdeXvBBg4XFpGhWSoiEmcSq9AL1sKHL3lTEEMH4PQecPNf4bTuYEYvoNcFbSv4Q0REYlPwCz10ENbMgiWvwKb/QHJN6NjbOymoxdl+pxMRiZjgFvrO9V6JLxsHRbuh8Wlwza/hvAG6j6eIBFKwCv3ALlg9HVZMhs3vQ1INOPMGyBoCmZdDkq4WLCLBFf+FHjoI/50HKybBugVQGoJmZ8HVj0Gn26F+C78TiohUi/gs9FARfLII1s6Gj9+AQ3ug3ilw0Xfh3L5wSked0SkiCSd+Cv3gF7DuTVj7mnckHtoPtRvCmdfDuX2g7RWQlOx3ShER38RHob/2gPfhZkkx1G3uFXiHmyDzMkhO8TudiEhMiI9CT2sDF34HzrrRu9KhjsRFRI4RH4V+2YN+JxARiXmaxyciEhAqdBGRgFChi4gEhApdRCQgVOgiIgGhQhcRCQgVuohIQKjQRUQCwpxz1ffDzHYAn53gy5sCOyMYJx7oPScGvefgO9n3e6pzrllFG1VroZ8MM1vsnMvyO0d10ntODHrPwVdd71dDLiIiAaFCFxEJiHgq9FF+B/CB3nNi0HsOvmp5v3Ezhi4iIt8sno7QRUTkG8RFoZvZtWaWa2brzewRv/NEmpm1NrNFZrbWzFab2f3h9Y3NbIGZrQs/NvI7a6SZWbKZLTWzOeHltmb2fvg9TzKzmn5njCQzSzOzqWb2cXh/Xxz0/WxmPwr/vV5lZhPMrHbQ9rOZvWxmBWa2qsy6cveref4Y7rMVZnZ+pHLEfKGbWTLwHHAd0AHob2Yd/E0VcYeBHzvnzgK6AveG3+MjwELnXDtgYXg5aO4H1pZZfgr4ffg97waG+pIqep4F5jnnzgQ64b33wO5nM8sAfghkOefOAZKBfgRvP48Grv3auuPt1+uAduGvYcDzkQoR84UOdAHWO+c+dc4VAxOBXj5niijn3Fbn3Efh53vx/pFn4L3PMeHNxgA5/iSMDjNrBdwAvBheNuAqYGp4k0C9ZzNrAFwOvATgnCt2zhUS8P2Md2e0VDOrAdQBthKw/eycexvY9bXVx9uvvYC/Oc97QJqZpUciRzwUegawuczylvC6QDKzTKAz8D7Qwjm3FbzSB5r7lywq/gA8BJSGl5sAhc65w+HloO3r04AdwCvhYaYXzawuAd7Pzrk84GlgE16R7wGWEOz9fMTx9mvUOi0eCt3KWRfIqTlmVg+YBjzgnPvC7zzRZGY9gQLn3JKyq8vZNEj7ugZwPvC8c64zsJ8ADa+UJzxu3AtoC7QE6uINOXxdkPZzRaL29zweCn0L0LrMcisg36csUWNmKXhlPs45Nz28evuRX8XCjwV+5YuCS4GbzGwj3jDaVXhH7GnhX80hePt6C7DFOfd+eHkqXsEHeT9fDWxwzu1wzoWA6cAlBHs/H3G8/Rq1TouHQv8QaBf+VLwm3gcqs33OFFHhseOXgLXOuWfKfGs2MCj8fBAwq7qzRYtzboRzrpVzLhNvn/7DOTcAWATcFt4saO95G7DZzNqHV/UA1hDg/YyV5xXJAAAAy0lEQVQ31NLVzOqE/54fec+B3c9lHG+/zgYGhme7dAX2HBmaOWnOuZj/Aq4H/gt8AvyP33mi8P664f3KtQJYFv66Hm9MeSGwLvzY2O+sUXr/3YE54eenAR8A64EpQC2/80X4vZ4HLA7v65lAo6DvZ+CXwMfAKmAsUCto+xmYgPcZQQjvCHzo8fYr3pDLc+E+W4k3AygiOXSmqIhIQMTDkIuIiFSCCl1EJCBU6CIiAaFCFxEJCBW6iEhAqNBFRAJChS4iEhAqdBGRgPh/qP21z2dL9a0AAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(data['x'], data[\"y\"], 'o')\n",
"x = np.linspace(0, 100)\n",
"plt.plot(x, model(x, *popt));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, this is where `MultiNorm` comes in!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## MultiNorm\n",
"\n",
"In the last section we fit our model `y = m x + b` with two parameters `(m, b)` to some data, and optained the best-fit parameter vector `popt` and a covariance matrix `pcov`. Using this information, we can create a `MultiNorm` object.\n",
"\n",
"### Create"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"from multinorm import MultiNorm"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"mn = MultiNorm(mean=popt, cov=pcov, names=[\"b1\", \"b2\", \"b3\"])"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"MultiNorm with n=3 parameters:\n",
" mean error\n",
"name \n",
"b1 72.462232 1.734039\n",
"b2 2.618077 0.088295\n",
"b3 0.067359 0.003447"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mn"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`MultiNorm` is very simple, it just stores the `mean` vector and `cov` matrix\n",
"and exposes them as pandas Series and DataFrame objects, with `name` as index.\n",
"\n",
"Of course the `mean` and `cov` are still available as Numpy arrays, using the `.values` property of the pandas objects."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"['b1', 'b2', 'b3']"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mn.names"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"name\n",
"b1 72.462232\n",
"b2 2.618077\n",
"b3 0.067359\n",
"Name: mean, dtype: float64"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mn.mean"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
name
\n",
"
b1
\n",
"
b2
\n",
"
b3
\n",
"
\n",
"
\n",
"
name
\n",
"
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
b1
\n",
"
3.006892
\n",
"
-0.069747
\n",
"
-0.005014
\n",
"
\n",
"
\n",
"
b2
\n",
"
-0.069747
\n",
"
0.007796
\n",
"
0.000250
\n",
"
\n",
"
\n",
"
b3
\n",
"
-0.005014
\n",
"
0.000250
\n",
"
0.000012
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"name b1 b2 b3\n",
"name \n",
"b1 3.006892 -0.069747 -0.005014\n",
"b2 -0.069747 0.007796 0.000250\n",
"b3 -0.005014 0.000250 0.000012"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mn.cov"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Errors\n",
"\n",
"In the printout of `mn` above, we saw that in addition to the `names`, `mean` and `cov` we passed in, another property `error` exists.\n",
"\n",
"This is the vector of \"1 sigma\" parameter error estimates:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"name\n",
"b1 1.734039\n",
"b2 0.088295\n",
"b3 0.003447\n",
"Name: error, dtype: float64"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mn.error"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Usually you would state your measurement giving the estimated value and error for each parameter:\n",
"- `b1 = 72.46 +/- 1.73`\n",
"- `b2 = 2.618 +/- 0.088`\n",
"- `b3 = 0.0673 +/- 0.0034`\n",
"\n",
"The relationship between the covariance matrix and the errors is very simple: the diagonal entries of the covariance matrix contains the errors squared, meaning that you can compute the errors as the square root of the diagonal entries:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1.73403926, 0.08829495, 0.00344656])"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.sqrt(pcov.diagonal())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Parameters\n",
"\n",
"Possibly the property you'll use most is `parameters`, which is a `pandas.DataFrame` with `name` as index, and columns that contain `mean` and `error` for all parameters. It has a nice table printout, and you can use the common pandas methods to select rows and columns and entries."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"