We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Small issue in the text of kernel in the heteroskedastic example
*Are there any mistakes in the docs?*
Hi,

In
https://gpflow.readthedocs.io/en/develop/notebooks/advanced/heteroskedastic.html
the explanation of the kernel has an issue when rendering. In my Chrome, it looks like
"withbothkernelsbeingmodeledasseparateandindependent:math:‘SquaredExponential‘kernels."

It does look ok here:
https://github.com/GPflow/docs/blob/develop/doc/source/notebooks/advanced/heteroskedastic.ipynb

Kind regards, 

Ivan


--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `doc/source/notebooks/advanced/heteroskedastic.pct.py`
Content:
```
1 # ---
2 # jupyter:
3 #   jupytext:
4 #     formats: ipynb,.pct.py:percent
5 #     text_representation:
6 #       extension: .py
7 #       format_name: percent
8 #       format_version: '1.3'
9 #       jupytext_version: 1.6.0
10 #   kernelspec:
11 #     display_name: Python 3
12 #     language: python
13 #     name: python3
14 # ---
15 
16 # %% [markdown]
17 # # Heteroskedastic Likelihood and Multi-Latent GP
18 
19 # %% [markdown]
20 # ## Standard (Homoskedastic) Regression
21 # In standard GP regression, the GP latent function is used to learn the location parameter of a likelihood distribution (usually a Gaussian) as a function of the input $x$, whereas the scale parameter is considered constant. This is a homoskedastic model, which is unable to capture variations of the noise distribution with the input $x$.
22 #
23 #
24 # ## Heteroskedastic Regression
25 # This notebooks shows how to construct a model which uses multiple (here two) GP latent functions to learn both the location and the scale of the Gaussian likelihood distribution. It does so by connecting a **Multi-Output Kernel**, which generates multiple GP latent functions, to a **Heteroskedastic Likelihood**, which maps the latent GPs into a single likelihood.
26 #
27 # The generative model is described as:
28 #
29 # $$ f_1(x) \sim \mathcal{GP}(0, k_1(\cdot, \cdot)) $$
30 # $$ f_2(x) \sim \mathcal{GP}(0, k_2(\cdot, \cdot)) $$
31 # $$ \text{loc}(x) = f_1(x) $$
32 # $$ \text{scale}(x) = \text{transform}(f_2(x)) $$
33 # $$ y_i|f_1, f_2, x_i \sim \mathcal{N}(\text{loc}(x_i),\;\text{scale}(x_i)^2)$$
34 #
35 # The function $\text{transform}$ is used to map from the unconstrained GP $f_2$ to **positive-only values**, which is required as it represents the $\text{scale}$ of a Gaussian likelihood. In this notebook, the $\exp$ function will be used as the $\text{transform}$. Other positive transforms such as the $\text{softplus}$ function can also be used.
36 
37 # %%
38 import matplotlib.pyplot as plt
39 import numpy as np
40 import tensorflow as tf
41 import tensorflow_probability as tfp
42 import gpflow as gpf
43 
44 
45 # %% [markdown]
46 # ## Data Generation
47 # We generate heteroskedastic data by substituting the random latent functions $f_1$ and $f_2$ of the generative model by deterministic $\sin$ and $\cos$ functions. The input $X$ is built with $N=1001$ uniformly spaced values in the interval $[0, 4\pi]$. The outputs $Y$ are still sampled from a Gaussian likelihood.
48 #
49 # $$ x_i \in [0, 4\pi], \quad i = 1,\dots,N $$
50 # $$ f_1(x) = \sin(x) $$
51 # $$ f_2(x) = \cos(x) $$
52 # $$ \text{loc}(x) = f_1(x) $$
53 # $$ \text{scale}(x) = \exp(f_2(x)) $$
54 # $$ y_i|x_i \sim \mathcal{N}(\text{loc}(x_i),\;\text{scale}(x_i)^2)$$
55 
56 # %%
57 N = 1001
58 
59 np.random.seed(0)
60 tf.random.set_seed(0)
61 
62 # Build inputs X
63 X = np.linspace(0, 4 * np.pi, N)[:, None]  # X must be of shape [N, 1]
64 
65 # Deterministic functions in place of latent ones
66 f1 = np.sin
67 f2 = np.cos
68 
69 # Use transform = exp to ensure positive-only scale values
70 transform = np.exp
71 
72 # Compute loc and scale as functions of input X
73 loc = f1(X)
74 scale = transform(f2(X))
75 
76 # Sample outputs Y from Gaussian Likelihood
77 Y = np.random.normal(loc, scale)
78 
79 # %% [markdown]
80 # ### Plot Data
81 # Note how the distribution density (shaded area) and the outputs $Y$ both change depending on the input $X$.
82 
83 # %%
84 def plot_distribution(X, Y, loc, scale):
85     plt.figure(figsize=(15, 5))
86     x = X.squeeze()
87     for k in (1, 2):
88         lb = (loc - k * scale).squeeze()
89         ub = (loc + k * scale).squeeze()
90         plt.fill_between(x, lb, ub, color="silver", alpha=1 - 0.05 * k ** 3)
91     plt.plot(x, lb, color="silver")
92     plt.plot(x, ub, color="silver")
93     plt.plot(X, loc, color="black")
94     plt.scatter(X, Y, color="gray", alpha=0.8)
95     plt.show()
96     plt.close()
97 
98 
99 plot_distribution(X, Y, loc, scale)
100 
101 
102 # %% [markdown]
103 # ## Build Model
104 
105 # %% [markdown]
106 # ### Likelihood
107 # This implements the following part of the generative model:
108 # $$ \text{loc}(x) = f_1(x) $$
109 # $$ \text{scale}(x) = \text{transform}(f_2(x)) $$
110 # $$ y_i|f_1, f_2, x_i \sim \mathcal{N}(\text{loc}(x_i),\;\text{scale}(x_i)^2)$$
111 
112 # %%
113 likelihood = gpf.likelihoods.HeteroskedasticTFPConditional(
114     distribution_class=tfp.distributions.Normal,  # Gaussian Likelihood
115     scale_transform=tfp.bijectors.Exp(),  # Exponential Transform
116 )
117 
118 print(f"Likelihood's expected latent_dim: {likelihood.latent_dim}")
119 
120 # %% [markdown]
121 # ### Kernel
122 # This implements the following part of the generative model:
123 # $$ f_1(x) \sim \mathcal{GP}(0, k_1(\cdot, \cdot)) $$
124 # $$ f_2(x) \sim \mathcal{GP}(0, k_2(\cdot, \cdot)) $$
125 # with both kernels being modeled as separate and independent $\text{SquaredExponential}$ kernels.
126 
127 # %%
128 kernel = gpf.kernels.SeparateIndependent(
129     [
130         gpf.kernels.SquaredExponential(),  # This is k1, the kernel of f1
131         gpf.kernels.SquaredExponential(),  # this is k2, the kernel of f2
132     ]
133 )
134 # The number of kernels contained in gpf.kernels.SeparateIndependent must be the same as likelihood.latent_dim
135 
136 # %% [markdown]
137 # ### Inducing Points
138 # Since we will use the **SVGP** model to perform inference, we need to implement the inducing variables $U_1$ and $U_2$, both with size $M=20$, which are used to approximate $f_1$ and $f_2$ respectively, and initialize the inducing points positions $Z_1$ and $Z_2$. This gives a total of $2M=40$ inducing variables and inducing points.
139 #
140 # The inducing variables and their corresponding inputs will be Separate and Independent, but both $Z_1$ and $Z_2$ will be initialized as $Z$, which are placed as $M=20$ equally spaced points in $[\min(X), \max(X)]$.
141 #
142 
143 # %%
144 M = 20  # Number of inducing variables for each f_i
145 
146 # Initial inducing points position Z
147 Z = np.linspace(X.min(), X.max(), M)[:, None]  # Z must be of shape [M, 1]
148 
149 inducing_variable = gpf.inducing_variables.SeparateIndependentInducingVariables(
150     [
151         gpf.inducing_variables.InducingPoints(Z),  # This is U1 = f1(Z1)
152         gpf.inducing_variables.InducingPoints(Z),  # This is U2 = f2(Z2)
153     ]
154 )
155 
156 # %% [markdown]
157 # ### SVGP Model
158 # Build the **SVGP** model by composing the **Kernel**, the **Likelihood** and the **Inducing Variables**.
159 #
160 # Note that the model needs to be instructed about the number of latent GPs by passing `num_latent_gps=likelihood.latent_dim`.
161 
162 # %%
163 model = gpf.models.SVGP(
164     kernel=kernel,
165     likelihood=likelihood,
166     inducing_variable=inducing_variable,
167     num_latent_gps=likelihood.latent_dim,
168 )
169 
170 model
171 
172 # %% [markdown]
173 # ## Model Optimization
174 #
175 # ### Build Optimizers (NatGrad + Adam)
176 
177 # %%
178 data = (X, Y)
179 loss_fn = model.training_loss_closure(data)
180 
181 gpf.utilities.set_trainable(model.q_mu, False)
182 gpf.utilities.set_trainable(model.q_sqrt, False)
183 
184 variational_vars = [(model.q_mu, model.q_sqrt)]
185 natgrad_opt = gpf.optimizers.NaturalGradient(gamma=0.1)
186 
187 adam_vars = model.trainable_variables
188 adam_opt = tf.optimizers.Adam(0.01)
189 
190 
191 @tf.function
192 def optimisation_step():
193     natgrad_opt.minimize(loss_fn, variational_vars)
194     adam_opt.minimize(loss_fn, adam_vars)
195 
196 
197 # %% [markdown]
198 # ### Run Optimization Loop
199 
200 # %%
201 epochs = 100
202 log_freq = 20
203 
204 for epoch in range(1, epochs + 1):
205     optimisation_step()
206 
207     # For every 'log_freq' epochs, print the epoch and plot the predictions against the data
208     if epoch % log_freq == 0 and epoch > 0:
209         print(f"Epoch {epoch} - Loss: {loss_fn().numpy() : .4f}")
210         Ymean, Yvar = model.predict_y(X)
211         Ymean = Ymean.numpy().squeeze()
212         Ystd = tf.sqrt(Yvar).numpy().squeeze()
213         plot_distribution(X, Y, Ymean, Ystd)
214 
215 model
216 
217 # %% [markdown]
218 # ## Further reading
219 #
220 # See [Chained Gaussian Processes](http://proceedings.mlr.press/v51/saul16.html) by Saul et al. (AISTATS 2016).
221
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/doc/source/notebooks/advanced/heteroskedastic.pct.py b/doc/source/notebooks/advanced/heteroskedastic.pct.py
--- a/doc/source/notebooks/advanced/heteroskedastic.pct.py
+++ b/doc/source/notebooks/advanced/heteroskedastic.pct.py
@@ -122,6 +122,7 @@
 # This implements the following part of the generative model:
 # $$ f_1(x) \sim \mathcal{GP}(0, k_1(\cdot, \cdot)) $$
 # $$ f_2(x) \sim \mathcal{GP}(0, k_2(\cdot, \cdot)) $$
+#
 # with both kernels being modeled as separate and independent $\text{SquaredExponential}$ kernels.
 
 # %%

{"golden_diff": "diff --git a/doc/source/notebooks/advanced/heteroskedastic.pct.py b/doc/source/notebooks/advanced/heteroskedastic.pct.py\n--- a/doc/source/notebooks/advanced/heteroskedastic.pct.py\n+++ b/doc/source/notebooks/advanced/heteroskedastic.pct.py\n@@ -122,6 +122,7 @@\n # This implements the following part of the generative model:\n # $$ f_1(x) \\sim \\mathcal{GP}(0, k_1(\\cdot, \\cdot)) $$\n # $$ f_2(x) \\sim \\mathcal{GP}(0, k_2(\\cdot, \\cdot)) $$\n+#\n # with both kernels being modeled as separate and independent $\\text{SquaredExponential}$ kernels.\n \n # %%\n", "issue": "Small issue in the text of kernel in the heteroskedastic example\n*Are there any mistakes in the docs?*\r\nHi,\r\n\r\nIn\r\nhttps://gpflow.readthedocs.io/en/develop/notebooks/advanced/heteroskedastic.html\r\nthe explanation of the kernel has an issue when rendering. In my Chrome, it looks like\r\n\"withbothkernelsbeingmodeledasseparateandindependent:math:\u2018SquaredExponential\u2018kernels.\"\r\n\r\nIt does look ok here:\r\nhttps://github.com/GPflow/docs/blob/develop/doc/source/notebooks/advanced/heteroskedastic.ipynb\r\n\r\nKind regards, \r\n\r\nIvan\r\n\n", "before_files": [{"content": "# ---\n# jupyter:\n#   jupytext:\n#     formats: ipynb,.pct.py:percent\n#     text_representation:\n#       extension: .py\n#       format_name: percent\n#       format_version: '1.3'\n#       jupytext_version: 1.6.0\n#   kernelspec:\n#     display_name: Python 3\n#     language: python\n#     name: python3\n# ---\n\n# %% [markdown]\n# # Heteroskedastic Likelihood and Multi-Latent GP\n\n# %% [markdown]\n# ## Standard (Homoskedastic) Regression\n# In standard GP regression, the GP latent function is used to learn the location parameter of a likelihood distribution (usually a Gaussian) as a function of the input $x$, whereas the scale parameter is considered constant. This is a homoskedastic model, which is unable to capture variations of the noise distribution with the input $x$.\n#\n#\n# ## Heteroskedastic Regression\n# This notebooks shows how to construct a model which uses multiple (here two) GP latent functions to learn both the location and the scale of the Gaussian likelihood distribution. It does so by connecting a **Multi-Output Kernel**, which generates multiple GP latent functions, to a **Heteroskedastic Likelihood**, which maps the latent GPs into a single likelihood.\n#\n# The generative model is described as:\n#\n# $$ f_1(x) \\sim \\mathcal{GP}(0, k_1(\\cdot, \\cdot)) $$\n# $$ f_2(x) \\sim \\mathcal{GP}(0, k_2(\\cdot, \\cdot)) $$\n# $$ \\text{loc}(x) = f_1(x) $$\n# $$ \\text{scale}(x) = \\text{transform}(f_2(x)) $$\n# $$ y_i|f_1, f_2, x_i \\sim \\mathcal{N}(\\text{loc}(x_i),\\;\\text{scale}(x_i)^2)$$\n#\n# The function $\\text{transform}$ is used to map from the unconstrained GP $f_2$ to **positive-only values**, which is required as it represents the $\\text{scale}$ of a Gaussian likelihood. In this notebook, the $\\exp$ function will be used as the $\\text{transform}$. Other positive transforms such as the $\\text{softplus}$ function can also be used.\n\n# %%\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport tensorflow as tf\nimport tensorflow_probability as tfp\nimport gpflow as gpf\n\n\n# %% [markdown]\n# ## Data Generation\n# We generate heteroskedastic data by substituting the random latent functions $f_1$ and $f_2$ of the generative model by deterministic $\\sin$ and $\\cos$ functions. The input $X$ is built with $N=1001$ uniformly spaced values in the interval $[0, 4\\pi]$. The outputs $Y$ are still sampled from a Gaussian likelihood.\n#\n# $$ x_i \\in [0, 4\\pi], \\quad i = 1,\\dots,N $$\n# $$ f_1(x) = \\sin(x) $$\n# $$ f_2(x) = \\cos(x) $$\n# $$ \\text{loc}(x) = f_1(x) $$\n# $$ \\text{scale}(x) = \\exp(f_2(x)) $$\n# $$ y_i|x_i \\sim \\mathcal{N}(\\text{loc}(x_i),\\;\\text{scale}(x_i)^2)$$\n\n# %%\nN = 1001\n\nnp.random.seed(0)\ntf.random.set_seed(0)\n\n# Build inputs X\nX = np.linspace(0, 4 * np.pi, N)[:, None]  # X must be of shape [N, 1]\n\n# Deterministic functions in place of latent ones\nf1 = np.sin\nf2 = np.cos\n\n# Use transform = exp to ensure positive-only scale values\ntransform = np.exp\n\n# Compute loc and scale as functions of input X\nloc = f1(X)\nscale = transform(f2(X))\n\n# Sample outputs Y from Gaussian Likelihood\nY = np.random.normal(loc, scale)\n\n# %% [markdown]\n# ### Plot Data\n# Note how the distribution density (shaded area) and the outputs $Y$ both change depending on the input $X$.\n\n# %%\ndef plot_distribution(X, Y, loc, scale):\n    plt.figure(figsize=(15, 5))\n    x = X.squeeze()\n    for k in (1, 2):\n        lb = (loc - k * scale).squeeze()\n        ub = (loc + k * scale).squeeze()\n        plt.fill_between(x, lb, ub, color=\"silver\", alpha=1 - 0.05 * k ** 3)\n    plt.plot(x, lb, color=\"silver\")\n    plt.plot(x, ub, color=\"silver\")\n    plt.plot(X, loc, color=\"black\")\n    plt.scatter(X, Y, color=\"gray\", alpha=0.8)\n    plt.show()\n    plt.close()\n\n\nplot_distribution(X, Y, loc, scale)\n\n\n# %% [markdown]\n# ## Build Model\n\n# %% [markdown]\n# ### Likelihood\n# This implements the following part of the generative model:\n# $$ \\text{loc}(x) = f_1(x) $$\n# $$ \\text{scale}(x) = \\text{transform}(f_2(x)) $$\n# $$ y_i|f_1, f_2, x_i \\sim \\mathcal{N}(\\text{loc}(x_i),\\;\\text{scale}(x_i)^2)$$\n\n# %%\nlikelihood = gpf.likelihoods.HeteroskedasticTFPConditional(\n    distribution_class=tfp.distributions.Normal,  # Gaussian Likelihood\n    scale_transform=tfp.bijectors.Exp(),  # Exponential Transform\n)\n\nprint(f\"Likelihood's expected latent_dim: {likelihood.latent_dim}\")\n\n# %% [markdown]\n# ### Kernel\n# This implements the following part of the generative model:\n# $$ f_1(x) \\sim \\mathcal{GP}(0, k_1(\\cdot, \\cdot)) $$\n# $$ f_2(x) \\sim \\mathcal{GP}(0, k_2(\\cdot, \\cdot)) $$\n# with both kernels being modeled as separate and independent $\\text{SquaredExponential}$ kernels.\n\n# %%\nkernel = gpf.kernels.SeparateIndependent(\n    [\n        gpf.kernels.SquaredExponential(),  # This is k1, the kernel of f1\n        gpf.kernels.SquaredExponential(),  # this is k2, the kernel of f2\n    ]\n)\n# The number of kernels contained in gpf.kernels.SeparateIndependent must be the same as likelihood.latent_dim\n\n# %% [markdown]\n# ### Inducing Points\n# Since we will use the **SVGP** model to perform inference, we need to implement the inducing variables $U_1$ and $U_2$, both with size $M=20$, which are used to approximate $f_1$ and $f_2$ respectively, and initialize the inducing points positions $Z_1$ and $Z_2$. This gives a total of $2M=40$ inducing variables and inducing points.\n#\n# The inducing variables and their corresponding inputs will be Separate and Independent, but both $Z_1$ and $Z_2$ will be initialized as $Z$, which are placed as $M=20$ equally spaced points in $[\\min(X), \\max(X)]$.\n#\n\n# %%\nM = 20  # Number of inducing variables for each f_i\n\n# Initial inducing points position Z\nZ = np.linspace(X.min(), X.max(), M)[:, None]  # Z must be of shape [M, 1]\n\ninducing_variable = gpf.inducing_variables.SeparateIndependentInducingVariables(\n    [\n        gpf.inducing_variables.InducingPoints(Z),  # This is U1 = f1(Z1)\n        gpf.inducing_variables.InducingPoints(Z),  # This is U2 = f2(Z2)\n    ]\n)\n\n# %% [markdown]\n# ### SVGP Model\n# Build the **SVGP** model by composing the **Kernel**, the **Likelihood** and the **Inducing Variables**.\n#\n# Note that the model needs to be instructed about the number of latent GPs by passing `num_latent_gps=likelihood.latent_dim`.\n\n# %%\nmodel = gpf.models.SVGP(\n    kernel=kernel,\n    likelihood=likelihood,\n    inducing_variable=inducing_variable,\n    num_latent_gps=likelihood.latent_dim,\n)\n\nmodel\n\n# %% [markdown]\n# ## Model Optimization\n#\n# ### Build Optimizers (NatGrad + Adam)\n\n# %%\ndata = (X, Y)\nloss_fn = model.training_loss_closure(data)\n\ngpf.utilities.set_trainable(model.q_mu, False)\ngpf.utilities.set_trainable(model.q_sqrt, False)\n\nvariational_vars = [(model.q_mu, model.q_sqrt)]\nnatgrad_opt = gpf.optimizers.NaturalGradient(gamma=0.1)\n\nadam_vars = model.trainable_variables\nadam_opt = tf.optimizers.Adam(0.01)\n\n\n@tf.function\ndef optimisation_step():\n    natgrad_opt.minimize(loss_fn, variational_vars)\n    adam_opt.minimize(loss_fn, adam_vars)\n\n\n# %% [markdown]\n# ### Run Optimization Loop\n\n# %%\nepochs = 100\nlog_freq = 20\n\nfor epoch in range(1, epochs + 1):\n    optimisation_step()\n\n    # For every 'log_freq' epochs, print the epoch and plot the predictions against the data\n    if epoch % log_freq == 0 and epoch > 0:\n        print(f\"Epoch {epoch} - Loss: {loss_fn().numpy() : .4f}\")\n        Ymean, Yvar = model.predict_y(X)\n        Ymean = Ymean.numpy().squeeze()\n        Ystd = tf.sqrt(Yvar).numpy().squeeze()\n        plot_distribution(X, Y, Ymean, Ystd)\n\nmodel\n\n# %% [markdown]\n# ## Further reading\n#\n# See [Chained Gaussian Processes](http://proceedings.mlr.press/v51/saul16.html) by Saul et al. (AISTATS 2016).\n", "path": "doc/source/notebooks/advanced/heteroskedastic.pct.py"}], "after_files": [{"content": "# ---\n# jupyter:\n#   jupytext:\n#     formats: ipynb,.pct.py:percent\n#     text_representation:\n#       extension: .py\n#       format_name: percent\n#       format_version: '1.3'\n#       jupytext_version: 1.6.0\n#   kernelspec:\n#     display_name: Python 3\n#     language: python\n#     name: python3\n# ---\n\n# %% [markdown]\n# # Heteroskedastic Likelihood and Multi-Latent GP\n\n# %% [markdown]\n# ## Standard (Homoskedastic) Regression\n# In standard GP regression, the GP latent function is used to learn the location parameter of a likelihood distribution (usually a Gaussian) as a function of the input $x$, whereas the scale parameter is considered constant. This is a homoskedastic model, which is unable to capture variations of the noise distribution with the input $x$.\n#\n#\n# ## Heteroskedastic Regression\n# This notebooks shows how to construct a model which uses multiple (here two) GP latent functions to learn both the location and the scale of the Gaussian likelihood distribution. It does so by connecting a **Multi-Output Kernel**, which generates multiple GP latent functions, to a **Heteroskedastic Likelihood**, which maps the latent GPs into a single likelihood.\n#\n# The generative model is described as:\n#\n# $$ f_1(x) \\sim \\mathcal{GP}(0, k_1(\\cdot, \\cdot)) $$\n# $$ f_2(x) \\sim \\mathcal{GP}(0, k_2(\\cdot, \\cdot)) $$\n# $$ \\text{loc}(x) = f_1(x) $$\n# $$ \\text{scale}(x) = \\text{transform}(f_2(x)) $$\n# $$ y_i|f_1, f_2, x_i \\sim \\mathcal{N}(\\text{loc}(x_i),\\;\\text{scale}(x_i)^2)$$\n#\n# The function $\\text{transform}$ is used to map from the unconstrained GP $f_2$ to **positive-only values**, which is required as it represents the $\\text{scale}$ of a Gaussian likelihood. In this notebook, the $\\exp$ function will be used as the $\\text{transform}$. Other positive transforms such as the $\\text{softplus}$ function can also be used.\n\n# %%\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport tensorflow as tf\nimport tensorflow_probability as tfp\nimport gpflow as gpf\n\n\n# %% [markdown]\n# ## Data Generation\n# We generate heteroskedastic data by substituting the random latent functions $f_1$ and $f_2$ of the generative model by deterministic $\\sin$ and $\\cos$ functions. The input $X$ is built with $N=1001$ uniformly spaced values in the interval $[0, 4\\pi]$. The outputs $Y$ are still sampled from a Gaussian likelihood.\n#\n# $$ x_i \\in [0, 4\\pi], \\quad i = 1,\\dots,N $$\n# $$ f_1(x) = \\sin(x) $$\n# $$ f_2(x) = \\cos(x) $$\n# $$ \\text{loc}(x) = f_1(x) $$\n# $$ \\text{scale}(x) = \\exp(f_2(x)) $$\n# $$ y_i|x_i \\sim \\mathcal{N}(\\text{loc}(x_i),\\;\\text{scale}(x_i)^2)$$\n\n# %%\nN = 1001\n\nnp.random.seed(0)\ntf.random.set_seed(0)\n\n# Build inputs X\nX = np.linspace(0, 4 * np.pi, N)[:, None]  # X must be of shape [N, 1]\n\n# Deterministic functions in place of latent ones\nf1 = np.sin\nf2 = np.cos\n\n# Use transform = exp to ensure positive-only scale values\ntransform = np.exp\n\n# Compute loc and scale as functions of input X\nloc = f1(X)\nscale = transform(f2(X))\n\n# Sample outputs Y from Gaussian Likelihood\nY = np.random.normal(loc, scale)\n\n# %% [markdown]\n# ### Plot Data\n# Note how the distribution density (shaded area) and the outputs $Y$ both change depending on the input $X$.\n\n# %%\ndef plot_distribution(X, Y, loc, scale):\n    plt.figure(figsize=(15, 5))\n    x = X.squeeze()\n    for k in (1, 2):\n        lb = (loc - k * scale).squeeze()\n        ub = (loc + k * scale).squeeze()\n        plt.fill_between(x, lb, ub, color=\"silver\", alpha=1 - 0.05 * k ** 3)\n    plt.plot(x, lb, color=\"silver\")\n    plt.plot(x, ub, color=\"silver\")\n    plt.plot(X, loc, color=\"black\")\n    plt.scatter(X, Y, color=\"gray\", alpha=0.8)\n    plt.show()\n    plt.close()\n\n\nplot_distribution(X, Y, loc, scale)\n\n\n# %% [markdown]\n# ## Build Model\n\n# %% [markdown]\n# ### Likelihood\n# This implements the following part of the generative model:\n# $$ \\text{loc}(x) = f_1(x) $$\n# $$ \\text{scale}(x) = \\text{transform}(f_2(x)) $$\n# $$ y_i|f_1, f_2, x_i \\sim \\mathcal{N}(\\text{loc}(x_i),\\;\\text{scale}(x_i)^2)$$\n\n# %%\nlikelihood = gpf.likelihoods.HeteroskedasticTFPConditional(\n    distribution_class=tfp.distributions.Normal,  # Gaussian Likelihood\n    scale_transform=tfp.bijectors.Exp(),  # Exponential Transform\n)\n\nprint(f\"Likelihood's expected latent_dim: {likelihood.latent_dim}\")\n\n# %% [markdown]\n# ### Kernel\n# This implements the following part of the generative model:\n# $$ f_1(x) \\sim \\mathcal{GP}(0, k_1(\\cdot, \\cdot)) $$\n# $$ f_2(x) \\sim \\mathcal{GP}(0, k_2(\\cdot, \\cdot)) $$\n#\n# with both kernels being modeled as separate and independent $\\text{SquaredExponential}$ kernels.\n\n# %%\nkernel = gpf.kernels.SeparateIndependent(\n    [\n        gpf.kernels.SquaredExponential(),  # This is k1, the kernel of f1\n        gpf.kernels.SquaredExponential(),  # this is k2, the kernel of f2\n    ]\n)\n# The number of kernels contained in gpf.kernels.SeparateIndependent must be the same as likelihood.latent_dim\n\n# %% [markdown]\n# ### Inducing Points\n# Since we will use the **SVGP** model to perform inference, we need to implement the inducing variables $U_1$ and $U_2$, both with size $M=20$, which are used to approximate $f_1$ and $f_2$ respectively, and initialize the inducing points positions $Z_1$ and $Z_2$. This gives a total of $2M=40$ inducing variables and inducing points.\n#\n# The inducing variables and their corresponding inputs will be Separate and Independent, but both $Z_1$ and $Z_2$ will be initialized as $Z$, which are placed as $M=20$ equally spaced points in $[\\min(X), \\max(X)]$.\n#\n\n# %%\nM = 20  # Number of inducing variables for each f_i\n\n# Initial inducing points position Z\nZ = np.linspace(X.min(), X.max(), M)[:, None]  # Z must be of shape [M, 1]\n\ninducing_variable = gpf.inducing_variables.SeparateIndependentInducingVariables(\n    [\n        gpf.inducing_variables.InducingPoints(Z),  # This is U1 = f1(Z1)\n        gpf.inducing_variables.InducingPoints(Z),  # This is U2 = f2(Z2)\n    ]\n)\n\n# %% [markdown]\n# ### SVGP Model\n# Build the **SVGP** model by composing the **Kernel**, the **Likelihood** and the **Inducing Variables**.\n#\n# Note that the model needs to be instructed about the number of latent GPs by passing `num_latent_gps=likelihood.latent_dim`.\n\n# %%\nmodel = gpf.models.SVGP(\n    kernel=kernel,\n    likelihood=likelihood,\n    inducing_variable=inducing_variable,\n    num_latent_gps=likelihood.latent_dim,\n)\n\nmodel\n\n# %% [markdown]\n# ## Model Optimization\n#\n# ### Build Optimizers (NatGrad + Adam)\n\n# %%\ndata = (X, Y)\nloss_fn = model.training_loss_closure(data)\n\ngpf.utilities.set_trainable(model.q_mu, False)\ngpf.utilities.set_trainable(model.q_sqrt, False)\n\nvariational_vars = [(model.q_mu, model.q_sqrt)]\nnatgrad_opt = gpf.optimizers.NaturalGradient(gamma=0.1)\n\nadam_vars = model.trainable_variables\nadam_opt = tf.optimizers.Adam(0.01)\n\n\n@tf.function\ndef optimisation_step():\n    natgrad_opt.minimize(loss_fn, variational_vars)\n    adam_opt.minimize(loss_fn, adam_vars)\n\n\n# %% [markdown]\n# ### Run Optimization Loop\n\n# %%\nepochs = 100\nlog_freq = 20\n\nfor epoch in range(1, epochs + 1):\n    optimisation_step()\n\n    # For every 'log_freq' epochs, print the epoch and plot the predictions against the data\n    if epoch % log_freq == 0 and epoch > 0:\n        print(f\"Epoch {epoch} - Loss: {loss_fn().numpy() : .4f}\")\n        Ymean, Yvar = model.predict_y(X)\n        Ymean = Ymean.numpy().squeeze()\n        Ystd = tf.sqrt(Yvar).numpy().squeeze()\n        plot_distribution(X, Y, Ymean, Ystd)\n\nmodel\n\n# %% [markdown]\n# ## Further reading\n#\n# See [Chained Gaussian Processes](http://proceedings.mlr.press/v51/saul16.html) by Saul et al. (AISTATS 2016).\n", "path": "doc/source/notebooks/advanced/heteroskedastic.pct.py"}]}