x = np.linspace(1, 1000, 500)
signal = (
    gaussian(x, 6, 180, 5)
    + gaussian(x, 8, 350, 10)
    + gaussian(x, 15, 400, 8)
    + gaussian(x, 13, 700, 12)
    + gaussian(x, 9, 800, 10)
)
baseline = 5 + 15 * np.exp(-x / 400)
noise = np.random.default_rng(0).normal(0, 0.2, x.size)
y = signal + baseline + noise
y2 = signal * 1.3 + baseline * 2 + noise
y3 = signal * 0.5 + baseline * 0.5 + noise * 3
y4 = signal + baseline * 3 - 2 + noise * 2
total_y = (y, y2, y3, y4)

lam = 1e5
fit_baselines = baseline_fitter.collab_pls(total_y, method='mpls', method_kwargs={'lam': lam})[0]

fig, ax = plt.subplots(tight_layout={'pad': 0.2})
for y_values in total_y:
    data_handle = ax.plot(y_values, 'C0')  # C0 is first item in color cycle
for baseline in fit_baselines:
    baseline_handle = ax.plot(baseline, 'g--')
for y_values in total_y:
    individual_fit_handle = ax.plot(baseline_fitter.mpls(y_values, lam=lam)[0], 'r:')
ax.set_yticks([])
ax.set_xticks([])
ax.legend(
    (data_handle[0], baseline_handle[0], individual_fit_handle[0]),
    ('data', 'collab_pls fits', 'individual fits'), frameon=False
)
plt.show()