import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib.dates import DateFormatter
import matplotlib.dates as mdates
import matplotlib.ticker as ticker
import seaborn as sns
sns.set(style="ticks", font_scale=1.5) # white graphs, with large and legible letters
# %matplotlib widget# dirty trick to have dates in the middle of the 24-hour period
# make minor ticks in the middle, put the labels there!
# from https://matplotlib.org/stable/gallery/ticks/centered_ticklabels.html
def centered_dates(ax):
date_form = DateFormatter("%d %b") # %d 3-letter-Month
# major ticks at midnight, every day
ax.xaxis.set_major_locator(mdates.DayLocator(interval=1))
ax.xaxis.set_major_formatter(date_form)
# minor ticks at noon, every day
ax.xaxis.set_minor_locator(mdates.HourLocator(byhour=[12]))
# erase major tick labels
ax.xaxis.set_major_formatter(ticker.NullFormatter())
# set minor tick labels as define above
ax.xaxis.set_minor_formatter(date_form)
# completely erase minor ticks, center tick labels
for tick in ax.xaxis.get_minor_ticks():
tick.tick1line.set_markersize(0)
tick.tick2line.set_markersize(0)
tick.label1.set_horizontalalignment('center')
# creating the dictionary with the desired settings
plot_settings = {
'ylim': [5, 17.5],
'xlim': [df.index[0], df.index[-1]],
'ylabel': 'Temperature (°C)',
'title': 'Yatir Forest, 2022',
'yticks': [5, 10, 15]
}Let’s see on a graph the average temperature, with an envelope of 1 standard deviation around it:
fig, ax = plt.subplots(figsize=(8,5))
plot_std = ax.fill_between(df.index,
df['mean'] + df['std'],
df['mean'] - df['std'],
color="xkcd:pink", alpha=0.5)
plot_data, = ax.plot(df['temp'], color='black')
plot_mean, =ax.plot(df['mean'], color='xkcd:hot pink')
ax.legend([plot_data, plot_mean, plot_std],
['data', '3-hour average', r"$\pm1$ std envelope"],
frameon=False)
# applying the settings to the ax object
ax.set(**plot_settings)
centered_dates(ax)
# fig.savefig("YF-temperature_2022_jan.png", dpi=300)
We can calculate anything we want inside the sliding window. One good example is the Confidence Interval of the Mean, given by:
CI(\alpha) = Z(\alpha) \cdot SE.
This is called “רווח בר-סֶמֶך” in hebrew.
Z(\alpha) is the Z-score corresponding to the chosen confidence level \alpha. The most commonly used confidence level is 95%, which corresponds to a Z-score of 1.96. What does this mean? For a normal distribution, 95% of the points are within \pm 1.96 standard deviations from the mean. But be careful: the confidence interval of the mean is not about where the points are, it is about how well we know the mean. It is 1.96 times the standard error, and it gets narrower when we have more points.
Source: Dhaval Raval’s Medium article
You can find the Z-score using the following python code:
from scipy.stats import norm
confidence_level = 0.95
# 5% outside
out = 1 - confidence_level
# 0.975 of points to the left of right boundary
p = 1 - out/2
# inverse of cdf: 0.975 of the points will be smaller than what distance (in sigma units)?
z_score = norm.ppf(p)
print(f"z-score = {z_score}")z-score = 1.959963984540054If you are still not convinced why we need 0.975 instead of 0.95, read this excellent response on stackoverflow.
SE is the standard error:
SE = \frac{\sigma}{ \sqrt{N} }.
We can write a function to calculate the confidence interval of the mean, and use it with the sliding window:
def std_error_of_the_mean(window):
return window.std() / np.sqrt(window.count())
def confidence_interval(window):
return z_score * std_error_of_the_mean(window)
df['std_error'] = (
df['temp'].rolling('3h',
center=True)
.apply(std_error_of_the_mean)
)
df['confidence_int'] = (
df['temp'].rolling('3h',
center=True)
.apply(confidence_interval)
)fig, ax = plt.subplots(figsize=(8,5))
plot_std = ax.fill_between(df.index,
df['mean'] + df['confidence_int'],
df['mean'] - df['confidence_int'],
color="xkcd:pink", alpha=0.5)
plot_data, = ax.plot(df['temp'], color='black', alpha=0.3)
plot_mean, =ax.plot(df['mean'], color='xkcd:hot pink')
ax.legend([plot_data, plot_mean, plot_std],
['data', '3-hour running average', r"95% confidence interval"],
frameon=False)
# applying the settings to the ax object
ax.set(**plot_settings)
centered_dates(ax)
# fig.savefig("YF-temperature_2022_jan.png", dpi=300)
When the time series has a regular sampling frequency, all positions of the running window will have the same number of data points in them. Because the Confidence Interval is proportional to the Standard Error, and the SE is proportional to the Standard Deviation (\sqrt{N} is constant), then the envelope created by the CI is identical to the envelope created by the standard deviation, up to a multiplying constant. Nice.
Here the constant is 1.96/\sqrt{N}, and with N=18 points in a 3-hour window it is about 0.46: the shaded region of the second graph is less than half as wide as the one of the first graph, and most of the data lies outside of it. This is not a mistake. Remember that the confidence interval says how well we know the mean, not where the points are.
One more warning: the standard error assumes that the points in the window are independent of each other, but for temperature each point is very similar to the previous one. So the real uncertainty of the mean is larger than what we calculated. We will come back to this when we talk about autocorrelation.
fig, ax = plt.subplots(figsize=(8,5))
plot_ci, = ax.plot(df['confidence_int'], color='tab:red')
plot_std, = ax.plot(df['std'], color="black")
ax.legend([plot_ci, plot_std],
['confidence interval', 'standard deviation'],
frameon=False)
# applying the settings to the ax object
# ax.set(**plot_settings)
ax.set(xlim=[df.index[0], df.index[-1]])
centered_dates(ax)