Merge pull request #340 from pep-dortmund/data_to_fitting_exercises

add polyfit exercise and data to uncertainties curve_fit wrapper

exercises-toolbox/4-scipy/3-polyfit/aufgabe.txt

@@ -0,0 +1,4 @@
+Führe einmal `python loesung.py` aus, um die Daten zu erzeugen.
+
+Fitte die Daten aus daten.txt (x, y, y_err) an ein Polynom zweiten Grades.
+Nutze hierfür np.polyfit.

exercises-toolbox/4-scipy/3-polyfit/loesung.py

@@ -0,0 +1,40 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Generate data
+rng = np.random.default_rng(210)
+N = 50
+data_x = np.linspace(-10, 10, N)
+error_y = rng.lognormal(-1, 0.2, size=N)
+data_y = rng.normal(data_x**2, error_y)
+np.savetxt("daten.txt", np.column_stack([data_x, data_y, error_y]), header="x y y_err")
+
+# The solution starts here
+x, y, e_y = np.genfromtxt("daten.txt", unpack=True)
+
+
+def f(x, a, b, c):
+    return a * (x + b)**2 + c
+
+
+parameters, covariance_matrix = np.polyfit(x, y, deg=2, cov=True)
+uncertainties = np.sqrt(np.diag(covariance_matrix))
+
+for name, value, unc in zip('abc', parameters, uncertainties):
+    print(f'{name} = {value:.3f} ± {unc:.3f}')
+
+fig = plt.figure(layout="constrained")
+ax = fig.add_subplot()
+
+ax.errorbar(x, y, yerr=e_y, fmt="k.", label="data")
+
+t = np.linspace(-10, 10, 500)
+ax.plot(t, f(t, *parameters), label="Fit")
+ax.plot(t, t**2, "--", label="Original")
+
+ax.set_xlim(t[0], t[-1])
+ax.set_xlabel(r"$t$")
+ax.set_ylabel(r"$f(t)$")
+ax.legend()
+
+plt.savefig("loesung.pdf")

exercises-toolbox/4-scipy/3-curve_fit/aufgabe.txt → exercises-toolbox/4-scipy/4-curve_fit/aufgabe.txt

@@ -1,8 +1,7 @@
-# curve_fit
+# curve_fit 
 
 Aufgabe:
-
-Führe einmal `python loesung.py` aus, um die daten zu erzeugen.
+Führe einmal `python loesung.py` aus, um die Daten zu erzeugen.
 
 Fitte die Daten aus daten.txt (x, y, y_err) mit der Funktion f(x) = a * sin(b * x + c) + d.
 Benutze dazu scipy.optimize.curve_fit.

exercises-toolbox/4-scipy/3-curve_fit/vorlage.py → exercises-toolbox/4-scipy/4-curve_fit/vorlage.py

@@ -19,5 +19,5 @@ def f(x, a, b, c, d):
 ax.set_xlim(t[0], t[-1])
 ax.set_xlabel(r"$t$")
 ax.set_ylabel(r"$f(t)$")
-ax.legend(loc="best")
+ax.legend()
 fig.savefig("loesung.pdf")

exercises-toolbox/4-scipy/4-peakdetect/loesung.py → exercises-toolbox/4-scipy/5-peakdetect/loesung.py

@@ -35,5 +35,5 @@ def e(x, a, b, c):
 ax.set_xlabel(r"$t \ / \ \mathrm{ms}$")
 ax.set_ylabel(r"$U \ / \ \mathrm{V}$")
 ax.set_xlim(0, 0.3)
-ax.legend(loc="best")
+ax.legend()
 fig.savefig("loesung.pdf")

exercises-toolbox/4-scipy/4-peakdetect/vorlage.py → exercises-toolbox/4-scipy/5-peakdetect/vorlage.py

@@ -21,5 +21,5 @@
 ax.set_xlabel(r"$t \ / \ \mathrm{ms}$")
 ax.set_ylabel(r"$U \ / \ \mathrm{V}$")
 ax.set_xlim(0, 0.3)
-ax.legend(loc="best")
+ax.legend()
 fig.savefig("loesung.pdf")

exercises-toolbox/5-uncertainties/3-curve_fit/aufgabe.txt

@@ -21,3 +21,7 @@ def ucurve_fit(f, x, y, **kwargs):
     …
     … = scipy.optimize.curve_fit(…, **kwargs)
     …
+
+Überprüfe deine Wrapper-Implementation an den Daten in daten.txt und fitte diese an
+A * cos(B * x) + C.
+Plotte den Fit zusammen mit den originalen Daten.

exercises-toolbox/5-uncertainties/3-curve_fit/loesung.py

@@ -1,4 +1,5 @@
 import numpy as np
+import matplotlib.pyplot as plt
 import scipy.optimize
 import uncertainties as unc
 import uncertainties.unumpy as unp
@@ -15,3 +16,34 @@ def ucurve_fit(f, x, y, **kwargs):
     )
 
     return unc.correlated_values(popt, pcov)
+
+
+def f(x, a, b, c):
+    return a * np.cos(x * b) + c
+
+# Generate data
+length = 100
+x = np.linspace(0, 3 * np.pi, length)
+rng = np.random.default_rng()
+y1 = rng.normal(0.0, 0.2, length)
+y2 = np.abs(rng.normal(0.0, 0.2, length))
+
+y = f(x, 1, 1, -3) + y1
+
+np.savetxt("daten.txt", np.column_stack([x, y, y2]), header="x\ty\ty_err")
+
+# Solution
+x, y_0, y_err = np.genfromtxt("daten.txt", unpack=True)
+y = unp.uarray(y_0, y_err)
+params = ucurve_fit(f, x, y)
+print("a * cos(x * b) + c")
+for char, p in zip("abc", params):
+    print(f"{char} = {p}")
+
+fig = plt.figure(layout="constrained")
+ax = fig.add_subplot()
+ax.errorbar(x, unp.nominal_values(y), yerr=y_err, fmt=".", label="Daten")
+ax.plot(x, f(x, *unp.nominal_values(params)), label="Fit")
+ax.set_xticks([0, np.pi, 2 * np.pi, 3 * np.pi], [0, "π", "2π", "3π"])
+ax.legend()
+plt.savefig("loesung.pdf")