二氧化碳高斯过程在冒纳罗亚#

此高斯过程 (GP) 示例展示了如何

设计协方差函数的组合
使用加性 GP，其各个组件可用于预测
执行最大后验 (MAP) 估计

自 1950 年代后期以来，冒纳罗亚天文台一直在定期测量大气中的 CO\(_2\)。在 1950 年代后期，查尔斯·基林发明了一种测量大气 CO\(_2\) 浓度的精确方法。从那时起，冒纳罗亚天文台几乎持续记录 CO\(_2\) 的测量值。查看最近几小时的测量结果此处。

在 1950 年代后期，人们对化石燃料燃烧如何影响气候知之甚少。最初几年的数据收集表明，CO\(_2\) 水平随着夏季和冬季而上升和下降，跟踪北半球植被的生长和衰败。随着多年过去，持续上升的趋势越来越受到关注。凭借 70 多年收集的数据，基林曲线已成为最重要的气候指标之一。

关于这些测量背后的历史及其对当今气候学的影响以及其他有趣的阅读资料

http://scrippsco2.ucsd.edu/history_legacy/early_keeling_curve#
https://scripps.ucsd.edu/programs/keelingcurve/2016/05/23/why-has-a-drop-in-global-co2-emissions-not-caused-co2-levels-in-the-atmosphere-to-stabilize/#more-1412

让我们加载数据，整理一下，然后看看。原始数据集位于此处。此笔记本使用 Bokeh 包进行从交互性中获益的绘图。

准备数据#

注意

此笔记本使用的库不是 PyMC 依赖项，因此需要专门安装才能运行此笔记本。打开下面的下拉菜单以获取更多指导。

import numpy as np
import pandas as pd
import pymc3 as pm

from bokeh.io import output_notebook
from bokeh.models import BoxAnnotation, Label, Legend, Span
from bokeh.palettes import brewer
from bokeh.plotting import figure, show

output_notebook()

正在加载 BokehJS ...

%config InlineBackend.figure_format = 'retina'
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)

# get data
try:
    data_monthly = pd.read_csv("../data/monthly_in_situ_co2_mlo.csv", header=56)
except FileNotFoundError:
    data_monthly = pd.read_csv(pm.get_data("monthly_in_situ_co2_mlo.csv"), header=56)

# replace -99.99 with NaN
data_monthly.replace(to_replace=-99.99, value=np.nan, inplace=True)

# fix column names
cols = [
    "year",
    "month",
    "--",
    "--",
    "CO2",
    "seasonaly_adjusted",
    "fit",
    "seasonally_adjusted_fit",
    "CO2_filled",
    "seasonally_adjusted_filled",
]
data_monthly.columns = cols
cols.remove("--")
cols.remove("--")
data_monthly = data_monthly[cols]

# drop rows with nan
data_monthly.dropna(inplace=True)

# fix time index
data_monthly["day"] = 15
data_monthly.index = pd.to_datetime(data_monthly[["year", "month", "day"]])
cols.remove("year")
cols.remove("month")
data_monthly = data_monthly[cols]

data_monthly.head(5)

	CO2	季节性调整	拟合	季节性调整拟合	CO2_filled	季节性调整填充
1958-03-15	315.71	314.44	316.19	314.91	315.71	314.44
1958-04-15	317.45	315.16	317.30	314.99	317.45	315.16
1958-05-15	317.51	314.70	317.87	315.07	317.51	314.70
1958-07-15	315.86	315.20	315.86	315.22	315.86	315.20
1958-08-15	314.93	316.20	313.98	315.29	314.93	316.20

# function to convert datetimes to indexed numbers that are useful for later prediction
def dates_to_idx(timelist):
    reference_time = pd.to_datetime("1958-03-15")
    t = (timelist - reference_time) / pd.Timedelta(365, "D")
    return np.asarray(t)


t = dates_to_idx(data_monthly.index)

# normalize CO2 levels
y = data_monthly["CO2"].values
first_co2 = y[0]
std_co2 = np.std(y)
y_n = (y - first_co2) / std_co2

data_monthly = data_monthly.assign(t=t)
data_monthly = data_monthly.assign(y_n=y_n)

您可能对这些数据很熟悉，因为它在 Gaussian Processes for Machine Learning 一书中被用作示例，作者是 Rasmussen 和 Williams [2006]。他们使用的数据集版本从 1950 年代后期开始，但在 2003 年底结束。为了使我们的 PyMC3 示例与他们的示例具有一定的可比性，我们使用 2004 年之前的数据段作为“训练”集。我们将使用 2004 年至 2022 年的数据来测试我们的预测。

# split into training and test set
sep_idx = data_monthly.index.searchsorted(pd.to_datetime("2003-12-15"))
data_early = data_monthly.iloc[: sep_idx + 1, :]
data_later = data_monthly.iloc[sep_idx:, :]

# plot training and test data
p = figure(
    x_axis_type="datetime",
    title="Monthly CO2 Readings from Mauna Loa",
    plot_width=550,
    plot_height=350,
)
p.yaxis.axis_label = "CO2 [ppm]"
p.xaxis.axis_label = "Date"
predict_region = BoxAnnotation(
    left=pd.to_datetime("2003-12-15"), fill_alpha=0.1, fill_color="firebrick"
)
p.add_layout(predict_region)
ppm400 = Span(location=400, dimension="width", line_color="red", line_dash="dashed", line_width=2)
p.add_layout(ppm400)

p.line(data_monthly.index, data_monthly["CO2"], line_width=2, line_color="black", alpha=0.5)
p.circle(data_monthly.index, data_monthly["CO2"], line_color="black", alpha=0.1, size=2)

train_label = Label(
    x=100,
    y=165,
    x_units="screen",
    y_units="screen",
    text="Training Set",
    render_mode="css",
    border_line_alpha=0.0,
    background_fill_alpha=0.0,
)
test_label = Label(
    x=585,
    y=80,
    x_units="screen",
    y_units="screen",
    text="Test Set",
    render_mode="css",
    border_line_alpha=0.0,
    background_fill_alpha=0.0,
)

p.add_layout(train_label)
p.add_layout(test_label)
show(p)

Bokeh 图是交互式的，因此可以使用右侧边栏进行平移和缩放。季节性的上升和下降以及上升趋势都显而易见。这是一个链接，指向此曲线在不同时间尺度上的绘图，以及在历史冰芯数据背景下的绘图。

400 ppm 水平用虚线突出显示。除了拟合描述性模型外，我们的目标是预测首次超过 400 ppm 阈值的月份，即 2013 年 5 月。在上面的数据集中，2013 年 5 月的 CO\(_2\) 平均读数约为 399.98，足够接近我们的正确目标日期。

使用 GP 对基林曲线建模#

作为起点，我们使用 Rasmussen 和 Williams [2006] 中描述的 GP 模型。与教科书中对协方差函数超参数使用平坦先验，然后最大化边际似然性不同，我们在超参数上放置了一些信息性先验，并使用优化来找到 MAP 点。我们使用 gp.Marginal，因为假设存在高斯噪声。

R&W [Rasmussen 和 Williams，2006] 模型是信号的三个 GP 和噪声的一个 GP 的总和。

由指数二次核表示的长期平滑上升趋势。
一个周期性项，其周期性从精确周期性衰减。这由 Periodic 和 Matern52 协方差函数的乘积表示。
具有有理二次核的中小型不规则性。
噪声建模为 Matern32 和白噪声核之和。

CO\(_2\) 作为时间函数的先验是，

\[ f(t) \sim \mathcal{GP}_{\text{slow}}(0,\, k_1(t, t')) + \mathcal{GP}_{\text{med}}(0,\, k_2(t, t')) + \mathcal{GP}_{\text{per}}(0,\, k_3(t, t')) + \mathcal{GP}_{\text{noise}}(0,\, k_4(t, t')) \]

超参数先验#

我们对协方差函数的尺度超参数使用相当无信息的先验，对长度尺度使用信息性 Gamma 参数。用于长度尺度先验的 PDF 如下所示

x = np.linspace(0, 150, 5000)
priors = [
    ("ℓ_pdecay", pm.Gamma.dist(alpha=10, beta=0.075)),
    ("ℓ_psmooth", pm.Gamma.dist(alpha=4, beta=3)),
    ("period", pm.Normal.dist(mu=1.0, sigma=0.05)),
    ("ℓ_med", pm.Gamma.dist(alpha=2, beta=0.75)),
    ("α", pm.Gamma.dist(alpha=5, beta=2)),
    ("ℓ_trend", pm.Gamma.dist(alpha=4, beta=0.1)),
    ("ℓ_noise", pm.Gamma.dist(alpha=2, beta=4)),
]

colors = brewer["Paired"][7]

p = figure(
    title="Lengthscale and period priors",
    plot_width=550,
    plot_height=350,
    x_range=(-1, 8),
    y_range=(0, 2),
)
p.yaxis.axis_label = "Probability"
p.xaxis.axis_label = "Years"

for i, prior in enumerate(priors):
    p.line(
        x,
        np.exp(prior[1].logp(x).eval()),
        legend_label=prior[0],
        line_width=3,
        line_color=colors[i],
    )
show(p)

ℓ_pdecay：周期性衰减。此参数越小，周期性消失得越快。我怀疑 CO\(_2\) 的季节性会在短期内消失（希望如此），并且数据中也没有证据表明这一点。大部分先验质量都在 60 年到 >140 年之间。
ℓ_psmooth：周期性成分的平滑度。它控制周期性的“正弦”程度。数据显示季节性不是精确的正弦波，但与正弦波没有太大区别。我们使用模式为 1 的 Gamma 分布，并且方差不太大，大部分先验质量在 0.5 到 2 左右。
period：周期。我们对 \(p\) 施加了非常强的先验，周期以 1 为中心。R&W 将 \(p=1\) 固定，因为周期是年度周期。
ℓ_med：这是短期到中期长期变化的长度尺度。此先验的大部分质量都在 6 年以下。
α：这是形状参数。此先验以 3 为中心，因为我们预计会存在比指数二次函数可以解释的更多的变化。
ℓ_trend：长期趋势的长度尺度。它具有广泛的先验，质量在十年尺度上。大部分质量在 10 年到 60 年之间。
ℓ_noise：噪声协方差的长度尺度。此噪声应非常快速，在几个月到最多一两年内。

我们事先知道哪些 GP 成分应该具有更大的幅度，因此我们将此信息包含在尺度参数中。

x = np.linspace(0, 4, 5000)
priors = [
    ("η_per", pm.HalfCauchy.dist(beta=2)),
    ("η_med", pm.HalfCauchy.dist(beta=1.0)),
    (
        "η_trend",
        pm.HalfCauchy.dist(beta=3),
    ),  # will use beta=2, but beta=3 is visible on plot
    ("σ", pm.HalfNormal.dist(sigma=0.25)),
    ("η_noise", pm.HalfNormal.dist(sigma=0.5)),
]

colors = brewer["Paired"][5]

p = figure(title="Scale priors", plot_width=550, plot_height=350)
p.yaxis.axis_label = "Probability"
p.xaxis.axis_label = "Years"

for i, prior in enumerate(priors):
    p.line(
        x,
        np.exp(prior[1].logp(x).eval()),
        legend_label=prior[0],
        line_width=3,
        line_color=colors[i],
    )
show(p)

对于所有尺度先验，我们都使用将尺度缩小为零的分布。季节性成分和长期趋势在零附近具有最少的质量，因为它们是数据中最大的影响因素。

η_per：周期性或季节性成分的尺度。
η_med：短期到中期成分的尺度。
η_trend：长期趋势的尺度。
σ：白噪声的尺度。
η_noise：相关的短期噪声的尺度。

PyMC3 中的模型#

下面是实际模型。三个组件 GP 中的每一个都是单独构建的。由于我们正在进行 MAP，因此我们使用 Marginal GP，最后调用 .marginal_likelihood 方法来指定边际后验。

# pull out normalized data
t = data_early["t"].values[:, None]
y = data_early["y_n"].values

with pm.Model() as model:
    # yearly periodic component x long term trend
    η_per = pm.HalfCauchy("η_per", beta=2, testval=1.0)
    ℓ_pdecay = pm.Gamma("ℓ_pdecay", alpha=10, beta=0.075)
    period = pm.Normal("period", mu=1, sigma=0.05)
    ℓ_psmooth = pm.Gamma("ℓ_psmooth ", alpha=4, beta=3)
    cov_seasonal = (
        η_per**2 * pm.gp.cov.Periodic(1, period, ℓ_psmooth) * pm.gp.cov.Matern52(1, ℓ_pdecay)
    )
    gp_seasonal = pm.gp.Marginal(cov_func=cov_seasonal)

    # small/medium term irregularities
    η_med = pm.HalfCauchy("η_med", beta=0.5, testval=0.1)
    ℓ_med = pm.Gamma("ℓ_med", alpha=2, beta=0.75)
    α = pm.Gamma("α", alpha=5, beta=2)
    cov_medium = η_med**2 * pm.gp.cov.RatQuad(1, ℓ_med, α)
    gp_medium = pm.gp.Marginal(cov_func=cov_medium)

    # long term trend
    η_trend = pm.HalfCauchy("η_trend", beta=2, testval=2.0)
    ℓ_trend = pm.Gamma("ℓ_trend", alpha=4, beta=0.1)
    cov_trend = η_trend**2 * pm.gp.cov.ExpQuad(1, ℓ_trend)
    gp_trend = pm.gp.Marginal(cov_func=cov_trend)

    # noise model
    η_noise = pm.HalfNormal("η_noise", sigma=0.5, testval=0.05)
    ℓ_noise = pm.Gamma("ℓ_noise", alpha=2, beta=4)
    σ = pm.HalfNormal("σ", sigma=0.25, testval=0.05)
    cov_noise = η_noise**2 * pm.gp.cov.Matern32(1, ℓ_noise) + pm.gp.cov.WhiteNoise(σ)

    # The Gaussian process is a sum of these three components
    gp = gp_seasonal + gp_medium + gp_trend

    # Since the normal noise model and the GP are conjugates, we use `Marginal` with the `.marginal_likelihood` method
    y_ = gp.marginal_likelihood("y", X=t, y=y, noise=cov_noise)

    # this line calls an optimizer to find the MAP
    mp = pm.find_MAP(include_transformed=True)

100.00% [300/300 00:52<00:00 logp = 1,734.4, ||grad|| = 14.264]

# display the results, dont show transformed parameter values
sorted([name + ":" + str(mp[name]) for name in mp.keys() if not name.endswith("_")])

['period:1.00038447904791',
 'α:1.1815982319316372',
 'η_med:0.021391276520692636',
 'η_noise:0.006819703084865676',
 'η_per:0.08803544888083512',
 'η_trend:1.712437847900024',
 'σ:0.0059989944162813155',
 'ℓ_med:1.507055362902511',
 'ℓ_noise:0.16455363060019101',
 'ℓ_pdecay:125.58803198611457',
 'ℓ_psmooth :0.7229240021307682',
 'ℓ_trend:52.58970494079544']

乍一看，结果看起来是合理的。决定季节性变化速度的长度尺度约为 126 年。这意味着根据数据，我们预计如此强的周期性在几个世纪后才会消失。趋势长度尺度也很长，约为 50 年。

检查每个加性 GP 组件的拟合#

下面的代码查看总 GP 以及每个组件的单独拟合。总拟合及其 \(2\sigma\) 不确定性以红色显示。

# predict at a 15 day granularity
dates = pd.date_range(start="3/15/1958", end="12/15/2003", freq="15D")
tnew = dates_to_idx(dates)[:, None]

print("Predicting with gp ...")
mu, var = gp.predict(tnew, point=mp, diag=True)
mean_pred = mu * std_co2 + first_co2
var_pred = var * std_co2**2

# make dataframe to store fit results
fit = pd.DataFrame(
    {"t": tnew.flatten(), "mu_total": mean_pred, "sd_total": np.sqrt(var_pred)},
    index=dates,
)

print("Predicting with gp_trend ...")
mu, var = gp_trend.predict(
    tnew, point=mp, given={"gp": gp, "X": t, "y": y, "noise": cov_noise}, diag=True
)
fit = fit.assign(mu_trend=mu * std_co2 + first_co2, sd_trend=np.sqrt(var * std_co2**2))

print("Predicting with gp_medium ...")
mu, var = gp_medium.predict(
    tnew, point=mp, given={"gp": gp, "X": t, "y": y, "noise": cov_noise}, diag=True
)
fit = fit.assign(mu_medium=mu * std_co2 + first_co2, sd_medium=np.sqrt(var * std_co2**2))

print("Predicting with gp_seasonal ...")
mu, var = gp_seasonal.predict(
    tnew, point=mp, given={"gp": gp, "X": t, "y": y, "noise": cov_noise}, diag=True
)
fit = fit.assign(mu_seasonal=mu * std_co2 + first_co2, sd_seasonal=np.sqrt(var * std_co2**2))
print("Done")

Predicting with gp ...
Predicting with gp_trend ...
Predicting with gp_medium ...
Predicting with gp_seasonal ...
Done

## plot the components
p = figure(
    title="Decomposition of the Mauna Loa Data",
    x_axis_type="datetime",
    plot_width=550,
    plot_height=350,
)
p.yaxis.axis_label = "CO2 [ppm]"
p.xaxis.axis_label = "Date"

# plot mean and 2σ region of total prediction
upper = fit.mu_total + 2 * fit.sd_total
lower = fit.mu_total - 2 * fit.sd_total
band_x = np.append(fit.index.values, fit.index.values[::-1])
band_y = np.append(lower, upper[::-1])

# total fit
p.line(
    fit.index,
    fit.mu_total,
    line_width=1,
    line_color="firebrick",
    legend_label="Total fit",
)
p.patch(band_x, band_y, color="firebrick", alpha=0.6, line_color="white")

# trend
p.line(
    fit.index,
    fit.mu_trend,
    line_width=1,
    line_color="blue",
    legend_label="Long term trend",
)

# medium
p.line(
    fit.index,
    fit.mu_medium,
    line_width=1,
    line_color="green",
    legend_label="Medium range variation",
)

# seasonal
p.line(
    fit.index,
    fit.mu_seasonal,
    line_width=1,
    line_color="orange",
    legend_label="Seasonal process",
)

# true value
p.circle(data_early.index, data_early["CO2"], color="black", legend_label="Observed data")
p.legend.location = "top_left"
show(p)

拟合与观察到的数据非常吻合。趋势、季节性和短期/中期效应也清晰地分离出来。如果放大以使季节性过程填充绘图窗口（从 x 等于 1955 到 2004，从 y 等于 310 到 320），则它似乎随着时间的推移而变宽。让我们绘制每十年第一年的情况

# plot several years

p = figure(title="Several years of the seasonal component", plot_width=550, plot_height=350)
p.yaxis.axis_label = "Δ CO2 [ppm]"
p.xaxis.axis_label = "Month"

colors = brewer["Paired"][5]
years = ["1960", "1970", "1980", "1990", "2000"]

for i, year in enumerate(years):
    dates = pd.date_range(start="1/1/" + year, end="12/31/" + year, freq="10D")
    tnew = dates_to_idx(dates)[:, None]

    print("Predicting year", year)
    mu, var = gp_seasonal.predict(
        tnew, point=mp, diag=True, given={"gp": gp, "X": t, "y": y, "noise": cov_noise}
    )
    mu_pred = mu * std_co2

    # plot mean
    x = np.asarray((dates - dates[0]) / pd.Timedelta(30, "D")) + 1
    p.line(x, mu_pred, line_width=1, line_color=colors[i], legend_label=year)

p.legend.location = "bottom_left"
show(p)

Predicting year 1960
Predicting year 1970
Predicting year 1980
Predicting year 1990
Predicting year 2000

此图清楚地表明随着时间的推移存在变宽。因此，似乎随着大气中 CO\(_2\) 的增加，由于北半球植被的生长和衰败而引起的吸收/释放周期变得稍微更加明显。

二氧化碳水平将在哪一天突破 400 ppm？#

我们的预测看起来如何？显然，观察到的数据呈上升趋势，季节性影响非常明显。我们的 GP 模型是否能够很好地捕捉到这一点以进行合理的推断？我们的“训练”集一直到 2003 年底，因此我们将预测从 2004 年 1 月到 2022 年底。

尽管此事件除了是一个漂亮的整数外没有任何特别的意义，但我们的次要目标是看看我们能否很好地预测首次突破 400 ppm 标记的日期。此事件首次发生在 2013 年 5 月，并且有一些关于其他重要里程碑的新闻文章。

dates = pd.date_range(start="11/15/2003", end="12/15/2022", freq="10D")
tnew = dates_to_idx(dates)[:, None]

print("Sampling gp predictions ...")
mu_pred, cov_pred = gp.predict(tnew, point=mp)

# draw samples, and rescale
n_samples = 2000
samples = pm.MvNormal.dist(mu=mu_pred, cov=cov_pred, shape=(n_samples, len(tnew))).random()
samples = samples * std_co2 + first_co2

Sampling gp predictions ...

# make plot
p = figure(x_axis_type="datetime", plot_width=700, plot_height=300)
p.yaxis.axis_label = "CO2 [ppm]"
p.xaxis.axis_label = "Date"

# plot mean and 2σ region of total prediction
# scale mean and var
mu_pred_sc = mu_pred * std_co2 + first_co2
sd_pred_sc = np.sqrt(np.diag(cov_pred) * std_co2**2)

upper = mu_pred_sc + 2 * sd_pred_sc
lower = mu_pred_sc - 2 * sd_pred_sc
band_x = np.append(dates, dates[::-1])
band_y = np.append(lower, upper[::-1])

p.line(dates, mu_pred_sc, line_width=2, line_color="firebrick", legend_label="Total fit")
p.patch(band_x, band_y, color="firebrick", alpha=0.6, line_color="white")

# some predictions
idx = np.random.randint(0, samples.shape[0], 10)
p.multi_line(
    [dates] * len(idx),
    [samples[i, :] for i in idx],
    color="firebrick",
    alpha=0.5,
    line_width=0.5,
)

# true value
p.circle(data_later.index, data_later["CO2"], color="black", legend_label="Observed data")

ppm400 = Span(
    location=400,
    dimension="width",
    line_color="black",
    line_dash="dashed",
    line_width=1,
)
p.add_layout(ppm400)
p.legend.location = "bottom_right"
show(p)

平均预测和 \(2\sigma\) 不确定性以红色显示。边际后验的几个样本也显示在那里。看起来我们的模型对二氧化碳的释放量有些过于乐观。 \(2\sigma\) 不确定性首次超过 400 ppm 阈值是在 2015 年 5 月，晚了两年。

发生这种情况的一个原因是我们的 GP 先验具有零均值。这意味着我们编码了先验信息，表明当我们远离观察到的数据时，该函数应趋于零。这个假设可能是不合理的。CO\(_2\) 趋势也可能比线性增长更快——这是准确预测的重要知识。另一种可能性是 MAP 估计。在不查看完整后验的情况下，我们估计中的不确定性被低估了。严重程度未知。

具有零均值 GP 先验会导致预测偏差很大。解决此问题的一些可能性是使用常数均值函数，其值可能会分配给历史的或工业革命前的 CO\(_2\) 平均值。但这可能不是未来 CO\(_2\) 水平的最佳指标。

此外，仅使用历史 CO\(_2\) 数据可能不是最佳预测指标。除了使用 GP 拟合来查看决定 CO\(_2\) 水平的潜在行为外，我们还可以纳入其他信息，例如化石燃料燃烧释放的 CO\(_2\) 量。

接下来，我们将了解如何使用 PyMC3 的 GP 功能来改进模型，查看完整后验，并纳入关于 CO\(_2\) 水平驱动因素的其他数据来源。

作者#

由 Bill Engels 于 2017 年 9 月编写 (pymc#2444)
由 Chris Fonnesbeck 于 2020 年 12 月更新
由 Danh Phan 于 2022 年 5 月重新执行 (pymc-examples#316)

参考文献#

[1] (1,2,3)

Carl Edward Rasmussen 和 Christopher K. I. Williams. Gaussian Processes for Machine Learning. The MIT Press, 2006. ISBN 026218253X. URL: https://gaussianprocess.org/gpml/.

%load_ext watermark
%watermark -n -u -v -iv -w -p bokeh

Last updated: Fri May 06 2022

Python implementation: CPython
Python version       : 3.9.5
IPython version      : 7.25.0

bokeh: 2.4.2

numpy : 1.21.0
pandas: 1.4.0
pymc3 : 3.11.4

Watermark: 2.3.0

许可声明#

此示例库中的所有笔记本均根据 MIT 许可证提供，该许可证允许修改和再分发以用于任何用途，前提是版权和许可声明得到保留。

引用 PyMC 示例#

要引用此笔记本，请使用 Zenodo 为 pymc-examples 存储库提供的 DOI。

重要提示

许多笔记本都改编自其他来源：博客、书籍……在这种情况下，您也应该引用原始来源。

另请记住引用您的代码使用的相关库。

这是一个 bibtex 中的引用模板

@incollection{citekey,
  author    = "<notebook authors, see above>",
  title     = "<notebook title>",
  editor    = "PyMC Team",
  booktitle = "PyMC examples",
  doi       = "10.5281/zenodo.5654871"
}

渲染后可能看起来像

分类

标签