8S50x_P04.ipynb

{"cells": [{"cell_type": "markdown", "id": "b6cb2f21", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<hr style=\"height: 1px;\">\n", "<i>This notebook was authored by the 8.S50x Course Team, Copyright 2022 MIT All Rights Reserved.</i>\n", "<hr style=\"height: 1px;\">\n", "<br>\n", "\n", "<h1>Guided Problem Set 4: Fitting with LMFIT</h1>\n"]}, {"cell_type": "markdown", "id": "63166c6c", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_4_0'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">P4.0 Overview</h2>\n"]}, {"cell_type": "markdown", "id": "0bd10115", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<h3>Navigation</h3>\n", "\n", "<table style=\"width:100%\">\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_4_1\">P4.1 Using LMFIT to Fit Data</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#problems_4_1\">P4.1 Problems</a></td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_4_2\">P4.2 Another LMFIT Example</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#problems_4_2\">P4.2 Problems</a></td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_4_3\">P4.3 A More Complicated Model</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#problems_4_3\">P4.3 Problems</a></td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_4_4\">P4.4 Fitting Data Containing Noise</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#problems_4_4\">P4.4 Problems</a></td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_4_5\">P5.5 Interpreting Common Errors with LMFIT</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#problems_4_5\">P5.5 Problems</a></td>\n", "    </tr>\n", "</table>\n", "\n"]}, {"cell_type": "markdown", "id": "2bea26bf", "metadata": {"tags": ["learner", "catsoop_00", "md"]}, "source": ["<h3>Learning Objectives</h3>\n", "\n", "In these problems, we will explore the following objectives:\n", "\n", "- Fitting with the software package `lmfit`, including some more complicated models\n", "- How to interpret errors that can appear when using `lmfit`\n", "\n"]}, {"cell_type": "markdown", "id": "7a74d690", "metadata": {"tags": ["learner", "md"]}, "source": ["<h3>Importing Libraries</h3>\n", "\n", "Before beginning, run the cell below to import the relevant libraries for this notebook."]}, {"cell_type": "code", "execution_count": null, "id": "91d60e99", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.0-runcell00\n", "\n", "!pip install lmfit"]}, {"cell_type": "code", "execution_count": null, "id": "37643ed7", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.0-runcell01\n", "\n", "import numpy as np                    #https://numpy.org/doc/stable/\n", "np.random.seed(421421) # ensure calls to np.random give repeatable results\n", "\n", "import matplotlib.pyplot as plt       #https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.html\n", "\n", "from lmfit.models import LinearModel  #https://lmfit.github.io/lmfit-py/builtin_models.html#lmfit.models.LinearModel\n", "\n", "import scipy.stats as stats           #https://docs.scipy.org/doc/scipy/reference/stats.html\n"]}, {"cell_type": "markdown", "id": "6ca4a6cd", "metadata": {"tags": ["learner", "md"]}, "source": ["<h3>Setting Default Figure Parameters</h3>\n", "\n", "The following code cell sets default values for figure parameters."]}, {"cell_type": "code", "execution_count": null, "id": "5e447b55", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.0-runcell02\n", "\n", "#set plot resolution\n", "%config InlineBackend.figure_format = 'retina'\n", "\n", "#set default figure parameters\n", "plt.rcParams['figure.figsize'] = (9,6)\n", "\n", "medium_size = 12\n", "large_size = 15\n", "\n", "plt.rc('font', size=medium_size)          # default text sizes\n", "plt.rc('xtick', labelsize=medium_size)    # xtick labels\n", "plt.rc('ytick', labelsize=medium_size)    # ytick labels\n", "plt.rc('legend', fontsize=medium_size)    # legend\n", "plt.rc('axes', titlesize=large_size)      # axes title\n", "plt.rc('axes', labelsize=large_size)      # x and y labels\n", "plt.rc('figure', titlesize=large_size)    # figure title\n"]}, {"cell_type": "markdown", "id": "96b69ee5", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_4_1'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">P4.1 Using LMFIT to Fit Data</h2>    \n", "\n", "| [Top](#section_4_0) | [Previous Section](#section_4_0) | [Problems](#problems_4_1) | [Next Section](#section_4_2) |\n"]}, {"cell_type": "markdown", "id": "53436383", "metadata": {"tags": ["learner", "catsoop_01", "md"]}, "source": ["<h3>Overview</h3>\n", "\n", "To get some practice fitting, suppose you have some data coming from the function $y=2x$. We'll fit this data using a more general model function of $y=mx+b$"]}, {"cell_type": "markdown", "id": "564d6507", "metadata": {"tags": ["learner", "catsoop_01", "md"]}, "source": ["<h3>Generating Sample Data</h3>\n", "\n", "Let's generate some example data using made-up uncertainties which are given in a table in the code shown below. These uncertainties are assumed to be the standard deviations of normal distributions (a common assumption). Therefore, we randomly select each data point $y_i$ from a normal distribution with a mean of $2 x_i$ and a standard deviation equal to the uncertainty of point $i$.\n", "\n", "Run the following cell and print the output."]}, {"cell_type": "code", "execution_count": null, "id": "c856134f", "metadata": {"tags": ["learner", "py", "catsoop_01", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.1-runcell01\n", "\n", "import numpy as np\n", "\n", "np.random.seed(421421)\n", "\n", "xi = np.array([2,3,4,5,6,7])\n", "yi = 2*xi\n", "y_unc = np.array([0.3, 0.4, 0.45, 0.35, 0.6, 0.5]) # uncertainties of point i\n", "\n", "# randn samples \"standard normal\" dist. (mean 0, std. dev 1). \n", "yi = yi + np.random.randn(len(xi))*y_unc # Multiplying by y_unc multiplies that std. dev.\n", "\n", "print('yi with uncertainties:', yi)"]}, {"cell_type": "markdown", "id": "c461f010", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='problems_4_1'></a>     \n", "\n", "| [Top](#section_4_0) | [Restart Section](#section_4_1) | [Next Section](#section_4_2) |\n"]}, {"cell_type": "markdown", "id": "f5064a50", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 4.1.1</span>\n", "\n", "For the random sample generated above, plot the data points and their standard deviations as dots with vertical error bars. The starting code shown below already plots the expected distribution $y_{i}=2x_{i}$.\n", "\n", "How many of the data points are within 1 standard deviation of the model function (i.e., how many error bars are touching the line?). Enter your answer as a number.\n", "\n", "Note: The answer crucially depends on the assignment `np.random.seed(421421)`. If you use a different random seed, your answer may not match ours.\n", "\n", "Hint: The documentation for functions that you might find useful can be found below:\n", "\n", "- <a href=\"https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.scatter.html\" target=\"_blank\">matplotlib.pyplot.scatter</a>\n", "- <a href=\"https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.errorbar.html\" target=\"_blank\">matplotlib.pyplot.errorbar</a>\n", "- <a href=\"https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html\" target=\"_blank\">matplotlib.pyplot.plot</a>"]}, {"cell_type": "code", "execution_count": null, "id": "7c20f8bf", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P4.1.1\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "#plot the model function\n", "plt.plot(xi, 2 * xi, label='Prediction')\n", "\n", "#plot the data and the error bars\n", "#YOUR CODE HERE\n", "\n", "\n", "plt.show()"]}, {"cell_type": "markdown", "id": "9c74272c", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": [">#### Follow-up 4.1.1a (ungraded)\n", ">   \n", ">If you didn't do it already, add a legend, title, and axis labels to your graph. (Hint: keep looking at the matplotlib documentation.)\n", "\n", "<br>"]}, {"cell_type": "markdown", "id": "c0e655ad", "metadata": {"tags": ["learner", "catsoop_01", "md"]}, "source": ["<h3>Fitting a Model</h3>\n", "\n", "Now, let's try to fit these points using a model in `lmfit` to represent the data. Models can either be selected from a <a href=\"https://lmfit.github.io/lmfit-py/builtin_models.html\" target=\"_blank\">large list of functions</a> already available in `lmfit`, or you can define them yourself. Here we use a preset linear model."]}, {"cell_type": "code", "execution_count": null, "id": "2d992fb2", "metadata": {"tags": ["learner", "py", "catsoop_01", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.1-runcell02\n", "\n", "from lmfit.models import LinearModel\n", "model = LinearModel()"]}, {"cell_type": "markdown", "id": "02a98422", "metadata": {"tags": ["learner", "catsoop_01", "md"]}, "source": ["In Lesson 4, you learned about the math behind the concept of fitting. But, for now, let's use `lmfit` as a black box to do the work for us. Just keep in mind that what `lmfit` is doing behind the scenes is a minimization algorithm.\n", "\n", "Below we fit the data `(xi,yi)` with the model that we defined, using the inverse uncertainties as weights. Thus, data with larger uncertainty will have less influence on the fit.\n", "\n", "Then we print the output of the `fit_report()`, which shows the fit statistics and parameters that were computed from the fitting process."]}, {"cell_type": "code", "execution_count": null, "id": "c34ed316", "metadata": {"tags": ["learner", "py", "catsoop_01", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.1-runcell03\n", "\n", "result = model.fit(yi, x=xi, weights=1/y_unc);\n", "print(result.fit_report())"]}, {"cell_type": "markdown", "id": "e6791f3c", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 4.1.2</span>\n", "\n", "What are the values of the slope and intercept? Report your answer as a list of numbers with precision 1e-2: `[slope,intercept]`"]}, {"cell_type": "markdown", "id": "bba54533", "metadata": {"tags": ["learner", "catsoop_01", "md"]}, "source": ["<h3>Reading the Fit Report</h3>\n", "\n", "Look in the `Variables` section of the fit results which show a slope of about 2 and an intercept of about zero, consistent with the true model!\n", "\n", "Very helpfully, `lmfit` also gives you uncertainties on the fit parameters. In fact, without these uncertainties, you would have no idea whether the deviations from the expected results of a slope of 2 and an intercept of 0 were important. We will discuss how to interpret this more quantitatively later. For now, just consider that finding deviations which are smaller than the uncertainties indicates that the fit results are reasonably close to the expectation.\n", "\n", "You'll learn how `lmfit` finds these fit parameter uncertainties in a later section, as well as what the `chi-square` and `reduced chi-square` entries are."]}, {"cell_type": "markdown", "id": "21697993", "metadata": {"tags": ["learner", "catsoop_01", "md"]}, "source": ["<h3>Plot the Fit Result</h3>\n", "\n", "There's one more thing to do: It's always a good idea to verify that your fit result does actually fit your data. This is particularly true in cases where the model is more complicated than in this case. So, let's plot the data together with the fitted model."]}, {"cell_type": "code", "execution_count": null, "id": "c8de3247", "metadata": {"tags": ["learner", "py", "catsoop_01", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.1-runcell04\n", "\n", "result.plot();"]}, {"cell_type": "markdown", "id": "a9fc1c1b", "metadata": {"tags": ["learner", "catsoop_01", "md"]}, "source": ["You can see that the model fits the data very well in the bottom plot, and the top plot demonstrates that deviations (residuals) of the data from the model are random; they don't seem correlated with $x$ nor with each other. This is a good thing, because `lmfit` assumed that the data points were uncorrelated with each other when performing the fit. That the residuals appear random suggests that the linear model is indeed a good match for the data. If, for example, the data actually followed an exponential function, then we would see a trend in the residuals."]}, {"cell_type": "markdown", "id": "049331a4", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_4_2'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">P4.2 Another LMFIT Example</h2>    \n", "\n", "| [Top](#section_4_0) | [Previous Section](#section_4_1) | [Problems](#problems_4_2) | [Next Section](#section_4_3) |\n"]}, {"cell_type": "markdown", "id": "33aa0877", "metadata": {"tags": ["learner", "catsoop_02", "md"]}, "source": ["<h3>Overview</h3>\n", "\n", "Fitting using a pre-determined model is all well and good, but how do you fit to data using your own model function?"]}, {"cell_type": "markdown", "id": "cc6f196c", "metadata": {"tags": ["learner", "catsoop_02", "md"]}, "source": ["Let's arbitrarily choose the function\n", "\n", "$$f(x) = \\frac{\\cos(kx)}{x^a}$$\n", "\n", "as our model, with free parameters $k > 0$ and $a > 0$.\n", "\n", "Run the cell below to define the function. In the following problems, we will use this function to fit some data that we will generate."]}, {"cell_type": "code", "execution_count": null, "id": "1dff3a7a", "metadata": {"tags": ["learner", "py", "catsoop_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.2-runcell01\n", "\n", "import numpy as np\n", "\n", "def model_fn(x, k, a):# independent variable must be first argument\n", "    return np.cos(k * x) / x**a"]}, {"cell_type": "markdown", "id": "9c9234ba", "metadata": {"tags": ["learner", "catsoop_02", "md"]}, "source": ["In the following problems, we will generate some \"data\" by adding random noise to the output of `model_fn`, using defined `k` and `a` values. We will call this `y_data`. We will also generate some random uncertainty for each data point in `y_data`. Finally, we will use `lmfit` to fit the function `model_fn` to the data `y_data`, using the uncertainties as weights. We will see how close the fit values are to our predefined `k` and `a` values!!"]}, {"cell_type": "markdown", "id": "985a5d40", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='problems_4_2'></a>     \n", "\n", "| [Top](#section_4_0) | [Restart Section](#section_4_2) | [Next Section](#section_4_3) |\n"]}, {"cell_type": "markdown", "id": "6848d45e", "metadata": {"tags": ["learner", "catsoop_02", "md"]}, "source": ["<h3>Generating Data</h3>\n", "\n", "In these first few questions, we will focus on generating data and uncertainties from our model function `model_fn`. We will ensure that the data have characteristics constistent with \"real\" observations."]}, {"cell_type": "markdown", "id": "cfd632f2", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 4.2.1</span>\n", "\n", "First, let's generate some \"true\" values without noise. To do this, we will choose an interval to sample over, some \"true\" parameters to use, and implement the function `model_fn` to generate some \"true\" values, called `y_true`.\n", "\n", "Define an array `x` with 20 values evenly spaced on the interval $[0.1, \\pi]$. To get `y_true`, we will call `model_fn`with the inputs `x`, `TRUE_K`, and `TRUE_A`, where we use $k=\\pi$ and $a=1$ as the true values of the parameters.\n", "\n", "Print the `y_true` values."]}, {"cell_type": "code", "execution_count": null, "id": "b4876b6e", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P4.2.1\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "import numpy as np\n", "np.random.seed(2345789)\n", "\n", "x = np.linspace(0.1, np.pi, 20)\n", "TRUE_K = #YOUR CODE HERE\n", "TRUE_A = #YOUR CODE HERE\n", "y_true = #YOUR CODE HERE\n", "\n", "\n", "print(\"y_true values \",y_true)"]}, {"cell_type": "markdown", "id": "ee598164", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 4.2.2</span>\n", "\n", "Now we will define an array, `y_data`, which will be equal to the \"true\" value of the fuction `model_fn`, plus some random error (or noise). In short, we will draw the random error for each point from a Gaussian distribution, whose width is related to the uncertainty in our measurement. Thus, we must first generate some values for the uncertainty in our data.\n", "\n", "First, let's define the uncertainty `y_unc`. In real data, there are ways of estimating the uncertainty of any given measurement, but here we're making up uncertainties arbitrarily. For this case, we will sample the uncertainties from a uniform distribution in the range $[0.1, 0.5]$ with `np.random.random()`. Why not?\n", "\n", "Next, we make the assumption that the random error for each point is drawn from a Gaussian distribution with a width (stdev) that depends on the uncertainty that we have defined. If we were measuring real data, we would be estimating the uncertainty in each measurement based on our observations (here we're kind of doing things backwards).\n", "\n", "Finally, to get the random error for our `y_data`, at each point we will sample a Gaussian distribution whose mean is equal to the corresponding data point in `y_true`, and whose width that is equal to the corresponding uncertainy in `y_unc`. **Use `np.random.randn()` to do this, and note that `sigma * np.random.randn(...) + mu` will produce samples from a Gaussian with `sigma` and `mu`.**\n", "\n", "Print the `y_data` values."]}, {"cell_type": "code", "execution_count": null, "id": "af96ed79", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P4.2.2\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "import numpy as np\n", "np.random.seed(2345789)\n", "\n", "y_unc = 0.1 + 0.4 * np.random.random(len(y_true))\n", "\n", "y_data = #YOUR CODE HERE\n", "\n", "\n", "print(\"y_data values \",y_data)"]}, {"cell_type": "markdown", "id": "f4712971", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": [">#### Follow-up 4.2.2a (ungraded)\n", ">   \n", ">Print out the array of sampled uncertainties, `y_unc`. Do these values make sense given how we generated them? (Think about the average value of the uniform distribution.)\n", "\n", "\n", "<br>"]}, {"cell_type": "markdown", "id": "4ead67b5", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 4.2.3</span>\n", "\n", "To really convince ourselves that our application of the uncertainties to the data is valid, we can plot the deviation of the `y_data` values from the `y_true` and then divide by the sampled uncertainty. This gives us a formula $(y_{data}-y_{true})/\\sigma_{data}$.\n", "\n", "A histogram of this distribution is called the pull distribution. The nice thing about pull distributions is that, for correct uncertainties, it will match exactly to a normal distribution with mean $\\mu=0$ and standard deviation $\\sigma=1.0$. \n", "\n", "Fill in the code below and make the plot. What are the mean and stdev of `y_pull`? Report your answer with precision 1e-3, in the following form: `[mean(y_pull), stdev(y_pull)]`.\n"]}, {"cell_type": "code", "execution_count": null, "id": "9dafe65e", "metadata": {"tags": ["py", "draft", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P4.2.3\n", "\n", "import scipy.stats as stats\n", "#Use your y_true, y_data, and y_unc values from before!\n", "\n", "#optionally plot the uncertainty\n", "#bins=np.arange(0,0.5,0.05)\n", "#plt.hist(y_unc,bins=bins)\n", "#plt.show()\n", "\n", "\n", "y_pull = #YOUR CODE HERE\n", "y_pull_mean = #YOUR CODE HERE\n", "y_pull_stdev = #YOUR CODE HERE\n", "\n", "print('mean:',y_pull_mean)\n", "print('stdev:',y_pull_stdev)\n", "\n", "#let's look at the pull distribution, plotted with normal distribution\n", "bins=np.arange(-2.0,2.5,0.5)\n", "plt.hist(y_pull,bins=bins,density=True)\n", "plt.plot(bins,stats.norm.pdf(bins))\n", "plt.xlabel('x-x$_{true}$/$\\sigma$')\n", "plt.show();"]}, {"cell_type": "markdown", "id": "6e594b79", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 4.2.4</span>\n", "\n", "Now plot your data with error bars and the true function to ensure they match. Which of the following statements is true of the data?\n", "\n", "A) None of the data points are within one uncertainty of the model\n", "\n", "B) The majority of the data points ARE NOT within one uncertainty of the model\n", "\n", "C) The majority of the data points ARE within one uncertainty of the model\n", "\n", "D) All of the data points are within one uncertainty of the model"]}, {"cell_type": "code", "execution_count": null, "id": "32ba3422", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P4.2.4\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "pass"]}, {"cell_type": "markdown", "id": "97ce6cd7", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": [">#### Follow-up 4.2.4a (ungraded)\n", ">   \n", ">Once again, try plotting this result with a legend and axis labels!\n", "\n", "<br>"]}, {"cell_type": "markdown", "id": "28d55e49", "metadata": {"tags": ["learner", "catsoop_02", "md"]}, "source": ["<h3>Making a Model with lmfit</h3>\n", "\n", "Now, we want to use the function `model_fn` as the model in a fit to the data using `lmfit`. This involves making a `Model` object, and giving it a `Parameters` object to describe the model parameters.\n", "\n", "Each parameter has `min`, `max`, and `value` arguments that specify the minimum allowable value, the maximum allowable value, and the initial value, respectively. None of these are required, but it's often a good idea to put them in if you expect your values to be within a certain range."]}, {"cell_type": "code", "execution_count": null, "id": "7b2d14f4", "metadata": {"tags": ["learner", "py", "catsoop_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.2-runcell02\n", "\n", "from lmfit import Model, Parameters\n", "\n", "model = Model(model_fn)\n", "\n", "params = Parameters()\n", "params.add('k', min=0, max=5, value=1)\n", "params.add('a', min=0, max=3, value=2)"]}, {"cell_type": "markdown", "id": "735bbf99", "metadata": {"tags": ["learner", "catsoop_02", "md"]}, "source": ["Finally, let's run the fit! We reuse the code from the previous example, but we have to pass `params` into the fit function this time.\n", "\n", "<!--\n", "#initial code\n", "result = model.fit(y, params, x=x, weights=1/yerr);\n", "\n", "print(result.fit_report())\n", "\n", "result.plot();\n", "-->"]}, {"cell_type": "code", "execution_count": null, "id": "4ba039cf", "metadata": {"tags": ["learner", "py", "catsoop_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.2-runcell03\n", "\n", "result = model.fit(y_data, params, x=x, weights=1/y_unc);\n", "\n", "print(result.fit_report())\n", "\n", "result.plot();"]}, {"cell_type": "markdown", "id": "f1217969", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 4.2.5</span>\n", "\n", "Do your fit results agree with the true values, within the uncertainty of the fit? Which was determined with higher precision: $k$ or $a$? Does this make sense given your knowledge of the model function and the uncertainties?"]}, {"cell_type": "markdown", "id": "4e743962", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_4_3'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">P4.3 A More Complicated Model</h2>    \n", "\n", "| [Top](#section_4_0) | [Previous Section](#section_4_2) | [Problems](#problems_4_3) | [Next Section](#section_4_4) |\n"]}, {"cell_type": "markdown", "id": "f6501f08", "metadata": {"tags": ["learner", "catsoop_03", "md"]}, "source": ["<h3>Overview</h3>\n", "\n", "Project 1, which you will do in Week 4, uses `lmfit` to fit data with a relatively complicated gravitational wave model. Here we'll grapple with some of the same challenges you'll encounter there.\n", "\n", "We'll look at a fit model which is similar to a black hole merger waveform."]}, {"cell_type": "markdown", "id": "0bec1871", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='problems_4_3'></a>     \n", "\n", "| [Top](#section_4_0) | [Restart Section](#section_4_3) | [Next Section](#section_4_4) |\n"]}, {"cell_type": "markdown", "id": "463ce2a5", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 4.3.1</span>\n", "\n", "The code below defines a complicated model function with six parameters. Add code to generate true data from the function and store it as `yi_true` (do not include any uncertainties). Then plot your data (this time, since we don't have any uncertainties, you should just use `plt.plot`) to see what the waveform looks like. The quantities $y_i$ represent a strain as measured in an experiment.\n", "\n", "You will see oscillations, with the data showing a number of maxima, not all of which are the same height. The term \"global\" is applied to the maximum with the largest amplitude, while the remaining maxima with smaller amplitudes are called \"local\". \n", "\n", "Within the data range $-15.5 \\leq x_i \\leq 4.5$, how many local maxima are reached before the global maximum is reached? Do not include the global maximum itself in your count.\n"]}, {"cell_type": "code", "execution_count": null, "id": "f89cfa7b", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P4.3.1\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "np.random.seed(0x98a09fe)\n", "\n", "def complicated_model_fn(x, lambda_plus, lambda_minus, max_amp, omega_0, omega_max, omega_sigma):\n", "    omega = (omega_max - omega_0) * (np.exp(-np.minimum(x, 0)**2 / omega_sigma)) + omega_0\n", "    lambdas = np.array([lambda_plus if xvalue > 0 else lambda_minus for xvalue in x])\n", "    amplitude = max_amp * np.exp(-abs(x) / lambdas)\n", "    return amplitude * np.cos(omega * x)\n", "\n", "\n", "LAMBDA_PLUS_TRUE = 1.0\n", "LAMBDA_MINUS_TRUE = 4\n", "MAX_AMP_TRUE = 1.2\n", "OMEGA_0_TRUE = 3.0\n", "OMEGA_MAX_TRUE = 6.0\n", "OMEGA_SIGMA_TRUE = 4.0\n", "\n", "\n", "#######FILL IN CODE BELOW#######\n", "xi = np.linspace(-15.5, 4.5, 200)\n", "yi_true = # your code here\n", "\n", "# add your plotting code here\n", "plt.xlabel(\"Time (s)\")\n", "plt.ylabel(\"Strain\");"]}, {"cell_type": "markdown", "id": "4dd98652", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": [">#### Follow-up 4.3.1a (ungraded)\n", ">   \n", ">Is this what you would expect a black hole merger to look like? Why or why not?\n", "\n", "<br>"]}, {"cell_type": "markdown", "id": "85aae519", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_4_4'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">P4.4 Fitting Data Containing Noise</h2>    \n", "\n", "| [Top](#section_4_0) | [Previous Section](#section_4_3) | [Problems](#problems_4_4) | [Next Section](#section_4_5) |\n"]}, {"cell_type": "markdown", "id": "ee4fac38", "metadata": {"tags": ["learner", "md", "catsoop_04"]}, "source": ["<h3>Adding Noise</h3>\n", "\n", "Let's add some noise to our data to reflect the kind of measurement we might make in the \"real world.\" Since we are trying to find an oscillating signal in the data, let's make our life particularly difficult by adding some sine functions as noise to the `yi_true` values.\n", "\n", "Specifically, we will generate some noise composed of 10 sines with frequencies randomly taken from a uniform distribution ranging from $0.5$ to $7$, phases taken randomly from a uniform distribution ranging from $0$ to $2\\pi$, and amplitudes set so that the sum of all the noise amplitudes is on average equal to the maximum amplitude of the signal."]}, {"cell_type": "code", "execution_count": null, "id": "006dc807", "metadata": {"tags": ["learner", "py", "catsoop_04", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.4-runcell01\n", "\n", "np.random.seed(0x98a09fe)\n", "\n", "MAX_AMP_TRUE = 1.2\n", "NUMBER_SINES_TO_ADD = 10\n", "\n", "noise_frequencies = 0.5 + 7 * np.random.random(NUMBER_SINES_TO_ADD)\n", "noise_phases = 2 * np.pi * np.random.random(NUMBER_SINES_TO_ADD)\n", "noise_amplitudes = 2 * MAX_AMP_TRUE / NUMBER_SINES_TO_ADD * np.random.random(NUMBER_SINES_TO_ADD)\n", "# The above line sets noise amplitudes so that the sum of all the noise amplitudes is on average\n", "# equal to the maximum amplitude of the signal.\n", "\n", "yi = yi_true.copy()# yi contains the original function\n", "\n", "for freq, phase, amplitude in zip(noise_frequencies, noise_phases, noise_amplitudes):\n", "    yi += amplitude * np.sin(phase + freq * xi)\n", "\n", "plt.plot(xi, yi, label='Data')\n", "plt.plot(xi, yi_true, label='True')\n", "plt.xlabel(\"Time (s)\")\n", "plt.ylabel(\"Strain\")\n", "plt.legend();"]}, {"cell_type": "markdown", "id": "760550e0", "metadata": {"tags": ["learner", "md", "catsoop_04"]}, "source": ["Now, we want to try to fit this noisy data using the same `complicated_model_fn` that we used to generate `yi_true`. As before, we need to create an `lmfit` model and define the parameters needed for this function."]}, {"cell_type": "code", "execution_count": null, "id": "073c58a9", "metadata": {"tags": ["learner", "py", "catsoop_04", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.4-runcell02\n", "\n", "from lmfit import Model, Parameters\n", "np.random.seed(0x98a09fe)\n", "\n", "model = Model(complicated_model_fn)\n", "\n", "params = Parameters()\n", "params.add('lambda_plus', min=0.1, max=5, value=1.1)\n", "params.add('lambda_minus', min=0.1, max=5, value=1)\n", "params.add('max_amp', min=0, max=2, value=1)\n", "params.add('omega_0', min=0, max=5, value=1)\n", "params.add('omega_max', min=0, max=10, value=1)\n", "params.add('omega_sigma', min=0, max=5, value=1)"]}, {"cell_type": "markdown", "id": "d40708eb", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='problems_4_4'></a>     \n", "\n", "| [Top](#section_4_0) | [Restart Section](#section_4_4) | [Next Section](#section_4_5) |"]}, {"cell_type": "markdown", "id": "a0a32a17", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 4.4.1</span>\n", "\n", "Use `lmfit` with the `model` and `params` variables created in the code shown above to try to fit the new signal with added noise. Do not include a weights argument, so that the uncertainties on all data points are equal. Remember to plot the result and get the fit report.\n", "\n", "What are the best fit value of `omega_0` and its uncertainty (the number after the \"$+/-$\") found by the fit? Enter your answer as a list of numbers `[omega_0, omega_unc]` with precision 1e-3."]}, {"cell_type": "code", "execution_count": null, "id": "fd74dbc1", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P4.4.1\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "np.random.seed(0x98a09fe)\n", "pass"]}, {"cell_type": "markdown", "id": "56d05e39", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": [">#### Follow-up 4.4.1a (ungraded)\n", ">   \n", ">Are you happy with the fit? Why or why not? What are the uncertainties in the other fit parameters?\n", "\n", "<br>"]}, {"cell_type": "markdown", "id": "0e527784", "metadata": {"tags": ["learner", "md", "catsoop_04"]}, "source": ["<h3>A Better Fit</h3>\n", "\n", "Let's try to get better fit results. We'll do this by running the fit many times with different initial values and taking the best fit among them."]}, {"cell_type": "markdown", "id": "2b5794b5", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": [">#### Follow-up 4.4.1b (ungraded)\n", ">   \n", ">First, write a function that generates a new `params` object with initial values chosen randomly in the ranges given in the `params_min_max` dictionary.\n", ">\n", ">Try to complete the code below, and click the button to reveal the answer.\n", "\n", "<br>"]}, {"cell_type": "code", "execution_count": null, "id": "4766cfab", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P4.4.1b (UNGRADED)\n", "# Use this cell for drafting your work\n", "\n", "from lmfit import Model, Parameters\n", "np.random.seed(0x98a09fe)\n", "\n", "params_min_max = {\n", "    'lambda_plus': (0.1, 5),\n", "    'lambda_minus': (0.1, 5),\n", "    'max_amp': (0, 2),\n", "    'omega_0': (0, 5),\n", "    'omega_max': (0, 10),\n", "    'omega_sigma': (0, 5)\n", "}\n", "params_trues = {\n", "    'lambda_plus': LAMBDA_PLUS_TRUE,\n", "    'lambda_minus': LAMBDA_MINUS_TRUE,\n", "    'max_amp': MAX_AMP_TRUE,\n", "    'omega_0': OMEGA_0_TRUE,\n", "    'omega_max': OMEGA_MAX_TRUE,\n", "    'omega_sigma': OMEGA_SIGMA_TRUE\n", "}\n", "\n", "def get_random_params():\n", "    for param_name, min_max_tuple in params_min_max.items():\n", "        # This loop gives param_name e.g. 'lambda_plus', min_max_tuple e.g. (0.1, 5)\n", "        #your code here\n", "        pass\n", "    \n", "#TESTING: should = [3.47943488]\n", "params = get_random_params()\n", "print('omega_max initial value:', params['omega_max'].value)"]}, {"cell_type": "markdown", "id": "d2693370", "metadata": {"tags": ["learner", "md", "catsoop_04"]}, "source": ["Now we write a function that fits the data using these random parameters, returning the chi squared value and the fit result. Run the solution code to the previous follow-up question to define the `get_random_params()` function."]}, {"cell_type": "code", "execution_count": null, "id": "f26b06f6", "metadata": {"tags": ["learner", "py", "catsoop_04", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.4-runcell03\n", "\n", "np.random.seed(0x98a09fe)\n", "\n", "#DEFINE get_random_params() if you have not done so already\n", "\n", "def fit(empty_arg):\n", "    model = Model(complicated_model_fn)\n", "    params = get_random_params()\n", "    result = model.fit(yi, params, x=xi)\n", "    return result.chisqr, result"]}, {"cell_type": "markdown", "id": "4c1582f8", "metadata": {"tags": ["learner", "md", "catsoop_04"]}, "source": ["Perform this fit 250 times (so that you get many different sets of random parameters) and store the results for all 250 in the results array."]}, {"cell_type": "code", "execution_count": null, "id": "06b0ffa2", "metadata": {"tags": ["learner", "py", "catsoop_04", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.4-runcell04\n", "\n", "np.random.seed(0x98a09fe)\n", "\n", "#from multiprocessing import Pool\n", "#with Pool() as pool:\n", "#    results = pool.map(fit, np.zeros(NUM_FITS))\n", "\n", "NUM_FITS = 250\n", "results=np.array([])\n", "for i in range(NUM_FITS):\n", "    tmpchi2,tmpresult = fit(empty_arg='')\n", "    results = np.append(results,[tmpchi2,tmpresult])\n", "    \n", "    #uncomment to indicate how long this takes\n", "    #print('done with: ', i +1, \" out of \", NUM_FITS)\n", "    \n", "results = results.reshape(NUM_FITS,2)"]}, {"cell_type": "code", "execution_count": null, "id": "d87a0ed9", "metadata": {"tags": ["learner", "py", "catsoop_04", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.4-runcell05\n", "\n", "#optionally run this cell to examine/understand the results\n", "\n", "print('shape of results:',results.shape)\n", "print('the first element of results is the chisq of the first fit:',results[0][0])\n", "print('the second element of results is lmfit result of the first fit:',results[0][1])\n", "#print('output of all chisq results:',results[:,0])"]}, {"cell_type": "markdown", "id": "0ce2178b", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 4.4.2</span>\n", "\n", "Sort the `results` array generated above by the chi squared value from lowest to highest (the chi squared value is the first element of every entry in `results`).\n", "\n", "The fit result with the lowest chi squared value will now be the first entry of the sorted `results`. Note that the second element of this first entry in `results` is the fit result object, so you can use `plot()` to show what it looks like.\n", "\n", "What is the value of the lowest Chi-Squared? Enter your answer as a number with precision 1e-3."]}, {"cell_type": "code", "execution_count": null, "id": "3801f318", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P4.4.2\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "np.random.seed(0x98a09fe)\n", "pass"]}, {"cell_type": "markdown", "id": "b7acbe0c", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": [">#### Follow-up 4.4.2a (ungraded)\n", ">   \n", ">First, write a function that generates a new `params` object with initial values chosen randomly in the ranges given in the `params_min_max` dictionary.\n", ">\n", ">How many iterations of the fitting process resulted in the same (approximately) `chisq` value?\n", "\n", "<br>"]}, {"cell_type": "markdown", "id": "53a56d63", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 4.4.3</span>\n", "\n", "For each parameter of the best fit result, display the true value (stored in `params_trues`), the fit value, the fit uncertainty, and the difference between the true and fit values divided by the fit uncertainty.\n", "\n", "What is the absolute value of the difference between the true and fit values divided by the fit uncertainty for the `omega_max` parameter? Enter your answer as a number with precision 1e-3."]}, {"cell_type": "code", "execution_count": null, "id": "d3603680", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P4.4.3\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "pass"]}, {"cell_type": "markdown", "id": "1d766899", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": [">#### Follow-up 4.4.3a (ungraded)\n", ">       \n", ">Are you happier with the fit now? What are some other things you might do to get an even better result?"]}, {"cell_type": "markdown", "id": "23e57cb4", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": [">#### Follow-up 4.4.3b (ungraded)\n", ">       \n", ">We have been trying to fit a known physical model to simulated data, assuming we don't know anything about the nature of the noise. Let's now consider trying to add terms to our model to ALSO model the noise (why not?). This will give you practice with changing the model.\n", ">\n", ">In the new function below, we have added sine terms to the original model, whose paramters we will now try to fit. Run the code and analyze the output.\n", ">\n", ">Note that this fit is much slower, so only 50 different attempts are made. You might want to try being patient and running more. How would you characterize the fit now?"]}, {"cell_type": "code", "execution_count": null, "id": "41a98192", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P4.4-runcell06\n", "\n", "def extra_complicated_model_fn(x, lambda_plus, lambda_minus, max_amp, omega_0, omega_max, omega_sigma,noise1,freq1,phase1,noise2,freq2,phase2,noise3,freq3,phase3,noise4,freq4,phase4):\n", "    omega = (omega_max - omega_0) * (np.exp(-np.minimum(x, 0)**2 / omega_sigma)) + omega_0\n", "    lambdas = np.array([lambda_plus if xvalue > 0 else lambda_minus for xvalue in x])\n", "    amplitude = max_amp * np.exp(-abs(x) / lambdas)\n", "    output = amplitude * np.cos(omega * x) + noise1*np.sin(freq1*x+phase1)\n", "    output = output + noise2*np.sin(freq2*x+phase2)\n", "    output = output + noise3*np.sin(freq3*x+phase3)\n", "    output = output + noise4*np.sin(freq4*x+phase4)\n", "    return output\n", "\n", "params_min_max = {\n", "    'lambda_plus': (0.1, 5),\n", "    'lambda_minus': (0.1, 5),\n", "    'max_amp': (0, 2),\n", "    'omega_0': (0, 5),\n", "    'omega_max': (0, 10),\n", "    'omega_sigma': (0, 5),\n", "    'noise1': (0, 1),\n", "    'noise2': (0, 1),\n", "    'noise3': (0, 1),\n", "    'noise4': (0, 1),\n", "    'freq1': (1, 7),\n", "    'freq2': (1, 7),\n", "    'freq3': (1, 7),\n", "    'freq4': (1, 7),\n", "    'phase1': (0, 2*np.pi),\n", "    'phase2': (0, 2*np.pi),\n", "    'phase3': (0, 2*np.pi),\n", "    'phase4': (0, 2*np.pi)\n", "}\n", "#params=model.make_params()\n", "#print(params)\n", "def get_random_params():\n", "    params = Parameters()\n", "    for p, (p_min, p_max) in params_min_max.items():\n", "        value = p_min + (p_max - p_min) * np.random.random(1)\n", "        params.add(p, min=p_min, max=p_max, value=value)\n", "    return params\n", "\n", "def fit(empty_arg):\n", "    model = Model(extra_complicated_model_fn)\n", "    params = get_random_params()\n", "    result = model.fit(yi, params, x=xi)\n", "    return result.chisqr, result\n", "\n", "NUM_FITS = 50\n", "results=np.array([])\n", "for i in range(NUM_FITS):\n", "    tmpchi2,tmpresult = fit(empty_arg='')\n", "    results = np.append(results,[tmpchi2,tmpresult])\n", "    #uncomment to indicate how long this takes\n", "    #print('done with: ', i+1 , \" out of \", NUM_FITS)\n", "    \n", "results = results.reshape(NUM_FITS,2)\n", "\n", "results = sorted(results, key=lambda x:x[0])\n", "\n", "print('omega_0:',results[0][1].params['omega_0'].value)\n", "print('chisq:',results[0][1].chisqr)\n", "print('redchisq:',results[0][1].redchi)\n", "\n", "\n", "results[0][1].plot();\n", "for param, info in results[0][1].params.items():\n", "    print(f\"{param}:\\tTrue: {float(params_trues[param])}\\tFit: {info.value} +/- {info.stderr} + \\t Sigmas: {abs(float(info.value) - float(params_trues[param])) / float(info.stderr)}\")\n"]}, {"cell_type": "markdown", "id": "7d3cf0dd", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_4_5'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">P4.5 Interpreting Common Errors with LMFIT</h2>   \n", "\n", "| [Top](#section_4_0) | [Previous Section](#section_4_4) | [Problems](#problems_4_5) |\n"]}, {"cell_type": "markdown", "id": "123a03e6", "metadata": {"tags": ["learner", "md", "catsoop_05"]}, "source": ["<h3>Overview</h3>\n", "\n", "Let's walk through a few common errors you might encounter when using `lmfit`, so you know what causes them and have options to address them."]}, {"cell_type": "markdown", "id": "22950215", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='problems_4_5'></a>   \n", "\n", "| [Top](#section_4_0) | [Restart Section](#section_4_5) |\n"]}, {"cell_type": "markdown", "id": "b2e5e622", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 4.5.1 Error 1</span>\n", "\n", "Let's run some code that we had previously, but take out the parameter limits on $k$. (Make sure you run the previous relevant blocks of code before running this one.)\n", "\n", "What caused the error? Choose the best option from the list below:\n", "\n", "- The model did not have a parameter `k` to start with\n", "- The fit did not converge because of a timeout error\n", "- The fit had a parameter go to infinity, which generated `NaN` values"]}, {"cell_type": "code", "execution_count": null, "id": "92a1c6bf", "metadata": {"tags": ["py", "learner", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P4.5.1\n", "\n", "from lmfit import Model, Parameters\n", "\n", "x = np.linspace(0.1, np.pi, 20)\n", "y = model_fn(x, TRUE_K, TRUE_A)\n", "y_unc = 0.1 + 0.4 * np.random.random(len(x))\n", "y = y + np.random.randn(len(x)) * y_unc\n", "\n", "model = Model(model_fn)\n", "params = Parameters()\n", "params.add('k')\n", "#with ranges\n", "#params.add('k', min=0, max=5, value=1)\n", "params.add('a', min=0, max=3, value=2)\n", "\n", "result = model.fit(y, params, x=x, weights=1/y_unc);\n", "\n", "print(result.fit_report())\n", "\n", "result.plot();\n", "\n", "#THROWS AN ERROR"]}, {"cell_type": "markdown", "id": "1774af5a", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 4.5.2 Error 2</span>\n", "\n", "Let's try fitting the wrong model to some data. We'll generate data according to the function $f(x)=x^2$, but fit a Gaussian model instead.\n", "\n", "Was an error thrown? Look more closely at the fit report. What is the warning line? Choose from the following:\n", "\n", "- Warning: uncertainties could not be estimated.\n", "- Warning: fit did not converge.\n", "- Warning: too many parameters.\n", "- Warning: model function mismatch.\n"]}, {"cell_type": "code", "execution_count": null, "id": "b21a4ecf", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P4.5.2\n", "\n", "import numpy as np\n", "from lmfit.models import GaussianModel\n", "\n", "np.random.seed(2)\n", "\n", "# Quadratic data\n", "xi = np.array([-2, -1, 0, 1, 2])\n", "yerr = np.array([0.3, 0.4, 0.45, 0.35, 0.6])\n", "yi = xi**2 +yerr*np.random.normal(xi.shape)\n", "\n", "# Gaussian model\n", "model = GaussianModel()\n", "\n", "results = model.fit(yi, x=xi, weights = 1/yerr);\n", "\n", "print(results.fit_report())\n", "results.plot();"]}, {"cell_type": "markdown", "id": "ac6dcdc0", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 4.5.3 Error 3</span>\n", "\n", "The fitting algorithms in `lmfit` rely on the fact that your model function needs to handle `np` arrays. What happens if yours doesn't?\n", "\n", "Run the code below. What was the error? Choose from the following:\n", "\n", "- \u2018numpy.ndarray object is not callable\u2019 \n", "- The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()\n", "- invalid syntax\n", "- setting an array element with a sequence"]}, {"cell_type": "code", "execution_count": null, "id": "4bd76c4d", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P4.5.3\n", "\n", "import numpy as np\n", "from lmfit import Model, Parameters\n", "\n", "TRUE_HEIGHT = 1.0\n", "\n", "def heaviside(x, height):\n", "    if x > 0:\n", "        return height\n", "    return 0.0\n", "\n", "xi = np.linspace(-5, 5, 10)\n", "try:\n", "    yi = heaviside(xi, TRUE_HEIGHT)\n", "except:\n", "    yi = np.array([heaviside(x, TRUE_HEIGHT) for x in xi])\n", "\n", "yerr = np.random.random(len(xi)) * 0.4 + 0.1\n", "yi += np.random.randn(len(xi)) * yerr\n", "\n", "model = Model(heaviside)\n", "params = Parameters()\n", "params.add('height', min=0.1, max=10, value=2)\n", "\n", "results = model.fit(yi, params, x=xi, weights = 1/yerr);\n", "\n", "print(results.fit_report())\n", "\n", "results.plot();\n", "\n", "#THROWS AN ERROR"]}, {"cell_type": "markdown", "id": "603cf20f", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": [">#### Follow-up 4.5.3a (ungraded)\n", ">       \n", ">Modify the model function so that it works for numpy arrays. There are several ways to do this, with varying levels of computational speed.\n", ">\n", ">Remember that, as a fallback, you can always explicitly turn a python list `list` into a numpy array by calling `np.array(list)`.\n", "\n"]}], "metadata": {"celltoolbar": "Tags", "kernelspec": {"display_name": "Python 3 (ipykernel)", "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.9.9"}}, "nbformat": 4, "nbformat_minor": 5}