bias_vs_variance_of_estimation.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit


def model1(x, k):
    return k * x

def model2(x, k, a):
    return k * x + x * np.sin(a * x)

x = np.linspace(0, 5, 30)
k = 1.5
a = 1.
y = k * x + x * np.sin(a * x)
yerr = 0.3

lin = list()
linsin = list()
for i in range(300):
    data = y + np.random.normal(0, yerr, 30)
    popt1, pcov1 = curve_fit(model1, x, data, sigma=np.repeat(yerr, 30))
    lin.append(popt1[0])
    popt2, pcov2 = curve_fit(model2, x, data, sigma=np.repeat(yerr, 30))
    linsin.append(popt2)


plt.figure()
plt.plot(x, data, '.k', label='data')
plt.xlabel(r'x')
plt.ylabel(r'y')
xx = np.linspace(0, 5, 1000)
yy = model2(xx, k, a)
plt.plot(xx, yy, label='true')
plt.legend()
for i in range(0, 300, 15):
    plt.plot(xx, model1(xx, lin[i]), alpha=0.15, color='r')
for i in range(0, 300, 15):
    plt.plot(xx, model2(xx, linsin[i][0], linsin[i][1]), alpha=0.15, color='g')