{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Week 6 Day 3: Fitting\n", "\n", "## Objectives\n", "\n", "* Learn how to interpolate using several methods\n", "* Learn how to perform a simple fit on data\n", "* Learn about various tools for data fitting" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's prepare a set of \"semi-continuous\" values to interpolate on." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = np.linspace(0, 200, 100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And, here's the data for a cross-section of resonant scattering of a neutron from a nucleus." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x_data = np.linspace(0, 200, 9)\n", "y_data = np.array([10.6, 16.0, 45.0, 83.5, 52.8, 19.9, 10.8, 8.25, 4.7])\n", "# e_data = np.array([9.34, 17.9, 41.5, 85.5, 51.5,\n", "# 21.5, 10.8, 6.29, 4.14])\n", "e_data = np.array([5, 7, 10, 12, 10, 7, 5, 4, 4])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.errorbar(x_data, y_data, e_data, fmt=\"o\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Linear interpolation\n", "\n", "We could numerically interpolate between values. The numpy function `interp` does this for us:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "y = np.interp(x, x_data, y_data)\n", "fig, ax = plt.subplots()\n", "ax.errorbar(x_data, y_data, e_data, fmt=\"o\")\n", "ax.plot(x, y, \".\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Lagrange interpolation\n", "\n", "We can build our own Lagrange interpolation function:\n", "\n", "$$\n", "g(x)\n", "=\n", "\\sum_{i=1}^{n}\n", "g_i \\lambda_i(x)\n", "$$\n", "\n", "$$\n", "\\lambda_i(x)\n", "=\n", "\\prod_{j(\\ne i)=1}^{n}\n", "\\frac{x - x_n}\n", " {x_i - x_n}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def interp_langrange(x_data, y_data, x):\n", " \"A custom Lagrange interpolate function.\"\n", "\n", " # The shape will be N x N x n, where N is the number of points\n", " # in x and y, and n is the number of points in the result.\n", " x = np.asarray(x)\n", " x = x.reshape(1, 1, -1)\n", "\n", " # Make grid. Have one \"throw away\" diminsion to make [N,N,1] shaped arrays\n", " x_n, x_i, _ = np.meshgrid(x_data, x_data, np.array([1]))\n", "\n", " # Using a mask here avoids numpy warnings and makes the product easy\n", " x_n = np.ma.array(x_n, mask=np.eye(x_data.shape[0]))\n", " x_i = np.ma.array(x_i, mask=np.eye(x_data.shape[0]))\n", "\n", " V = (x - x_n) / (x_i - x_n)\n", " v = np.prod(V, axis=1).data # Convert back to non-masked arrays\n", "\n", " return y_data @ v\n", "\n", "\n", "# This was added to test interp_langrange\n", "x_d2 = np.array([0, 1, 2, 4])\n", "y_d2 = np.array([-12, -12, -24, -60])\n", "\n", "interp_langrange(x_d2, y_d2, [0.5, 4]) # should be -10.125 and -60" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.errorbar(x_data, y_data, e_data, fmt=\"o\")\n", "ax.plot(x, interp_langrange(x_data, y_data, x), \".\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Of course, you shouldn't be writing your own algorithms for something that is generally useful - it's included in SciPy. Let's see what that would have looked like instead:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy import interpolate" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "p = interpolate.lagrange(x_data, y_data) # Makes a callable object \"p\"" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.errorbar(x_data, y_data, e_data, fmt=\"o\")\n", "ax.plot(x, p(x), \".\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cubic splines\n", "\n", "This is a visually appealing method for interpolation. It makes small cubic polynomial fits, with the further requirement that they smoothly connect to the next segment (thus looks better than the similar integration tool we saw earlier)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "spline = interpolate.CubicSpline(x_data, y_data)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.plot(x, spline(x))\n", "ax.plot(x_data, y_data, \"o\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Curve fitting\n", "\n", "We'll start with using the scipy optimize package. We'll look at how to do this ourselves with a simple function later." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import curve_fit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our fit function is a theoretical model: a Breit-Wigner resonance:\n", "\n", "$$\n", "f(E) =\n", "\\frac{f_r}{\\left(E-E_r\\right)^2 + \\Gamma^2/4}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def f(E, f_r, E_r, Γ):\n", " return f_r / ((E - E_r) ** 2 + Γ ** 2 / 4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we look up the definition of curve_fit (using shift-tab, `curve_fit?`, or `help(curve_fit)`), and then feed it the parameters it wants:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "popt, pcov = curve_fit(\n", " f, x_data, y_data, p0=[1.0, 1.0, 1.0], sigma=e_data, absolute_sigma=True\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It returns two objects; the array of fit value, and a covariance matrix." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.errorbar(x_data, y_data, e_data, fmt=\".\")\n", "ax.plot(x, f(x, *popt))\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Advanced usage: The Uncertainties package\n", "\n", "Let's use that correlation matrix that `curve_fit` returned! If you have the uncertainties package, you can make uncertain values - even ones related through an covariance matrix." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from uncertainties import correlated_values, unumpy" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "copt = correlated_values(popt, pcov)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, since we can compute the value + uncertainty for a calculation, including the correct relationship between variables, we can compute the function again, and plot a shaded uncertainty band." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_with_unc(x, y, ax=None):\n", " if ax is None:\n", " fig, ax = plt.subplots()\n", "\n", " y_nv = unumpy.nominal_values(y)\n", " y_sd = unumpy.std_devs(y)\n", "\n", " (line,) = ax.plot(x, y_nv)\n", " ax.fill_between(x, y_nv - y_sd, y_nv + y_sd, color=line.get_color(), alpha=0.2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.errorbar(x_data, y_data, e_data, fmt=\".\")\n", "plot_with_unc(x, f(x, *copt), ax=ax)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Linear least squares\n", "\n", "Let's do a least squares fit ourselves. As long as we have a linear dependence *on the parameters*, we can solve this with linear algebra.\n", "\n", "Some data:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = np.linspace(1, 2)\n", "x = np.array([1.0, 1.1, 1.24, 1.35, 1.451, 1.5, 1.92])\n", "y = np.array([0.52, 0.8, 0.7, 1.8, 2.9, 2.9, 3.6])\n", "sig = np.array([0.1, 0.1, 0.2, 0.3, 0.2, 0.1, 0.1])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.errorbar(x, y, sig, fmt=\".\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Manual least squares fitting:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fitting function:\n", "\n", "$$\n", "g(x_i) = a_0 x_i^2 + a_1 x_i + a_2\n", "\\tag{1}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you compute the minimum of $\\chi^2$, you get:\n", "\n", "$$\n", "\\sum _i\n", "\\frac{g(x_i)}{\\sigma_i^2}\n", "\\frac{\\partial g ( x_i)}{\\partial a_m}\n", "=\n", "\\sum _i\n", "\\frac{y_i}{\\sigma_i^2}\n", "\\frac{\\partial g ( x_i)}{\\partial a_m}\n", "\\tag{2}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can rewrite this as a matrix expression, expanding (2) in terms of (1):\n", "\n", "$$\n", "A x = b\n", "\\tag{3}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "A = \\left[\n", "\\begin{matrix}\n", " \\sum_i \\frac{x_i^4}{\\sigma_i^2} & \\sum_i \\frac{x_i^3}{\\sigma_i^2} & \\sum_i \\frac{x_i^2}{\\sigma_i^2} \\\\\n", " \\sum_i \\frac{x_i^3}{\\sigma_i^2} & \\sum_i \\frac{x_i^2}{\\sigma_i^2} & \\sum_i \\frac{x_i^1}{\\sigma_i^2} \\\\\n", " \\sum_i \\frac{x_i^2}{\\sigma_i^2} & \\sum_i \\frac{x_i }{\\sigma_i^2} & \\sum_i \\frac{1 }{\\sigma_i^2} \\\\\n", "\\end{matrix}\n", "\\right]\n", "\\tag{3.1}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "b = \\left[\n", "\\begin{matrix}\n", " \\sum_i \\frac{y_i x_i^2}{\\sigma_i^2} \\\\\n", " \\sum_i \\frac{y_i x_i }{\\sigma_i^2} \\\\\n", " \\sum_i \\frac{y_i }{\\sigma_i^2} \\\\\n", "\\end{matrix}\n", "\\right]\n", "\\tag{3.2}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x = \\left[\n", "\\begin{matrix}\n", " a_2 \\\\\n", " a_1 \\\\\n", " a_0 \\\\\n", "\\end{matrix}\n", "\\right]\n", "\\tag{3.3}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sig2 = sig ** 2\n", "ss = np.array(\n", " [\n", " np.sum(x ** 4 / sig2),\n", " np.sum(x ** 3 / sig2),\n", " np.sum(x ** 2 / sig2),\n", " np.sum(x / sig2),\n", " np.sum(1 / sig2),\n", " ]\n", ")\n", "A = np.stack([ss[:3], ss[1:4], ss[2:]])\n", "\n", "b = np.array([np.sum(x ** 2 * y / sig2), np.sum(x * y / sig2), np.sum(y / sig2)])\n", "print(b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could have been much more clever. Looking at the above expression for M, we could instead have created the powers, then raise to a matrix of powers!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "power = sum(np.mgrid[2:-1:-1, 2:-1:-1])\n", "Ap = np.sum(x[:, None, None] ** power / sig[:, None, None] ** 2, axis=0)\n", "assert np.all(A == Ap)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Exercise for the reader: try the same trick on `b`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xvec = np.linalg.solve(A, b)\n", "print(xvec)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "p = np.poly1d(xvec)\n", "fig, ax = plt.subplots()\n", "ax.errorbar(x, y, sig, fmt=\".\")\n", "ax.plot(t, p(t))\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Numpy Polyfit\n", "\n", "I'm still not advocating that real fits be done this way. We could have just done this with the polyfit function in numpy:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xvec = np.polyfit(x, y, 2, w=1 / sig)\n", "p = np.poly1d(xvec)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.errorbar(x, y, sig, fmt=\".\")\n", "ax.plot(t, p(t))\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Scipy curve fit\n", "\n", "Or the more general scipy `curve_fit`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def new_p(x, *values):\n", " return np.poly1d(values)(x)\n", "\n", "\n", "xvec, _ = curve_fit(new_p, x, y, p0=[1, 1, 1], sigma=sig)\n", "p = np.poly1d(xvec)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.errorbar(x, y, sig, fmt=\".\")\n", "ax.plot(t, p(t))\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A 2D fit\n", "\n", "Let's try in 2D. Try to fit the following data with a function of the form:\n", "\n", "$$\n", "f(x,y) =\n", "a \\sin(b x) + c \\cos(d y)\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "state = np.random.RandomState(42)\n", "x, y = np.ogrid[-np.pi : np.pi : 100j, -np.pi : np.pi : 100j]\n", "z = 1.2 * np.sin(2.3 * x) + 0.6 * np.cos(3.1 * y) + state.rand(100, 100) * 0.3" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.imshow(z)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To get this into `curve_fit`, we have to jump through some hoops. We need to give x and y together (it will be used by the function exactly as we give it), and we need to return a 1D array of values to compare:\n", "\n", "\n", "\n", "> Note: some older examples will show you a really odd syntax:\n", ">\n", "> ```python\n", "> def f((x,y), a, b, c, d):\n", "> ```\n", ">\n", "> This was more trouble than it was worth and was dropped in Python 3. Just unpack yourself if you want to unpack.\n", "\n", "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def f(xy, a, b, c, d):\n", " x, y = xy\n", " z = a * np.sin(b * x) + c * np.cos(d * y)\n", " return z.ravel()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We had to return a flattened array, since `curve_fit` only likes flat (1D) arrays. At this point, when computing a least squared fit, the result shape does not matter, so it's okay to do. (The fact we had to build it into our fitting function is really irritating, though). We'll need to do the same thing with `z`, as well, since it's going to be compared with the return value:\n", "\n", "> Note: `flatten()`, `ravel()`, and `reshape(-1)` all do basically the same thing, but `flatten` always makes a copy, so `ravel` is usually better." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll need to broadcast our 100x1 and 1x100 arrays to 100x100 so we can flatten them:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xfull, yfull = np.broadcast_arrays(x, y)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "popt, pcov = curve_fit(f, [xfull.ravel(), yfull.ravel()], z.ravel(), p0=[1, 2, 0.5, 3])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's use the uncertainties package again, this time to provide nicer printouts of the fit values:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "copt = correlated_values(popt, pcov)\n", "for i in range(4):\n", " print(\"abcd\"[i], \"=\", format(copt[i], \"P\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Just for comparison, this is what it would have looked like if we calculated that ourselves:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for i in range(4):\n", " print(\"abcd\"[i], \"=\", popt[i], \"±\", np.sqrt(pcov[i, i]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, the result:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.imshow(f([x, y], *popt).reshape(100, 100))\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Try it yourself:\n", "\n", "* What happens if your initial guess is not as good? Why?" ] } ], "metadata": { "kernelspec": { "display_name": "compclass", "language": "python", "name": "compclass" }, "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.8.2" } }, "nbformat": 4, "nbformat_minor": 4 }