fourier series: more consistent with standard definition

[mathause]

• Closes #xxx
• Tests added
• Fully documented, including CHANGELOG.rst

Make our fourier series more consistent with the standard usage (e.g. on wikipedia for the fourier series or in numpy for fft¹), defined as:

$$s_N(x) = A_0 + \sum_{n=1}^N \left(A_n \cos \left(2 \pi \tfrac{n}{P} x \right) + B_n \sin \left(2 \pi \tfrac{n}{P} x \right) \right)$$

So this PR does two things:

1. switch the order of $\cos$ and $\sin$ so that $A_n$ is associated with $\cos$
2. use fractional years from $0...\frac{11}{12}$ instead of $\frac{1}{12}...1$. Note that this is nothing more than a phase shift with an angle of -30°. This phase shift leads to different values for $A_n$ and $B_n$ but to the same generated series.

¹ well fft uses the exponential form but that's equivalent

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 49.76%. Comparing base (0e15d4e) to head (f45bdb1).
Report is 53 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main     #512   +/-   ##
=======================================
  Coverage   49.76%   49.76%           
=======================================
  Files          50       50           
  Lines        3563     3563           
=======================================
  Hits         1773     1773           
  Misses       1790     1790           

Flag	Coverage Δ
unittests	49.76% <100.00%> (ø)

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[mathause]

To mention the motivation behind this - if we disregard the covariates we can use fft to estimate the coefficients with is much faster. With this change the coeffs are much easier to compare:

import numpy as np

import mesmer

# create test data
n_years = 5
n_months = n_years * 12

months = np.tile(np.arange(12), n_years)
alpha = 2 * np.pi * months / 12

monthly_target = 3.14 * np.cos(alpha) + 2.72 * np.sin(alpha) + np.random.randn(n_months)

# use mesmer to estimate coeffs w/o covariates
yearly_predictor = np.zeros(n_months)

coeffs, mse = mesmer.stats._harmonic_model._fit_fourier_coeffs_np(
    yearly_predictor,
    monthly_target,
    first_guess=np.zeros(4)
)
print(coeffs)

# use fft to estimate coeffs
freq = np.fft.rfftfreq(n_years * 12, 1/12)
sel = freq == 1.

coeffs_rfft = np.fft.rfft(monthly_target) / n_months

a_1 = 2 * coeffs_rfft[sel].real
b_1 = -2 * coeffs_rfft[sel].imag

# compare them
np.testing.assert_allclose(coeffs[1], a_1, rtol=1e-5)
np.testing.assert_allclose(coeffs[3], b_1, rtol=1e-5)

[veni-vidi-vici-dormivi]

Okay, but can we disregard the covariates?

[mathause]

Okay, but can we disregard the covariates?

I don't know¹. It is a bit academic but helped me to understand the whole thing better. We could use fft to compute the first_guess potentially for all frequencies at once (as the fft is - well - fast, in contrast to my idea in #492). Although I am not sure if covariates for the first frequency influences the others.

For me this is a bit cleaner in any case & not very costly as it's not released in this form.

Actually, I just wondered if we should reorganize coeffs even a bit more?

state	$a_0$	$a_1$	$b_0$	$b_1$
current	$\sin \cdot \mathrm{T}$	$\sin$	$\cos \cdot \mathrm{T}$	$\cos$
this PR	$\cos \cdot \mathrm{T}$	$\cos$	$\sin \cdot \mathrm{T}$	$\sin$
even better?	$\cos$	$\sin$	$\cos \cdot \mathrm{T}$	$\sin \cdot \mathrm{T}$

---

¹We could find out if the covariates can be disregarded, also computing the BIC or deviance with and without the covariates.

[mathause]

Ok let's keep the order of coeffs as suggested here (cosT, cos, sinT, sin). Is the PR ok for you, @veni-vidi-vici-dormivi?

[veni-vidi-vici-dormivi]

Yes definitely thank you, but maybe we should open an issue with all the information of this PR and #516?

[mathause]

I thought if there was a better order of the coeffs for point 4 in #519 - but I did not come to any conclusion. So I suggest we merge as is and leave that for the future (if ever...).