{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Curve-Fitting comparison: Python and Mathematica\n",
"\n",
"+ Implementation of curve-fitting in Python.
\n",
"+ Starting point: version of mlr.py written for PHYS 310 Jack Gallimore.
\n",
"+ Comments based on treatment in Numerical Recipes (Chapter 14.3 in my copy).
\n",
"+ Compare with results of Mathematica for same data sets: see pythonTest.nb.\n",
"\n",
"Define function that is a linear combination of \"Basis Functions,\" $X_k(x)$, evaluated at all values of independent variable $x_j$ -- see Eq. (14.3.2) in Numerical Recipes. These are combined into the single matrix that gets weighted for the illustrated matrix in Fig. 14.3.1 of Numerical Recipes.\n",
"\n",
"Covariance Matrix; See Eq. (14.3.8), and following discussion:\n",
"\n",
"$$\\alpha_{kj} = \\sum_{i=1}^N\\frac{X_j(x_i)X_k(x_i)}{\\sigma_i^2} \\hspace{0.2in}\\mbox{and}\n",
" \\hspace{0.2in}C_{jk} = [\\alpha]_{jk}^{-1}$$\n",
"\n",
"This result gives covariance matrix if uncertainties are measurement uncertainties with well-known values. To get covariance matrix from uncertainties estimated from observed spread of values, multiply $C_{jk}$ by $\\chi_R^2$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Linear fitting "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy import linalg\n",
"from scipy import stats\n",
"\n",
"import matplotlib as mpl # As of July 2017 Bucknell computers use v. 2.x \n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Following is an Ipython magic command that puts figures in the notebook.\n",
"# For figures in separate windows, comment out following line and uncomment\n",
"# the next line\n",
"# Must come before defaults are changed.\n",
"%matplotlib notebook\n",
"#%matplotlib\n",
"\n",
"# As of Aug. 2017 reverting to 1.x defaults.\n",
"# In 2.x text.ustex requires dvipng, texlive-latex-extra, and texlive-fonts-recommended, \n",
"# which don't seem to be universal\n",
"# See https://stackoverflow.com/questions/38906356/error-running-matplotlib-in-latex-type1cm?\n",
"mpl.style.use('classic')\n",
" \n",
"# M.L. modifications of matplotlib defaults using syntax of v.2.0 \n",
"# More info at http://matplotlib.org/2.0.0/users/deflt_style_changes.html\n",
"# Changes can also be put in matplotlibrc file, or effected using mpl.rcParams[]\n",
"plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures\n",
"plt.rc('axes', labelsize=16, titlesize=14)\n",
"plt.rc('figure', autolayout = True) # Adjusts supblot parameters for new size"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### THE FOLLOWING CELL DEFINING BASIS FUNCTIONS IS ONLY FUNCTION CELL THAT SHOULD REQUIRE MODIFICATION FOR NEW MODEL"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Basis functions for linear model: func = a0*X0 + a1*X1 + a2*X2 + ...\n",
"def basis(x):\n",
" ''' Returns basis functions for linear model\n",
" \n",
" The function to be fit is assumed to be of the form\n",
" f(x) = a0*X0 + a1*X1 + a2*X2 + a3*X3 + ...\n",
" where a0, a1, a2, ... are constants, and X0, X1, X2, ... are defined below.\n",
" '''\n",
" X2 = x**2 \n",
" X1 = x\n",
" X0 = 0.*X1 + 1. # Need array of len(x), thus the 0.*X1\n",
" return np.array([X0,X1,X2])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Evaluate function that is a linear combination of basis functions, with cofficients contained in array `a`"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def func(x, a):\n",
" '''Given basis functions and coefficients, returns value of linear function'''\n",
" return np.dot(basis(x).T, a)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Function that does fitting, and formats output:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def LinearModelFit(x, y, u):\n",
" '''Performs linear least squares fit to a set of 2-d data with uncertainties\n",
" \n",
" x = array of x values [x0, x1, x2, ...]\n",
" y = array of values of dependent variable [y0, y1, y2, ...]\n",
" u = array of uncertainties for dependent variable [u0, u1, u2, ...]\n",
" '''\n",
" X = basis(x).T # Basis functions evaluated at all x (the X_j(x_i)) of N.R.)\n",
" W = np.diag(1/u) # Matrix with uncertainties on diagonal\n",
" Xw = np.dot(W,X) # A_ij of Eq. (14.3.4)\n",
" Yw = np.dot(y,W) # b_i of Eq. (14.3.5)\n",
" fit = np.linalg.lstsq(Xw,Yw,rcond=1) # lstq returns: best values, chi2, ....\n",
" covariance = np.linalg.inv(np.dot(Xw.T,Xw))\n",
" uncertainty = np.sqrt(np.diag(covariance))\n",
" return(fit[0], uncertainty,fit[1], covariance)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Generate some noisy data (x,y) with uncertainties u, or import sample data file."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# Cell for generating data; overwritten by following cell if data is coming from file.\n",
"x = np.linspace(0, 10, 11)\n",
"sigma = 5.0\n",
"a = np.array([1.2, 3.4, -0.9])\n",
"y = func(x,a) + stats.norm.rvs(0, sigma, size=len(x))\n",
"u = sigma*np.ones(len(y))\n",
"#sp.savetxt(\"sample.dat\",sp.array([x,y,u]).T)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"data = np.loadtxt(\"sample.dat\") #Each line in file corresponds single data point: x,y,u"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"x = data.T[0]\n",
"y = data.T[1]\n",
"u = data.T[2]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Two previous cells would be more \"pythonic\" as\n",
"\n",
" x, y, u = sp.loadtxt(\"sample.dat\", unpack=True)\n",
"\n",
"The \"unpack = True\" reads columns."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"x, y, u = np.loadtxt(\"sample.dat\", unpack=True)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"application/javascript": [
"/* Put everything inside the global mpl namespace */\n",
"window.mpl = {};\n",
"\n",
"\n",
"mpl.get_websocket_type = function() {\n",
" if (typeof(WebSocket) !== 'undefined') {\n",
" return WebSocket;\n",
" } else if (typeof(MozWebSocket) !== 'undefined') {\n",
" return MozWebSocket;\n",
" } else {\n",
" alert('Your browser does not have WebSocket support.' +\n",
" 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n",
" 'Firefox 4 and 5 are also supported but you ' +\n",
" 'have to enable WebSockets in about:config.');\n",
" };\n",
"}\n",
"\n",
"mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n",
" this.id = figure_id;\n",
"\n",
" this.ws = websocket;\n",
"\n",
" this.supports_binary = (this.ws.binaryType != undefined);\n",
"\n",
" if (!this.supports_binary) {\n",
" var warnings = document.getElementById(\"mpl-warnings\");\n",
" if (warnings) {\n",
" warnings.style.display = 'block';\n",
" warnings.textContent = (\n",
" \"This browser does not support binary websocket messages. \" +\n",
" \"Performance may be slow.\");\n",
" }\n",
" }\n",
"\n",
" this.imageObj = new Image();\n",
"\n",
" this.context = undefined;\n",
" this.message = undefined;\n",
" this.canvas = undefined;\n",
" this.rubberband_canvas = undefined;\n",
" this.rubberband_context = undefined;\n",
" this.format_dropdown = undefined;\n",
"\n",
" this.image_mode = 'full';\n",
"\n",
" this.root = $('