8S50x_P06_runcells.ipynb

{"cells": [{"cell_type": "markdown", "id": "e9aea549", "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 6: Matched Filtering Part I - Time Domain</h1>\n"]}, {"cell_type": "markdown", "id": "4ee4bd26", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_6_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\">P6.0 Overview</h2>\n"]}, {"cell_type": "markdown", "id": "36ce3aeb", "metadata": {"tags": ["md", "learner", "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_6_1\">P6.1 What is Matched Filtering?</a>\n", "        </td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#problems_6_1\">P6.1 Problems</a>\n", "        </td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_6_2\">P6.2 Fitting in the Time Domain: Part I</a>\n", "        </td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#problems_6_2\">P6.2 Problems</a>\n", "        </td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_6_3\">P6.3 Fitting in the Time Domain: Part II</a>\n", "        </td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#problems_6_3\">P6.3 Problems</a>\n", "        </td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_6_4\">P6.4 Sweeping the Time Window</a>\n", "        </td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#problems_6_4\">P6.4 Problems</a></td>\n", "    </tr>\n", "</table>"]}, {"cell_type": "code", "execution_count": null, "id": "409350fe", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P6.0-runcell00\n", "\n", "!pip install lmfit"]}, {"cell_type": "code", "execution_count": null, "id": "c2897cb9", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P6.0-runcell01\n", "\n", "import numpy as np                 #https://numpy.org/doc/stable/\n", "import matplotlib.pyplot as plt    #https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.html\n", "from mpl_toolkits import mplot3d   #https://matplotlib.org/2.0.2/mpl_toolkits/mplot3d/tutorial.html\n", "\n", "import scipy.stats as stats        #https://docs.scipy.org/doc/scipy/reference/stats.html\n", "from scipy.stats import chisquare  #https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chisquare.html\n", "\n", "import lmfit\n", "from lmfit import Model,Parameters #https://lmfit.github.io/lmfit-py/parameters.html\n", "                                     #https://lmfit.github.io/lmfit-py/model.html#lmfit.model.Model\n", "\n", "from scipy.io.wavfile import write   #https://docs.scipy.org/doc/scipy/reference/generated/scipy.io.wavfile.write.html\n"]}, {"cell_type": "code", "execution_count": null, "id": "e9d56913", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P6.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"]}, {"cell_type": "markdown", "id": "e394c30f", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_6_1'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">P6.1 What is Matched Filtering?</h2>    \n", "\n", "| [Top](#section_6_0) | [Previous Section](#section_6_0) | [Problems](#problems_6_1) | [Next Section](#section_6_2) |\n"]}, {"cell_type": "code", "execution_count": null, "id": "03e6d4b5", "metadata": {"tags": ["learner", "py", "catsoop_01", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P6.1-runcell01\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, time, lambda_plus, lambda_minus, max_amp, omega_0, omega_max, omega_sigma):\n", "    omega = (omega_max - omega_0) * (np.exp(-np.minimum(x - time, 0)**2 / omega_sigma)) + omega_0\n", "    lambdas = np.array([lambda_plus if xvalue > time else lambda_minus for xvalue in x])\n", "    amplitude = max_amp * np.exp(-abs(x - time) / lambdas)\n", "    return amplitude * np.cos(omega * (x-time))\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", "TIME_TRUE = 50.0\n", "\n", "xi = np.linspace(TIME_TRUE-15, TIME_TRUE+5, 200)\n", "yi_true = complicated_model_fn(xi, TIME_TRUE, LAMBDA_PLUS_TRUE, LAMBDA_MINUS_TRUE, MAX_AMP_TRUE,\n", "                               OMEGA_0_TRUE, OMEGA_MAX_TRUE, OMEGA_SIGMA_TRUE)\n", "\n", "plt.plot(xi, yi_true)\n", "plt.xlabel(\"Time (s)\")\n", "plt.ylabel(\"Strain\");"]}, {"cell_type": "code", "execution_count": null, "id": "786b3a30", "metadata": {"tags": ["learner", "py", "catsoop_01", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P6.1-runcell02\n", "\n", "np.random.seed(0x98a09fe)\n", "#np.random.seed(908)\n", "\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", "sample_spacing = 0.1\n", "xi = np.arange(-128, 128, sample_spacing)#times\n", "yi = np.zeros_like(xi)#data\n", "\n", "#Adding Noise\n", "for freq, phase, amplitude in zip(noise_frequencies, noise_phases, noise_amplitudes):\n", "    yi += amplitude * np.sin(phase + freq * xi)\n", "   \n", "#Adding Data\n", "signal = complicated_model_fn(xi, TIME_TRUE, LAMBDA_PLUS_TRUE, LAMBDA_MINUS_TRUE, MAX_AMP_TRUE,\n", "                               OMEGA_0_TRUE, OMEGA_MAX_TRUE, OMEGA_SIGMA_TRUE)\n", "\n", "yi+=signal\n", "\n", "plt.figure(figsize=(16, 5))\n", "plt.plot(xi, yi)\n", "plt.title(\"Signal plus noise\")\n", "plt.xlabel(\"Time (s)\")\n", "plt.ylabel(\"Strain\")\n", "plt.show()\n", "\n", "plt.figure(figsize=(16, 5))\n", "plt.plot(xi, yi)\n", "plt.plot(xi, signal)\n", "plt.title(\"Signal plus noise\")\n", "plt.xlabel(\"Time (s)\")\n", "plt.ylabel(\"Strain\")\n", "plt.show()\n", "\n", "plt.plot(xi, yi)\n", "plt.plot(xi, signal)\n", "plt.title(\"Signal plus noise\")\n", "plt.xlim(35,55)\n", "plt.xlabel(\"Time (s)\")\n", "plt.ylabel(\"Strain\")\n", "plt.show()"]}, {"cell_type": "markdown", "id": "54bcd7cd", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='problems_6_1'></a>     \n", "\n", "| [Top](#section_6_0) | [Restart Section](#section_6_1) | [Next Section](#section_6_2) |\n"]}, {"cell_type": "markdown", "id": "925de8e1", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 6.1.1</span>\n", "\n", "In the next four problems, we will generate some noise and create a SNR (signal to noise ratio) plot in order to identify the time at which the signal exists. Since we already know the signal and the noise separately, we can implement a naive approach to finding the signal time. We will calculate the SNR and then simply take the time at which its maximum occurs as the time for the signal event. Your goal is to explore how well this crude method works.\n", "\n", "First, generate some noise composed of 1,000 sines (100 times as many as above!) with frequencies randomly taken from a normal distribution with mean at $\\mu=0.8$ and standard deviation of $\\sigma =5$, phases taken from a random 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.\n", "\n", "In a previous assignment, we did something similar with only 10 sines, but we sampled from slightly different distributions. You can refer back to that previous example, if necessary."]}, {"cell_type": "code", "execution_count": null, "id": "65eba598", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P6.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", "MAX_AMP_TRUE = 1.2\n", "SAMPLE_SPACING = 0.1\n", "NUMBER_SINES_TO_ADD = 1000\n", "\n", "xi = np.arange(0, 128, SAMPLE_SPACING)#times\n", "\n", "def generate_noise(xi):\n", "  np.random.seed(908)\n", "  yi_noise = np.zeros_like(xi)\n", "\n", "  noise_frequencies = 0 #YOUR CODE HERE\n", "  noise_phases = 0 #YOUR CODE HERE\n", "  noise_amplitudes = 0 #YOUR CODE HERE\n", "\n", "  #Adding Noise\n", "  for freq, phase, amplitude in zip(noise_frequencies, noise_phases, noise_amplitudes):\n", "      yi_noise += amplitude * np.sin( phase + freq * xi)\n", "\n", "  return yi_noise\n", "\n", "plt.figure(figsize=(16, 5))\n", "plt.xlabel(\"Time (s)\")\n", "plt.ylabel(\"Strain\")\n", "plt.plot(xi, generate_noise(xi))\n", "plt.show()\n"]}, {"cell_type": "markdown", "id": "12dc2f2e", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 6.1.2</span>\n", "\n", "Now, we want to inject a signal waveform on top of the noise signal we just generated. Then, we will use a very simple approach to try to \"find\" that signal, namely looking for the location in the total amplitude which has the maximum value. Recall that the maximum amplitude of the merger signal is the parameter `TRUE_TIME`. So,  we can determine how well we \"found\" the injected signal by checking if the time we found in the total data is consistent with the `TRUE_TIME` of the injected waveform. \n", "\n", "To begin, create a set of 50 signals of the form shown earlier. Except for the `TRUE_TIME`, each waveform will have identical parameters as shown in the code below. The `TRUE_TIME` parameters of the injected signals should range over every whole second in the range [50, 100) (i.e., you should generate signals with `TRUE_TIME`=50, `TRUE_TIME`=51, ... `TRUE_TIME`=99).\n", "\n", "For each merger signal, inject it on top of the noise generated earlier and try to find the time at which the injection happens by looking for the maximum amplitude in the combined data. Lastly, make a plot of the true (injected) time, vs the time at which you found the maximum. Notice the `scale` parameter in the function `get_max_times`. This option can be used to rescale the size of the merger signal before it gets injected (in this case, `scale=1.0`, so we used the unmodified merger waveform). \n", "\n", "HINT: Take a look at the `np.argmax()` function.\n"]}, {"cell_type": "code", "execution_count": null, "id": "71058141", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P6.1.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", "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", "def complicated_model_fn(x, time, lambda_plus, lambda_minus, max_amp, omega_0, omega_max, omega_sigma):\n", "    omega = (omega_max - omega_0) * (np.exp(-np.minimum(x - time, 0)**2 / omega_sigma)) + omega_0\n", "    lambdas = np.array([lambda_plus if xvalue > time else lambda_minus for xvalue in x])\n", "    amplitude = max_amp * np.exp(-abs(x - time) / lambdas)\n", "    return amplitude * np.cos(omega * (x-time))\n", "\n", "def get_max_times(xi, yi_noise, true_times,scale=1.0,iCheck=False):\n", "    time_of_maximums = []\n", "\n", "    for t in true_times:\n", "\n", "        yi_signal = #YOUR CODE HERE\n", "        yi_test_noise = #YOUR CODE HERE\n", "        SNR = #YOUR CODE HERE\n", "        \n", "        time_of_maximums.append(xi[np.argmax(SNR)])\n", "        \n", "        if int(t) == 75 and iCheck:\n", "            plt.xlabel(\"Time (s)\")\n", "            plt.ylabel(\"Strain\")\n", "            plt.plot(xi,yi_test_noise)\n", "            plt.plot(xi,yi_signal)\n", "            plt.show()\n", "        \n", "    return time_of_maximums\n", "\n", "true_times = np.linspace(50, 100, 50)\n", "xi = np.arange(0, 128, SAMPLE_SPACING)\n", "yi_noise = generate_noise(xi)\n", "\n", "plt.plot(true_times, get_max_times(xi, yi_noise, true_times, iCheck=True), label = 'naive model')\n", "plt.xlabel('True Times (s)')\n", "plt.ylabel('Found Times (s)')\n", "plt.legend()\n", "plt.show()"]}, {"cell_type": "markdown", "id": "c5d5c590", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 6.1.3</span>\n", "\n", "The code for the previous problem plotted the found time versus the true time (true time on the x-axis, estimated time for y) and that result looked quite good. Let's check this comparison more carefully. \n", "\n", "First, modify the code below to add to the plot what a perfect algorithm would look like, namely the line y=x. You might want to give this line a different color so you can distinguish it from the data result, and plot it *after* the other line to put it on top. \n", "\n", "Then, make a second plot showing the fractional difference of the two times (found-true/true) versus the true time.\n", "\n", "How well does the results of this crude method compare to a perfect algorithm? Select the best answer below:\n", "\n", "- The crude method is an almost exact match to the ideal algorithm.\n", "\n", "- The crude method mostly gets the time at which the signal event occurs and occasionally overestimates/underestimates.\n", "\n", "- The crude method largely underestimates the time at which the signal event occurs but occasionally gets close or overestimates.\n", "\n", "- The crude method largely overestimates the time at which the signal event occurs but occasionally gets close or underestimates.\n", "\n", "- The crude method always underestimates the time at which the signal event occurs.\n", "\n", "- The crude method always overestimates the time at which the signal event occurs.\n", "\n", "- The crude method almost never gets the right answer. \n", "\n", "- The crude method doesn't work at all.\n"]}, {"cell_type": "code", "execution_count": null, "id": "db1e08fa", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P6.1.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", "true_times = np.linspace(50, 100, 50)\n", "xi = np.arange(0, 128, SAMPLE_SPACING)\n", "yi_noise = generate_noise(xi)\n", "\n", "plt.plot(true_times, get_max_times(xi, yi_noise, true_times,iCheck=True), label = 'naive model')\n", "\n", "#YOUR CODE HERE\n", "\n", "plt.xlabel('True Times (s)')\n", "plt.ylabel('Found Times (s)')\n", "plt.legend()\n", "plt.show()\n", "\n", "#YOUR CODE HERE"]}, {"cell_type": "markdown", "id": "cb0e8f8b", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 6.1.4</span>\n", "\n", "Finally, let's make the problem more challenging by scaling the injected merger signal down by a factor of 20. What happens to the crude method in this case? Select the best answer below:\n", "\n", "- The crude method is an almost exact match to the ideal algorithm.\n", "\n", "- The crude method largely underestimates the time at which the signal event occurs but occasionally gets close or overestimates.\n", "\n", "- The crude method largely overestimates the time at which the signal event occurs but occasionally gets close or underestimates.\n", "\n", "- The crude method always underestimates the time at which the signal event occurs.\n", "\n", "- The crude method always overestimates the time at which the signal event occurs.\n", "\n", "- The crude method almost never gets the right answer. \n", "\n", "- The crude method doesn't work at all."]}, {"cell_type": "code", "execution_count": null, "id": "eb318d7a", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P6.1.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", "true_times = np.linspace(50, 100, 50)\n", "xi = np.arange(0, 128, SAMPLE_SPACING)\n", "yi_noise = generate_noise(xi)\n", "\n", "# Modify the following to scale the signal down by a factor of 20, compare the the ideal model\n", "#plt.plot(true_times, get_max_times(xi, yi_noise, true_times), label = 'naive model')\n", "plt.plot(true_times, true_times, color='red', label='ideal model')\n", "plt.xlabel('True Times (s)')\n", "plt.ylabel('Found Times (s)')\n", "plt.legend()\n", "plt.show()"]}, {"cell_type": "markdown", "id": "886eded2", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_6_2'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">P6.2 Fitting in the Time Domain: Part I</h2>    \n", "\n", "| [Top](#section_6_0) | [Previous Section](#section_6_1) | [Problems](#problems_6_2) | [Next Section](#section_6_3) |\n"]}, {"cell_type": "code", "execution_count": null, "id": "daaa99f9", "metadata": {"tags": ["learner", "py", "catsoop_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P6.2-runcell01\n", "\n", "np.random.seed(0x98a09fe)\n", "\n", "def complicated_model_fn(x, time, lambda_plus, lambda_minus, max_amp, omega_0, omega_max, omega_sigma):\n", "    omega = (omega_max - omega_0) * (np.exp(-np.minimum(x - time, 0)**2 / omega_sigma)) + omega_0\n", "    lambdas = np.array([lambda_plus if xvalue > time else lambda_minus for xvalue in x])\n", "    amplitude = max_amp * np.exp(-abs(x - time) / lambdas)\n", "    return amplitude * np.cos(omega * (x-time))\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", "TIME_TRUE = 50.0\n", "\n", "xi = np.linspace(TIME_TRUE-15, TIME_TRUE+5, 200)\n", "yi_true = complicated_model_fn(xi, TIME_TRUE, LAMBDA_PLUS_TRUE, LAMBDA_MINUS_TRUE, MAX_AMP_TRUE,\n", "                               OMEGA_0_TRUE, OMEGA_MAX_TRUE, OMEGA_SIGMA_TRUE)\n", "\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", "plt.plot(xi, yi_true)\n", "plt.title(\"True Signal\")\n", "plt.xlabel(\"Time (s)\")\n", "plt.ylabel(\"Strain\")\n", "plt.show()\n", "\n", "sample_spacing = 0.1\n", "xi = np.arange(-128, 128, sample_spacing)#times\n", "yi = np.zeros_like(xi)#data\n", "\n", "#Adding Noise\n", "for freq, phase, amplitude in zip(noise_frequencies, noise_phases, noise_amplitudes):\n", "    yi += amplitude * np.sin(phase + freq * xi)\n", "\n", "#Adding Data\n", "signal= complicated_model_fn(xi, TIME_TRUE, LAMBDA_PLUS_TRUE, LAMBDA_MINUS_TRUE, MAX_AMP_TRUE,\n", "                               OMEGA_0_TRUE, OMEGA_MAX_TRUE, OMEGA_SIGMA_TRUE)\n", "yi+=signal\n", "\n", "\n", "plt.figure(figsize=(16, 5))\n", "plt.plot(xi, yi)\n", "plt.plot(xi, signal)\n", "plt.title(\"Signal plus noise\")\n", "plt.xlabel(\"Time (s)\")\n", "plt.ylabel(\"Strain\")\n", "plt.show()\n", "\n", "plt.plot(xi, yi)\n", "plt.plot(xi, signal)\n", "plt.title(\"Signal plus noise\")\n", "plt.xlim(35,55)\n", "plt.xlabel(\"Time (s)\")\n", "plt.ylabel(\"Strain\")\n", "plt.show()"]}, {"cell_type": "markdown", "id": "2b988556", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='problems_6_2'></a>     \n", "\n", "| [Top](#section_6_0) | [Restart Section](#section_6_2) | [Next Section](#section_6_3) |\n"]}, {"cell_type": "markdown", "id": "0253de1b", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 6.2.1</span>\n", "\n", "Since the signal waveform has a limited extent before and after the time of maximum amplitude, it's clear that fitting all of the data with just the merger signal function will fail miserably. Rather, we'll need to select a subset of the data to perform the fit. How much time before and after the time of maximum signal $t$ do you think the fit should include? We really only need to consider the region where the signal is rather large and therefore will be more easily separated from the noise. In practice this is something we could find systematically, for example by analyzing the data close to and far from the maximum signal point. However, for now, read the guidance below to choose an appropriate window.\n", "\n", "In what follows, only consider a 7-9 second long window, as this will include enough data to make our fits converge, but will still give few enough data points that the fits converge relatively fast. Furthermore, this window need not be symmetric around $t$, as much of the signal lies before the true time with only a short period of signal after $t$. Therefore, we want the number of earlier seconds to include (`t_before`) to be larger than the following time span (`t_after`).\n", "\n", "With these conditions, one possible choice for `[t_before, t_after]` is `[5,2]`. What is another acceptable answer, given the constraints that we outlined?\n", "\n", "Enter your answer as a list, formatted as `[t_before, t_after]`."]}, {"cell_type": "markdown", "id": "4b8a6c64", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 6.2.2</span>\n", "\n", "In order to fit the signal waveform to a set of data, we need to write a function `model_and_random_parameters(t)` that creates an `lmfit` `Model` and associated `Parameters` for  `complicated_model_fn`. As discussed above, we want to force the time of the signal to appear at an input time `t`.  To limit the range that the fit considers, we also want to constrain the parameters, with limits in the dictionary `params_min_max`. the initial values of the parameters should be chosen randomly within those given ranges.\n", "\n", "Complete the function `get_param_random_value` in the code below so that it choses random starting points for the parameters which are uniformly distributed between `p_min` and `p_max`.\n", "\n", "To check if you've done this correctly, the lines at the end of the code print out the value of one of the parameters. With the given starting seed, the answer should be about 1.51. Test this in your notebook."]}, {"cell_type": "code", "execution_count": null, "id": "4d2bbd85", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P6.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", "from lmfit import Model, Parameters\n", "    \n", "    \n", "def get_param_random_value(p_min,p_max):\n", "    #get a uniformly distributed random value between p_min and p_max\n", "    #return a float\n", "    return #YOUR CODE HERE\n", "\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", "\n", "def model_and_random_parameters(t):\n", "    model = Model(complicated_model_fn)\n", "    params = Parameters()\n", "    params.add('time', value=t, vary=False)\n", "    for p, (p_min, p_max) in params_min_max.items():\n", "        value = get_param_random_value(p_min,p_max)\n", "        params.add(p, min=p_min, max=p_max, value=value)\n", "    return model, params\n", "\n", "\n", "#TEST EXAMPLE: SHOULD = 1.51166\n", "t=0.1\n", "np.random.seed(1)\n", "print(model_and_random_parameters(t)[1].get('omega_0').value)"]}, {"cell_type": "markdown", "id": "fa313837", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 6.2.3</span>\n", "\n", "The code below includes a function `fit_once` that fits the model created in the previous problem and outputs the fit result. Remember that we only want to look at a limited subset of the data when we fit, namely the range `(t-t_before, t+t_after)`, where `t` is the specific time at which we want to look for the signal. This version uses the values `t_before = 5` and `t_after = 2`. \n", "\n", "You need to write a function `get_signal_indices` which returns the indices in the data arrays which correspond to this range of times. HINT: You may find the function `np.where()` useful.\n", "\n", "Optionally print the $\\chi^2$ result or the fit report. Note, that since a random seed is not defined, you may get different results each time you run. This is the point of choosing random initial parameters."]}, {"cell_type": "code", "execution_count": null, "id": "5a98bfdc", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P6.2.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", "import lmfit\n", "\n", "#THE WINDOW MUST BE [5,2] FOR YOUR ANSWER TO MATCH EXPECTED VALUES\n", "t_before = 5\n", "t_after = 2\n", "\n", "\n", "def get_signal_indices(xi, t, t_before, t_after):\n", "    #use np.where() to return a 1D the relevant indices\n", "    #note, the result of np.where() will be a tuple\n", "    return #YOUR CODE HERE\n", "\n", "def fit_once(xi, yi, t, t_before, t_after):\n", "    data_indices = get_signal_indices(xi, t, t_before, t_after)\n", "    data_x = xi[data_indices]\n", "    data_y = yi[data_indices]\n", "    model, params = model_and_random_parameters(t)    \n", "    result = model.fit(data_y, params, x=data_x)\n", "    return result\n", "\n", "result = fit_once(xi, yi, TIME_TRUE, t_before, t_after)\n", "result.plot();\n", "\n", "#print(\"Fit chi2 value: \", result.chisqr)\n", "#print(\"Fit chi2 probability: \",1-stats.chi2.cdf(result.chisqr,result.nfree))\n", "#print(result.fit_report())"]}, {"cell_type": "markdown", "id": "812323f3", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 6.2.4</span>\n", "\n", "Run the fit multiple times and print the $\\chi^{2}$ value and $\\chi^{2}$ probability using the following lines of code.\n", "\n", "<pre>\n", "print(\"Fit chi2 value: \", result.chisqr)\n", "print(\"Fit chi2 probability: \",1-stats.chi2.cdf(result.chisqr,result.nfree))\n", "</pre>\n", "\n", "What is the lowest $\\chi^{2}$ value that you obtain, and its corresponding $\\chi^{2}$ probability?\n", "\n", "Enter your answer as a list of numbers `[chi2, chi2_prob]` with precision 1e-3."]}, {"cell_type": "markdown", "id": "9378f2e9", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 6.2.5</span>\n", "\n", "\n", "Let's consider whether the best $\\chi^{2}$ value and its probability are reasonable numbers. In particular, what does the $\\chi^{2}$ probability say about the fit? Choose the best answer from the following options:\n", "\n", "- The fit is perfect! This is because our model is perfect. Our job is done.\n", "- The fit is too perfect, which means we should carefully consider the assumptions we have made.\n", "- The fit is okay, and we can do no better.\n", "- The fit is terrible, so we should adjust our model or the range of data that we are fitting.\n"]}, {"cell_type": "code", "execution_count": null, "id": "8a178b2b", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P6.2.5\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "import scipy.stats as stats\n", "\n", "print(\"Fit chi2 probability: \",1-stats.chi2.cdf(result.chisqr,result.nfree))"]}, {"cell_type": "markdown", "id": "98b43d78", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_6_3'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">P6.3 Fitting in the Time Domain: Part II</h2>    \n", "\n", "| [Top](#section_6_0) | [Previous Section](#section_6_2) | [Problems](#problems_6_3) | [Next Section](#section_6_4) |\n"]}, {"cell_type": "markdown", "id": "2249bf51", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='problems_6_3'></a>     \n", "\n", "| [Top](#section_6_0) | [Restart Section](#section_6_3) | [Next Section](#section_6_4) |\n"]}, {"cell_type": "markdown", "id": "5571fc1a", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 6.3.1</span>\n", "\n", "We want to modify the previous fit to use uncertainties with a value of $\\sigma=0.2$. To do this, you will run a \"weighted\" fit with `lmfit` by setting an array of weights. The code below is similar to what you saw above, but now with a function `fit_once_weighted`. You need to complete that function to use $\\sigma=0.2$.\n", "\n", "Note, the weights in `lmfit` are designed so that $w=1/\\sigma$, leading to the following:\n", "\n", "$$\\chi^{2} = \\sum_{i}\\frac{(f(x_{i})-f(x))^{2}}{\\sigma_{i}^{2}}$$\n", "\n", "With this new version of the fit, what is the $\\chi^{2}$ probability corresponding to the lowest $\\chi^{2}$ value ? Enter your answer as a number with precision 1e-3.\n", "\n", "Hint: As before, you can run the code multiple times or use a `for` loop.\n"]}, {"cell_type": "code", "execution_count": null, "id": "c2dbf2c5", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P6.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", "import lmfit\n", "\n", "def fit_once_weighted(xi, yi, t, t_before, t_after, weight=1.0):\n", "    data_indices = get_signal_indices(xi, t, t_before, t_after)\n", "    data_x = xi[data_indices]\n", "    data_y = yi[data_indices]\n", "    \n", "    weights = #YOUR CODE HERE\n", "    \n", "    model, params = model_and_random_parameters(t)\n", "    result = model.fit(data_y, params, x=data_x,weights=weights)\n", "    return result\n", "\n", "#The window must be [5,2] for your answer to match expected values\n", "t_before = 5\n", "t_after = 2\n", "\n", "unc=0.2\n", "result = fit_once_weighted(xi, yi, TIME_TRUE, t_before, t_after,1./unc)\n", "\n", "result.plot();\n", "print(\"Fit chi2 value: \", result.chisqr)\n", "print(\"Fit chi2 probability: \",1-stats.chi2.cdf(result.chisqr,result.nfree))"]}, {"cell_type": "markdown", "id": "c5879b89", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 6.3.2</span>\n", "\n", "The choice of $\\sigma=0.2$ in the previous problem was totally arbitrary. Can we come up with a strategy to compute a more realistic uncertainty for our data?  With LIGO data, this is a difficult question, since much of the wiggles from the \"Noise\" are actually understood as oscillations at certain frequencies. \n", "\n", "For now, we won't model the noise, but instead, we'll calculate uncertainties that reflect the average RMS of our noisy data by using a  signal-free region of time. \n", "\n", "The code below computes the standard deviation of the signal-free noise, but it needs to know what time range to use for that calculation. You need to complete the function `get_noise_indices` and then repeat the fit. In your noise calculation, only exclude the region that is considered in the fit.\n", "\n", "What $\\chi^{2}$ probability do you get, corresponding to the lowest $\\chi^{2}$ value? Again, run your code multiple times to obtain the min value. Enter your answer **for the $\\chi^{2}$ probability** as a number with precision 1e-3."]}, {"cell_type": "code", "execution_count": null, "id": "2c951269", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P6.3.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", "import lmfit\n", "\n", "def get_noise_indices(xi, t, t_before, t_after):\n", "    return #YOUR CODE HERE\n", "    \n", "\n", "def noise(xi, yi, t, t_before, t_after):\n", "    data_indices = get_noise_indices(xi, t, t_before, t_after)\n", "    data_y = yi[data_indices]\n", "    return np.std(data_y)\n", "\n", "t_before = 5\n", "t_after = 2\n", "\n", "unc=noise(xi, yi, TIME_TRUE, t_before, t_after)\n", "result = fit_once_weighted(xi, yi, TIME_TRUE, t_before, t_after,1./unc)\n", "result.plot();\n", "print(\"unc value: \", unc)\n", "print(\"Fit chi2 value: \", result.chisqr)\n", "print(\"Fit chi2 probability: \",1-stats.chi2.cdf(result.chisqr,result.nfree))"]}, {"cell_type": "code", "execution_count": null, "id": "27b55c56", "metadata": {"tags": ["learner", "py", "learner_chopped", "catsoop_03"]}, "outputs": [], "source": ["#>>>RUN: P6.3-runcell01\n", "\n", "#Computing uncertainty: merging bins\n", "\n", "xi_old = xi.copy()\n", "yi_old = yi.copy()\n", "xi_new = np.array([ 0.5*(xi_old[2*i]+xi_old[2*i+1]) for i in range(len(xi_old)//2) ])\n", "yi_new = np.array([ 0.5*(yi_old[2*i]+yi_old[2*i+1]) for i in range(len(yi_old)//2) ])\n", "\n", "t_before = 5\n", "t_after = 2\n", "\n", "uncout=noise(xi_new, yi_new, TIME_TRUE, t_before, t_after)\n", "result = fit_once_weighted(xi_new, yi_new, TIME_TRUE, t_before, t_after,1./uncout)\n", "result.plot();\n", "print(\"unc value: \", uncout)\n", "print(\"Fit chi2 value: \", result.chisqr)\n", "print(\"Fit chi2 probability: \",1-stats.chi2.cdf(result.chisqr,result.nfree))"]}, {"cell_type": "code", "execution_count": null, "id": "b1dde408", "metadata": {"tags": ["learner", "py", "learner_chopped", "catsoop_03"]}, "outputs": [], "source": ["#>>>RUN: P6.3-runcell02\n", "\n", "#Computing uncertainty: points 2 samples away\n", "\n", "def noise_deltat(xi, yi, t, t_before, t_after, dt=2):#dt is the size distance of the samples\n", "    data_indices = get_noise_indices(xi, t, t_before, t_after)\n", "    #print(data_indices[0][:-dt],data_indices[0][dt:])\n", "    data_y = yi[data_indices[0][:-dt]]-yi[data_indices[0][dt:]]\n", "    return np.std(data_y)\n", "\n", "t_before = 5\n", "t_after = 2\n", "\n", "uncout=noise_deltat(xi_old, yi_old, TIME_TRUE, t_before, t_after)\n", "result = fit_once_weighted(xi_old, yi_old, TIME_TRUE, t_before, t_after,1./uncout)\n", "result.plot();\n", "print(\"unc value: \", uncout)\n", "print(\"Fit chi2 value: \", result.chisqr)\n", "print(\"Fit chi2 probability: \",1-stats.chi2.cdf(result.chisqr,result.nfree))"]}, {"cell_type": "markdown", "id": "49d18a81", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_6_4'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">P6.4 Sweeping the Time Window</h2>   \n", "\n", "| [Top](#section_6_0) | [Previous Section](#section_6_3) | [Problems](#problems_6_4) | [Next Section](#section_6_5) |\n"]}, {"cell_type": "code", "execution_count": null, "id": "7e740a89", "metadata": {"scrolled": false, "tags": ["py", "learner", "catsoop_04", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P6.4-runcell01\n", "\n", "#From the above analysis, we see that an uncertainty of 0.18 is more reasonable.\n", "#This was found for deltat = 2, which should limit that assumptions we are making\n", "#about the nature of the background noise\n", "\n", "#Let's try using this uncertainty\n", "\n", "\n", "def fit_once_new(t, t_before, t_after, weight=1.0):\n", "    data_indices = get_signal_indices(xi, t, t_before, t_after)\n", "    data_x = xi[data_indices]\n", "    data_y = yi[data_indices]\n", "    weights = np.ones(len(data_x))*weight\n", "    model, params = model_and_random_parameters(t)\n", "    result = model.fit(data_y, params, x=data_x,weights=weights)\n", "    return result\n", "\n", "t_before = 5\n", "t_after = 2\n", "\n", "unc=0.18\n", "result = fit_once_new(TIME_TRUE, t_before, t_after, 1.0/unc)\n", "result.plot();\n", "print(\"Fit chi2 value: \", result.chisqr)\n", "print(\"Fit chi2 probability: \",1-stats.chi2.cdf(result.chisqr,result.nfree))"]}, {"cell_type": "markdown", "id": "e41f50bb", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='problems_6_4'></a>     \n", "\n", "| [Top](#section_6_0) | [Restart Section](#section_6_4) |\n"]}, {"cell_type": "markdown", "id": "80e66bb1", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 6.4.1</span>\n", "\n", "In the code below, the function `fit_once_new` from above is run multiple times. Complete the code to find the best result in the list and return it. Specifically, modify the `get_min_result` function. \n", "\n", "What is the lowest $\\chi^2$ value and its corresponding probability? Enter your answer as a list of numbers `[chi2, chi2_prob]` with precision 1e-3."]}, {"cell_type": "code", "execution_count": null, "id": "22a7ac66", "metadata": {"scrolled": false, "tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P6.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", "def get_min_result(results):\n", "    min_result = None\n", "    min_chisq = None\n", "    #for each result in results, set a new min_result and min_chisq\n", "    #if result.chisqr is less than the currently stored value\n", "    \n", "    #YOUR CODE HERE\n", "    \n", "    return min_result\n", "\n", "NUM_FITS = 10\n", "\n", "def fit_many(t, t_before, t_after, weight, num):\n", "  results=[]\n", "  for i in range(num):\n", "    results.append(fit_once_new(t, t_before, t_after, weight))\n", "\n", "  min_result = get_min_result(results)\n", "  return min_result\n", "\n", "t_before = 5\n", "t_after = 2\n", "unc=0.18\n", "\n", "min_result = fit_many(TIME_TRUE, t_before, t_after, 1.0/unc, NUM_FITS)\n", "\n", "min_result.plot();\n", "print(\"Fit chi2 value: \", min_result.chisqr)\n", "print(\"Fit chi2 probability: \",1-stats.chi2.cdf(min_result.chisqr,min_result.nfree))"]}, {"cell_type": "markdown", "id": "aec4f0ce", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 6.4.2</span>\n", "\n", "Next, we want to see what the fits look like for different `t` values. The code cell below calls `fit_many` for values of $t \\in [-100, 100]$, where the $t$ values are separated by $\\Delta t = 1 \\mathrm{s}$ and stores all of the fit outputs in an array named `results`. The code is complete, you just need to run it.\n", "\n", "How long does it take to do this? (pick the closest answer)\n", "\n", "A. .01 seconds\n", "\n", "B. 1 second\n", "\n", "C. 5 minutes\n", "\n", "D. 5 hours (if this is the answer you pick **something is wrong**)\n", "\n", "E. 10 days (if this is the answer you pick **something is wrong**)"]}, {"cell_type": "code", "execution_count": null, "id": "f25836f4", "metadata": {"tags": ["py", "learner", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: P6.4.2\n", "%%time\n", "\n", "results = []\n", "delta_t = 1\n", "ts = np.arange(-100, 100, delta_t)\n", "\n", "t_before = 5\n", "t_after = 2\n", "unc=0.18\n", "\n", "NUM_FITS = 6\n", "\n", "for t in ts:\n", "  if t % 20 == 0: print(t)\n", "  results.append(fit_many(t, t_before, t_after, 1.0/unc, NUM_FITS))"]}, {"cell_type": "markdown", "id": "a32b8600", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Problem 6.4.3</span>\n", "\n", "Now we need to find out for which $t$ value we get a fit that is most likely to be our signal. One way of figuring this out is by looking for which fit has the largest `max_amp` parameter, as the signal will have a higher max amplitude than the surrounding noise.\n", "\n", "Plot `max_amp` as a function of $t$ given the `results` you just calculated. Find the value of $t$ which has the largest `max_amp` and plot the corresponding fit result.\n", "\n", "Does the fit look like it could be the signal we're looking for? If yes, then enter below at what value of $t$ this was. If not, keep searching through the next highest `max_amp` values till you get something that may be signal and answer that $t$ value below."]}, {"cell_type": "code", "execution_count": null, "id": "e33c5e8b", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>PROBLEM: P6.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", "amps = #LIST OF MAXIMUM AMPLITUDES\n", "result_max_amp = #RESULT CORRESPONDING TO MAX AMP\n", "\n", "result_max_amp.plot()\n", "print(\"Time of best fit result: \", ts[np.argmax(amps)])\n", "plt.show()\n", "\n", "plt.plot(ts, amps)\n", "plt.xlabel(\"Time(s)\")\n", "plt.ylabel(\"Wave amplitudes\")\n", "plt.show()"]}, {"cell_type": "code", "execution_count": null, "id": "bb5bcbe8", "metadata": {"scrolled": false, "tags": ["py", "learner", "learner_chopped", "catsoop_04"]}, "outputs": [], "source": ["#>>>RUN: P6.4-runcell02\n", "\n", "chi2 = [r.chisqr for r in results]\n", "min_result=results[np.argmin(chi2)]\n", "min_result.plot()\n", "print(\"Time of best fit result: \", ts[np.argmin(chi2)],\" with lowest chi^2: \",min_result.chisqr)\n", "plt.show()\n", "\n", "plt.plot(ts, chi2)\n", "plt.xlabel(\"Time(s)\")\n", "plt.ylabel(\"Wave $\\chi^{2}$\")\n", "plt.show()"]}], "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}