# SciPy: Python tools for science
*Monica Rizzo, 2019*

Among the many useful python packages for data analysis, SciPy is one of the most commonly used and widely applicable. In this notebook we're going to be going over some of the functions in Scipy you can use to fit and model data. 

First things first, lets import all of the important stuff.


```python
import numpy as np
from math import factorial
from scipy.optimize import curve_fit
from scipy.interpolate import interp1d
from scipy.optimize import brentq

from scipy.integrate import quad
from scipy.misc import derivative


import matplotlib.pyplot as plt


#this makes the plot labels bigger (looks nicer)
plt.rcParams.update({'font.size': 18})
```

## 1. Regression

There are a lot of things you can do with data. Often times, we will want to either 1) see how well data agrees with a given model or 2) without assuming anything, see what sort of model agrees the best with some data.
One way to test this is by using regression techniques. "Regression" is just the math word for fitting data to a certain function. You may have done something similar to this before with Excel, or by hand using graph paper. 

Lucky for us, SciPy can fit things for you! First, lets generate some data that follows a known function.
1. **Using numpy, create an array of n equally spaced x values within some range - doesn't really matter what range**
2. **Then choose two values for the slope m and y intercept b**
3. **Finally, create a numpy array of y values which follows the equation for a line (*hint: take advantage of numpy array addition/multiplication*)**


```python
#number of data points to generate
n = #complete

#create a numpy array of x values   
x_values = #complete

#define slope and y intercept
m = #complete
b = #complete

#make an array of y values that follow the equation of a line
y_values = #complete
```

Now let's plot the data to see what it looks like (always a good first step).
1. **Make a scatter plot of the x and y values you created above**


```python
plt.figure(figsize=(10, 8))
#complete
plt.grid(True)
#ALWAYS LABEL YOUR AXES
plt.xlabel(#complete
plt.ylabel(#complete
plt.show()
```

Obviously, we already know what this function looks like, so we already know what Scipy should return as the answer. 

The Scipy function we are going to use to do the fitting is called *curve_fit*. Using a user-defined python function that describes the fit (in this case a line), *curve_fit* will attempt to find the best fitting parameters for some input data. So, in this case, we want to create a fit function that returns a line and takes 2 free parameters: the slope and the y intercept. Then, we want to hand *curve_fit* that fit function, the x data, and the y data, and it's going to tell us what it thinks the best fitting parameters are, along with the uncertainties on those parameters (which we're not going to worry about right now).

1. **Define a function which takes as inputs x, m, and b, and returns the y values corresponding to a line**
2. **Pass *curve_fit* the name of your fit function, your x data, and your y data (see documentation here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html)**
3. **Print the result of *curve_fit* (note: result should be a 2 element array. The values are in the same order as the arguments of your fit function)**


```python
def fit_function():
    #complete
    
```


```python
result, cov = curve_fit(#complete

#print the result

```

Nice.



Now let's look at a more realistic example. In real life, data is noisy and usually has some associated uncertainty. Lets add a little bit of noise to our existing data, and assign some uncertainty.

1. **Add the noise values to your array of y values**
2. **Create an array - which will be your y error - of length n where all of the values are 3**


```python
#numpy has built-in capabilities to generate random numbers :)
#this will generate random numbers drawn from a bell curve
y_noise = np.random.normal(0, 3, n)

#array of y values plus noise
y_noisy_values = #complete

#array uncertainty values
y_error = #complete
```

Now, let's plot the data again. 

***Hint: in order to plot the errors, instead of using the 'plot' function, use 'errorbar': https://matplotlib.org/3.1.0/api/_as_gen/matplotlib.pyplot.errorbar.html***


```python
plt.figure(figsize=(10, 8))
#try using fmt='o' and markersize=6.0
plt.errorbar(#complete
#complete
```

Not so pretty anymore, but maybe still looks like a line? Let's see what SciPy thinks. 
1. **Using the same syntax as before, pass *curve_fit* your fit function, your x data, your noisy y data, and also the uncertainty on your y values - this can be done using the 'sigma' option.**
2. **As before, print the result**


```python
#use curve fit to find the fit parameters for your noisy data
result, cov = curve_fit(#complete

#print the results
```

As you can see, this answer is pretty different from the actual values of the parameters. We can plot what the fit looks like over our data to visually assess how bad/good it is. 
1. **Create a numpy array of *n_test* x values to plot, that has a range that is slightly bigger that the range of your x values**
2. **Pass this array, along with the best fit m and b values, to your fit function, and store the result as a separate array**
3. **Copy the cell you used to plot the data with errorbars, and add a line to plot the test x and y values**


```python
#number of test values
n_test = 500

#array of x test values
x_test_values = #complete

#array of y test values 
y_test_values = fit_function(#complete
```


```python
plt.figure(figsize=(10, 8))
#complete
```

Keep in mind that these answers are just guesses - they also have some uncertainty to them. The other thing that *curve_fit* returns - which we're calling *cov* - describes the uncertainites on the fit parameters. Let's see what those look like.


```python
#cov is a matrix, and we only want the diagonal values
unc = np.sqrt(np.diag(cov))

print("param1_err = %0.2f and param2_err = %0.2f" % (unc[0], unc[1]))
```

**Given what your data looks like, do you think these numbers are reasonable? Do your fit parameters to the noisy data agree with the real parameters within this uncertainty?**

## 2. Interpolation

Now let's switch gears a little bit. Suppose we have some data where instead of the function parameters (slopes, y-intercepts, etc.), we want to approximate what the values between our data points are. So in this case, we don't care about the specifics of the fit, we just want to connect the dots. This can be achieved using *interpolation* techniques. 

First, lets again create some data that follows a known function:
1. **Using numpy arrays, create an array of 10 x values that ranges from 0 to 2$\pi$**
2. **Create an array of y values which is the sine of the x values**


```python
#array of x values
x_values = #complete

#array of y values
y_values = #complete
```

As per *always*, plot your data before doing anything else.


```python
plt.figure(figsize=(10, 8))
#complete
```

This could be sine curve - and I'm sure if you connected the dots by hand you would draw something that resembles a sine function. But let's see how well we can approximate this by numerically interpolating. To start, let's use the default settings of *interp_1d*. When passed x and y values, *interp_1d* will return a function which interpolates the data at every point between the maximum and minimum x value. We can then plot the values of this interpolation function over a range of test x values. 

1. **Pass *interp1d* your x and y values (see documentation here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html#scipy.interpolate.interp1d)**
2. **Create an array of *n_test* test x values between 0 and 2$\pi$**
3. **Pass your interpolation function your array of test x values and store the result as test y values**
4. **Plot your test values as a line over a scatter plot of your data**


```python
#create an interpolation function
interp_func = interp1d(#complete

#array of test x values
x_test_values = #complete

#recall that interp_func is a function, so it can be used like any other function
y_test_values = #complete
```


```python
plt.figure(figsize=(10, 8))
#complete
```

As you can see, *interp_1d* connects the data points with straight lines by default. This isn't always the best way to interpolate data; we can tell just by looking at what we've plotted that our data most likely has some curvature to it. It's also generally better to interpolate data smoothly, because this makes performing operations like differentiation easier/more accurate. Lets try out some different types of interpolation.

To change the way *interp_1d* connects our data points, we can use the 'kind' option.
1. **Try setting kind to each of the following: 'zero', 'quadratic', 'cubic'. For each different kind, create a new variable name to store the interpolating function in.**
2. **Using the same method as above, plot the data and the interpolated approximation on top of each other. You can use the same test x values, but will need to create new test y values for each plot**


```python
interp_func_zero = interp1d(x_values, y_values, kind='zero')
#repeat for all interpolation kinds

```


```python
#test y values

#zero plot
```


```python
#test y values

#quadratic plot
```


```python
#test y values

#cubic plot
```

How can we judge which of these kinds of interpolation is the best? We know what the values of the function *should* be, so we can just pick a point and see how different the approximate answer is from the real answer. In the cell below, check to see how different each of the approximate y values is from the real value where x = 4. 

1. **Compute the real value - sin(4) - and store in a variable**
2. **For each of the interpolation functions, compute and store in different variables the approximate y value at an x value of 4**
3. **Calculate and print the percent error of each of the approximate y values**

*Hint: calculate percent error using: \begin{equation}\frac{|approximation - real|}{|real|} \times 100\end{equation}*

*Hint: python has a built in absolute value function called 'abs'*


```python
real_value = np.sin(4)

approx_value_linear = interp_func(4)
#repeat for all interpolation functions

#calculate percent error for each interpolation function
```

Based on this, which appears to be the best kind of interpolation? This happens to be the best kind in general. 

*Optional*: We can also check to make sure the area under the curve is what we expect it to be. For a sine curve in the range zero to 2$\pi$, the total area should be zero - this is because half the the curve is above zero and half is below. In order to find the area under the curve, we can *integrate* it - SciPy also has functions to do this. The particular algorithm we're going to use is called *quad* (which is all around the most accurate in SciPy). *quad* takes as arguments the name of the function you want to integrate, the upper bound of the integration, and the lower bound of the integration. 

1. **For each of the interpolation functions, use quad to integrate from 0 to 2$\pi$ (see documentation here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html#scipy.integrate.quad)**
2. **Print out the result of each. Why do you think the results are what they are?**

*Note: when dealing with computers, numbers smaller than or close to $10^{-16}$ are considered zero*


```python
res_linear, err = quad(interp_func, 0, 2*np.pi)
#repeat for all interpolation functions

```

## 3. Root Finding 

Now let's talk - briefly - about something completely different. In math classes, you may have needed to find the zeros of functions before (set y equal to zero and solve for x). Finding zeros is very important in the area of numerical optimization - which we won't talk about explicity - but can often be more complicated than just setting y equal to zero and solving for x. In fact, when you have just data points, you often can't do this at all. 

And that's where SciPy comes in. Suppose we have some interpolating function which passes through zero, and we want to know approximately where. One of the best algorithms to do this is something called *Brent's method* (which you may or may not learn about in college). *Brent's method* will attempt to find where a function crosses zero in a given interval; it needs to know the endpoints to look between in order to work. 

First, lets create some data to try this out on:
1. **Using numpy, create an array of 10 equally spaced x values between -2 and 2**
2. **Create an array of y values that is the x values cubed plus 3**


```python
#array of x values
x_values = #complete

#array of y values
y_values = #complete
```

Plot your data.


```python
plt.figure(figsize=(10, 8))
#complete
```

Based on the plot, at what x value do you expect the zero point to be? 

Let's see what SciPy thinks. We're going to use a function called *brentq*. This function takes the name of the function you want the zero of, the upper bound of that function, and the lower bound of that function. So first, we need to interpolate the data, and then pass that interpolant to *brentq*.

1. **Using the best interpolation kind from the previous section, create and store an interpolating function for your x and y data**
2. **Pass *brentq* the name of your interpolating function, the lower bound of the function, and the upper bound of the function (see documentation here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.html#scipy.optimize.brentq)**
3. **Print the result. Does this answer seem correct?**


```python
#interpolating function
interp_func = interp1d(#complete

#find zero using brentq
x0 = brentq(#complete

#print result
    
```

## 4. Fun With Real Data

Now let's tackle a (sort of) real-life example. Using the provided data file, we will attempt to fit the data it contains with a more complicated fit function. First, we will load and plot the data.

1. **Load the data and plot it (with error bars)**


```python
#load data and store it to an array
#columns: x values, y values, y error
data = #complete

#split the data into separate x, y, and y_err arrays using numpy array slicing
x_values = #complete
y_values = #complete
y_error = #complete
```


```python
#plot. that. data
plt.figure(figsize=(10, 8))
#complete
```

What do you think this function looks like? Using your knowledge of regression, try to think of a fitting function which might model this data, and then fit that model to the data.

1. **Define a fitting function, and fit it to the data using *curve_fit***
2. **Print the best fit parameters and their uncertainties**
3. **Plot the data and the fit on top of each other**

*Hint: try to fit a higher order polynomial*


```python
#create a fitting function
def :

#use curve fit to fit the function to the data
result, cov = #complete

#print the result and the uncertainties
#you can probably just copy and paste code from the first section, with minor adjustments

```


```python
#create test x values, and test y values from the fit

#plot the data and the fit 
plt.figure(figsize=(10, 8))

#complete

plt.show()
```

This data actually has a function which fits it best, for science reasons. In order to fit that function we're going to cheat a little and find the location of the peak, i.e. where the slope of the function equals zero. To do this, we need to first interpolate our data, and then find where the slope (or *derivative*) equals zero. 

1. **Interpolate the data using the best kind of interpolation**
2. **Define a function that takes x values and returns the derivative of your interpolation function at those x values (already completed for you)**
3. **Pass the derivative function to *brentq* (remember you need to pass it a lower and upper bound), and store the result. Print this result. Does it make sense?**


```python
#interpolation function 
interp_func = interp1d(#complete

#function which takes x values and returns derivative
def deriv_func(x):
    
    deriv = derivative(interp_func, x, dx = 1e-6)
    
    return deriv

#find root of derivative function 
x0 = brentq(#complete

#print result
```

The function we want to fit to is:

\begin{equation}f(x) = a ((x/x0)^{b1/-0.8} + (x/x0)^{b2/-0.8})^{-0.8},\end{equation}

where a, b1, and b2 are the three fit parameters. Now that we know x0, we can use *curve_fit* to find the remaining 3 parameters.

1. **Define a function called 'broken_power_law' which takes x values, a, b1, and b2 as input and returns the values corresponding to the above expression**
2. **Use curve fit to fit this function to the data, *including the error***
3. **Plot the data and the fit on top of each other. Would you say this fit is good or bad? Why or why not?**


```python
#create a function to return the expression above
def broken_power_law:
    #complete
```


```python
#fit the above function using curve fit
result, cov = curve_fit(broken_power_law, x_values, y_values, sigma=y_error)
```


```python
#create test x and y values 

#plot the data with error and the fit
plt.figure(figsize=(10, 8))
#complete
```

## More Fun Stuff: Bonus

As far as cool and useful tools go, SciPy has a *ton* of other nifty functions (with decent documentation), for almost every application. The data above is from a course I took on radiative processes in astrophysics (which is a complicated way of saying "light curves: the class"). You can learn various things about the object which produced this light curve by finding the slope of each side of the logarithm of the data. This is just one idea of something you can do with this data besides just fitting it.   

See what else you can do with it!

Some ideas:

* try a different fitting function
* interpolate the data
* find the point where the slope of the data is zero
* integrate/differentiate the data


```python
##play around with the data and see what you can do with it
##if you need help or suggestions, don't hesitate to ask!









```

## Resources 

A set of SciPy tutorials and documentation can be found here: https://docs.scipy.org/doc/scipy/reference/index.html

The Scipy Lectures site also has a pretty comprehensive guide to all of the modules: https://scipy-lectures.org/intro/scipy.html

Here is a nifty explanation of the difference between regression and interpolation (which I found useful): https://stats.stackexchange.com/questions/33659/how-is-interpolation-related-to-the-concept-of-regression

Here is a lecture (which is honestly pretty heavy, but if you want to look through it you can) on numerical root finding methods: http://butler.cc.tut.fi/~piche/numa/lecture0506.pdf



```python

```

026f3982f4f74cf9b655b3f821835fb29f9c08cb
ipynb
Jupyter Notebook
DataAnalysisModule/Intro_to_Scipy.ipynb
psheehan/CIERA-HS-Program
76f7f0ff994e74e646fa34bbb41c314bf7526e9b
Naumen Condor-1.1 MS-PL
2019-06-25T02:36:49.000Z
2020-06-09T21:44:41.000Z
DataAnalysisModule/Intro_to_Scipy.ipynb
psheehan/CIERA-HS-Program
76f7f0ff994e74e646fa34bbb41c314bf7526e9b
Naumen Condor-1.1 MS-PL
DataAnalysisModule/Intro_to_Scipy.ipynb
psheehan/CIERA-HS-Program
76f7f0ff994e74e646fa34bbb41c314bf7526e9b
Naumen Condor-1.1 MS-PL
2019-06-25T15:33:10.000Z
2021-05-12T18:04:36.000Z
Qwen/Qwen-72B
1. YES
2. YES