<jupyter_start><jupyter_text>Fitting curves to point clouds using numpy and scipy     Note for contributors: Remember to run Kernel > Restart & Clear output before adding any changes to git! **In this tutorial, we make our lives a bit easier and use some pre-defined fitting functions.**<jupyter_code># Numerical operations:
import numpy as np

# Plotting library:
import matplotlib.pyplot as plt

# Minimizing/fitting library:
import scipy.optimize<jupyter_output><empty_output><jupyter_text>Remember: You can always get help about any function, e.g.     np.array by typing``np.array?``  Fitting lines or polynomial Fitting a polynomial to a point cloud.For this we use the ``polyfit`` function of numpy (https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html). Source: https://stackoverflow.com/questions/19165259/ Let's start by defining some points and their x and y projections:<jupyter_code>points = np.array([(1, 1), (2, 4), (3, 1), (9, 3), (10, 10)])
x = points[:,0]
y = points[:,1]<jupyter_output><empty_output><jupyter_text>Plot them: Exercise 1 [super-easy]: Plot the points Fitting a polynomial of degree ``deg`` is easy as calling:<jupyter_code>coeffs = np.polyfit(x, y, deg=3)<jupyter_output><empty_output><jupyter_text>The vector ``coeffs`` now contains the coefficients of the polynomial as a vector of length ``deg+1``:<jupyter_code>coeffs<jupyter_output><empty_output><jupyter_text>It's easy to also get the corresponding function:<jupyter_code>f = np.poly1d(coeffs)<jupyter_output><empty_output><jupyter_text>Exercise 2 [super-easy]: Evaluate the difference $f(x)-y$ and interpret them. Let's finally have a look at the fit:<jupyter_code># Some points where we evaluate our new function
x_new = np.linspace(x[0], x[-1], 50)
y_new = f(x_new)

# Plot the datapoints
plt.plot(x, y, 'ko', label="data")

# Plot our fitted polynomial
plt.plot(x_new, y_new, 'r-', label="fit")

# Add legend etc
plt.legend()
plt.xlim([x[0]-1, x[-1] + 1 ])
plt.show()<jupyter_output><empty_output><jupyter_text>Exercise 3 [easy]: Fit a straight line to the following dataset and determine the line parameters. Bonus points: Also plot the fit as above.<jupyter_code>x_exercise3 = np.linspace(0, 1, 20)
y_exercise3 = x_exercise3 + np.random.random(len(x_exercise3))
plt.plot(x_exercise3, y_exercise3, ".")<jupyter_output><empty_output><jupyter_text>Fitting Gaussian (or any arbitrary function) From https://stackoverflow.com/questions/19206332 Again let's generate some points:<jupyter_code>x = np.asarray(range(10))
y = 0.3*x + np.asarray([0,1,2,3,4,5,4,3,2,1])<jupyter_output><empty_output><jupyter_text>Here, we define our own function that we want to fit to our data points:<jupyter_code>def gaus(x, norm, mean, sigma):
    # Note that this function takes a whole vector x of data!
    return norm * np.exp(-(x-mean)**2/(2*sigma**2))<jupyter_output><empty_output><jupyter_text>This function has 3 free parameters, ``norm``, ``mean``, ``sigma``, which we will now fit using ``scipy.optimize.curve_fit`` (https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html)<jupyter_code>popt, pcov = scipy.optimize.curve_fit(gaus, x, y)<jupyter_output><empty_output><jupyter_text>The variable ``popt`` holds the parameters:<jupyter_code>popt<jupyter_output><empty_output><jupyter_text>The other variable, ``pcov`` holds the covariance matrix of the fitted values. This gets relevant once you want to give errors on your fitted quantities, but we'll ignore this for now.<jupyter_code># Plot points
plt.plot(x, y,'ko',label='data')

# Plot our gaussian
# Define some points on the xaxis
x_fine = np.linspace(min(x), max(x), 200)
y_fine = gaus(x_fine, *popt)
# The *  unpacks the values from popt and uses them as
# parameters.

plt.plot(x_fine, y_fine, 'r-', label='fit')

# ...
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Exercise 4a [easy]: Modify the above example so that the sigma above is     fixed to 1.    Exercise 4b [medium]: Let's try to fit the same dataset with a slightly more complicated model that consists of a linear part plus a Gaussian, i.e.     ```f(x) = x + b + gaus(x, norm, mean, sigma)```Hint:Create a new function ``mymodel`` which takes the parameters ``x``, as well as all fitted parameters, i.e. ``b``, ``norm``, ``mean``, ``sigma``. In the definition you can also use the function ``gaus`` from above!<jupyter_code>def mymodel(x, b, norm, mean, sigma):
    return x + b + gaus(x, norm, mean, sigma)<jupyter_output><empty_output><jupyter_text>Exercise 4c [hard]: Fit the above points as two lines, i.e. $$f(x) = \begin{cases}  a_1 x + b_1 & \text{if}\ x<c\\  a_2 x + b_2 & \text{if}\ x\geq c\\\end{cases}$$    Such that $a_1 c + b_1 = a_2 c + b_2$.Hint: Solve the contraint for, say, b2Define a function line(x, a1, b1, a2, c)Note that x is a vector! Thus if x&lt;c: ... won't work. Rather, take a look at the function np.where.<jupyter_code>def line(x, a1, b1, a2, c):
    return np.where(x<c, a1*x+b1, a2*x+a1*c+b1-a2*c)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13541459534698747, 'fraction_numerical': 0.01988466892026248, 'mean_word_length': 5.558018252933508, 'pattern_counts': {'":': 0, '<': 59, '<?xml version=': 0, '>': 58, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '517149', 'license': 'MIT', 'n_tokens_mistral': 1881, 'n_tokens_neox': 1783, 'n_words': 549, 'path': 'examples/jupyter_notebooks/002_fitting_curves_to_points.ipynb', 'repo_name': 'klieret/HEPFittingTutorial', 'size': '5029'}