<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Review from the previous lecture
In yesterday's Lecture 2, you learned how to use the numpy module, how to make your own functions, and how to import and export data. Below is a quick review before we move on to Lecture 3.
Remember, to use the numpy module, first it must be imported
Step1:  You can do a lot with the numpy module. Below is an example to jog your memory
Step2:  Do you remember loops? Let's use a while loop to make an array of 10 numbers. Let's have each element be increased by 2 compared with the previous element. Let's also have the first element of the array be 1.
Step3:  There's your quick review of numpy and a while loop. Now we can move on to the content of Lecture 3.
A. Data Introspection
One of the most common tasks in experimental physics is trying to model experimental data with a function. This lecture will walk through how to accomplish this with python. First, let's import some data.
Step4:  Let's plot it to see what it looks like!
In python there are multiple plotting libraries, but the moset common one is matplotlib, and that is the one we will be using today.
Step6:  Usually, you will want to have some theory-based motivation for the function you choose to model some set of data, but for this example, we don't know anything about the data other than the points themselves. In this type of situation, trying to fit a simple function to the data is not a bad first step in trying to understand it. What function do you think might fit this data based on how it looks in the plot?
<details>
<summary>Answer</summary>
It looks like the data is shaped like a normal (gaussian) distribution, so let's try to fit it to that! First, let's define a gaussian function for fitting.
The equation for a Gaussian curve is the following
Step7:  What values of these parameters do you think will match the data above?
Step8:  Let's try plotting the function the function with these parameters next to the data! For this, we should define some evenly-spaced x-values to calculate the function at using np.linspace
Step9:  plt.plot is the standard do-it-all plotting function in matplotlib. Everything about how the series looks can be modified.
B. Goodness-of-fit
How good was your guess? How do you even answer that question?
<details>
<summary>Answer</summary>
Let's use something called the $L_2$ norm
Step10:  Try changing the parameters to something bad and see what happens to the value of l2_norm. Since this definition of the $L_2$ norm is not normalized by something like a standard deviation of the data, it can't tell us in absolute terms how good our funciton fits, but it can at least tell us if one set of parameters fits better than another. This is really helpful!
Step11:  But, how do we know when we have a best fit? How would you try to figure it out?
Thankfully, we don't have to create our own method to do this. The smart people working on the scipy package have already built an optimized tool for us to use! It's called the curve_fit function as is part of the scipy.optimize sub-package.
C. Fitting
Step12:  scipy.optimize.curve_fit is a type of minimization function. In this case, the function finds the parameters of another given function that minimize the $L_2$ norm between the data points and what our gauss function thinks the data points should be at a given x.
A quick, useful way to see what a function does without having to google it is to use the built-in python help function.
Step13:  That gave us a lot of information about the curve_fit function! As you can see, curve_fit takes a function as its first parameter, and it tells us exactly how to arrange the parameters of that function (thankfully, our gauss function should already have this form). The next two parameters curve_fit takes are xdata and ydata (x_data and y_data as we defined them). The rest are optional and will be talked about briefly at the end of this lecture.
Let's try calling help on the gauss function we defined above.
What do you think the help function does?
<details>
<summary>Answer</summary>
You can see help just returns the function "signature" and the string at the start of the function (called the "docstring").
</details>
Now, let's use curve_fit to fit the data to our function!
Step14:  Do you know what popt is? How would you find out?
<details>
<summary>Answer</summary>
If you look back at the `help` output from `curve_fit`, `popt` is a list of the best-fit parameters of our gauss function for this data. The parameters in the list are in the order that the parameters are listed in our `gauss` function (`mean`, `std`, `amp`). Let's try plotting the data, our guess, and the best fit from `curve_fit`!
</details>
Step15:  How close was your guess to the best fit?
D. Interpreting Fitting Errors
pcov is a little more complicated. pcov is what's called the "covariance matrix" of the best fit parameters. As shown in the help output, the standard deviations of the parameters can be recovered from this matrix in the following way
Step16:  The covariance matrix looks like this for the parameters in our gauss function
\begin{bmatrix}
s_{\mu}^2 & cov(\mu, \sigma) & cov(\mu, A)\
cov(\sigma, \mu) & s_{\sigma}^2 & cov(\sigma, A)\
cov(A, \mu) & cov(A, \sigma) & s_A^2
\end{bmatrix}
Where $s_x^2$ is the variance of a parameter $x$ and $s_x$ is its estimated standard deviation, and $cov(x, y)$ is the covariance between parameters $x$ and $y$.
Can you guess what np.sqrt(np.diag(pcov) does now?
Covariance can be difficult to visualize. It's often much easier to look at something called the "correlation coefficient" instead. The correlation coefficients can be easily found from the covariance matrix by using this transformation
<ASSISTANT_TASK:>
Python Code:
import numpy as np
Explanation: Review from the previous lecture
In yesterday's Lecture 2, you learned how to use the numpy module, how to make your own functions, and how to import and export data. Below is a quick review before we move on to Lecture 3.
Remember, to use the numpy module, first it must be imported:
End of explanation
np.linspace(0,10,11)
Explanation: You can do a lot with the numpy module. Below is an example to jog your memory:
End of explanation
#your code here
# start by defining the length of the array
arrayLength = 10
# make a numpy array of 10 zeros
# let's define the first element of the array
# make your while loop (don't forget to initialize i!) and print out the new array 
Explanation: Do you remember loops? Let's use a while loop to make an array of 10 numbers. Let's have each element be increased by 2 compared with the previous element. Let's also have the first element of the array be 1.
End of explanation
# This data was generated using the methods described in the preview notebook. Feel free to read through the preview to see how it was done.
data_filename = "https://raw.githubusercontent.com/astroumd/GradMap/master/notebooks/Lectures2021/Lecture3/Data/photopeak.txt"
x_data, y_data = np.loadtxt(data_filename, usecols=(0, 1), unpack=True)
print(x_data, y_data)
Explanation: There's your quick review of numpy and a while loop. Now we can move on to the content of Lecture 3.
A. Data Introspection
One of the most common tasks in experimental physics is trying to model experimental data with a function. This lecture will walk through how to accomplish this with python. First, let's import some data.
End of explanation
# Import the matplotlib package. Really, we only need the pyplot sub-package, so we import it this way. Renaming it to 'plt' is common practice.
import matplotlib.pyplot as plt
# The most basic kind of scatter plot you can plot using matplotlib is done like this:
plt.scatter(x_data, y_data)
# At the end of a cell where you are plotting things, this line tells python that you want do display the plots you defined in the cell.
plt.show()
Explanation: Let's plot it to see what it looks like!
In python there are multiple plotting libraries, but the moset common one is matplotlib, and that is the one we will be using today.
End of explanation
# Change the body of this function so that is returns the value at x of the gaussian defined by the parameters.
def gauss(x, mean, std, amp):
    Write what the function does here!
    raise NotImplementedError
Explanation: Usually, you will want to have some theory-based motivation for the function you choose to model some set of data, but for this example, we don't know anything about the data other than the points themselves. In this type of situation, trying to fit a simple function to the data is not a bad first step in trying to understand it. What function do you think might fit this data based on how it looks in the plot?
<details>
<summary>Answer</summary>
It looks like the data is shaped like a normal (gaussian) distribution, so let's try to fit it to that! First, let's define a gaussian function for fitting.
The equation for a Gaussian curve is the following:
$f(x) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{-(x-\mu)^2}{2\sigma^2}}$
where $\mu$ is the mean and $\sigma$ is the standard deviation. We also want to be able to scale our function to fit the scale of the data, so we should multiply the equation by some amplitude, A.
</details>
End of explanation
mu = 
sigma = 
amp = 
Explanation: What values of these parameters do you think will match the data above?
End of explanation
x_points = # Use np.linspace here.
plt.plot(x_data, y_data, linestyle='', marker='.', markersize=12) # This line makes the same plot as plt.scatter, but avoids some quirks in matplotlib.
# Try replacing this with plt.scatter(x_data, y_data) and see what happens!
# plt.scatter(x_data, y_data)
y_gauss = # Finish this line to get the y-values using the function you made above.
plt.plot(x_points, y_gauss)
plt.show()
Explanation: Let's try plotting the function the function with these parameters next to the data! For this, we should define some evenly-spaced x-values to calculate the function at using np.linspace:
End of explanation
y_data_gauss = # Finish this line to get y_gauss at the x_data.
l2_norm = # Finish this line to calculate the L_2 norm.
print(l2_norm)
Explanation: plt.plot is the standard do-it-all plotting function in matplotlib. Everything about how the series looks can be modified.
B. Goodness-of-fit
How good was your guess? How do you even answer that question?
<details>
<summary>Answer</summary>
Let's use something called the $L_2$ norm: $||y - f(x)||^2$ to get a metric of the difference between our data and our function. This may sound and look fancy, but all it's doing is calculating the distance at each x-value in the data between the data y-value and the function y-value. Then, it squares those distances and adds them all together.
This is defined by `np.sum((y_data - y_data_gauss)**2)` for our setup.
</details>
End of explanation
mu = 
sigma = 
amp = 
l2_norm = np.sum((y_data - gauss(x_data, mu, sigma, amp))**2)
print(l2_norm)
Explanation: Try changing the parameters to something bad and see what happens to the value of l2_norm. Since this definition of the $L_2$ norm is not normalized by something like a standard deviation of the data, it can't tell us in absolute terms how good our funciton fits, but it can at least tell us if one set of parameters fits better than another. This is really helpful!
End of explanation
import scipy.optimize
Explanation: But, how do we know when we have a best fit? How would you try to figure it out?
Thankfully, we don't have to create our own method to do this. The smart people working on the scipy package have already built an optimized tool for us to use! It's called the curve_fit function as is part of the scipy.optimize sub-package.
C. Fitting
End of explanation
help(scipy.optimize.curve_fit)
Explanation: scipy.optimize.curve_fit is a type of minimization function. In this case, the function finds the parameters of another given function that minimize the $L_2$ norm between the data points and what our gauss function thinks the data points should be at a given x.
A quick, useful way to see what a function does without having to google it is to use the built-in python help function.
End of explanation
# The variable names popt and pcov come from the curve_fit function. We will get into what they mean soon!
popt, pcov = scipy.optimize.curve_fit(gauss, x_data, y_data)
print(popt)
print(pcov)
Explanation: That gave us a lot of information about the curve_fit function! As you can see, curve_fit takes a function as its first parameter, and it tells us exactly how to arrange the parameters of that function (thankfully, our gauss function should already have this form). The next two parameters curve_fit takes are xdata and ydata (x_data and y_data as we defined them). The rest are optional and will be talked about briefly at the end of this lecture.
Let's try calling help on the gauss function we defined above.
What do you think the help function does?
<details>
<summary>Answer</summary>
You can see help just returns the function "signature" and the string at the start of the function (called the "docstring").
</details>
Now, let's use curve_fit to fit the data to our function!
End of explanation
x_points = 
plt.plot(x_data, y_data, linestyle='', marker='.', markersize=12, label='data')  # The parameter "label" gives a name to the data series.
# Try adding labels to these plots!
plt.plot() # Finish this line to plot your hand-picked best-guess parameters from earlier.
plt.plot()  # Finish this line to plot the best fit parameters from curve_fit.
# This function produces a legend of the plots and their labels. You can leave this line as is.
plt.legend()
plt.show()
Explanation: Do you know what popt is? How would you find out?
<details>
<summary>Answer</summary>
If you look back at the `help` output from `curve_fit`, `popt` is a list of the best-fit parameters of our gauss function for this data. The parameters in the list are in the order that the parameters are listed in our `gauss` function (`mean`, `std`, `amp`). Let's try plotting the data, our guess, and the best fit from `curve_fit`!
</details>
End of explanation
perr = np.sqrt(np.diag(pcov)) # This comes straight from the curve_fit help output.
# This is an example of string formatting in python.
# Each set of {} corresponds to a parameter passed to .format().
# .3e means "format this number in scientific notation with 3 digits after the decimal".
# using "*perr" is a nice python trick that expands perr into the three individual values it contains.
# All of this is out of the scope of this lecture, but it's good to get exposure to these things.
# Just look at how nice the out put looks!
print('s_mu = {:.3e}, s_sigma = {:.3e}, s_amp = {:.3e}'.format(*perr))
Explanation: How close was your guess to the best fit?
D. Interpreting Fitting Errors
pcov is a little more complicated. pcov is what's called the "covariance matrix" of the best fit parameters. As shown in the help output, the standard deviations of the parameters can be recovered from this matrix in the following way: perr = np.sqrt(np.diag(pcov))
Since we aren't teaching linear algebra here, all of the matrix manipulations will be given to you. All that needs to be taken away from this is how to read this specific matrix.
End of explanation
perr_transpose = np.atleast_2d(perr).T
pcor = pcov / perr / perr_transpose
print(pcor)
Explanation: The covariance matrix looks like this for the parameters in our gauss function
\begin{bmatrix}
s_{\mu}^2 & cov(\mu, \sigma) & cov(\mu, A)\
cov(\sigma, \mu) & s_{\sigma}^2 & cov(\sigma, A)\
cov(A, \mu) & cov(A, \sigma) & s_A^2
\end{bmatrix}
Where $s_x^2$ is the variance of a parameter $x$ and $s_x$ is its estimated standard deviation, and $cov(x, y)$ is the covariance between parameters $x$ and $y$.
Can you guess what np.sqrt(np.diag(pcov) does now?
Covariance can be difficult to visualize. It's often much easier to look at something called the "correlation coefficient" instead. The correlation coefficients can be easily found from the covariance matrix by using this transformation:
$cor(x, y) = \frac{cov(x, y)}{s_x s_y}$
Performing this transformation on the covariance matrix gives the correlation matrix, which looks like this:
\begin{bmatrix}
1 & cor(\mu, \sigma) & cor(\mu, A)\
cor(\sigma, \mu) & 1 & cor(\sigma, A)\
cor(A, \mu) & cor(A, \sigma) & 1
\end{bmatrix}
The code for this is a bit beyond the scope of this lecture, so it has been done for you below.
End of explanation
<END_TASK>