チュートリアル 2: 連続隠れ状態を持つベイズ推論と意思決定
第3週、第2日目: ベイズ意思決定
Neuromatch Academyによる
コンテンツ制作者: Eric DeWitt, Xaq Pitkow, Saeed Salehi, Ella Batty
コンテンツレビュアー: Hyosub Kim, Zahra Arjmandi, Anoop Kulkarni
制作編集: Ella Batty, Gagana B, Spiros Chavlis
チュートリアルの目的
チュートリアルの推定所要時間: 1時間35分
このノートブックでは、連続分布に対するガウス分布とベイズの定理を紹介し、事前情報と新しい測定値のシンプルで強力な組み合わせをモデル化できるようにします。このノートブックでは、二値状態のチュートリアルで探求した同じアイデアを扱いますが、新しい問題「アストロキャットを見つけて導く」ことを紹介します。この問題は、今後の日々でより複雑な形で再び登場します。
このノートブックで学ぶこと:
- ガウス分布とその優れた性質について学ぶ
- 二値隠れ状態のチュートリアルのアイデアを連続分布に拡張する方法を探る
- 異なる事前分布がより複雑な事後分布を生み出す様子を探る
- 推論でよく使われる損失関数や複雑な効用関数を探る
# @title Tutorial slides
# @markdown These are the slides for the videos in this tutorial
from IPython.display import IFrame
link_id = "uz9rc"
print(f"If you want to download the slides: https://osf.io/download/{link_id}/")
IFrame(src=f"https://mfr.ca-1.osf.io/render?url=https://osf.io/{link_id}/?direct%26mode=render%26action=download%26mode=render", width=854, height=480)
セットアップ
# @title Install and import feedback gadget
from vibecheck import DatatopsContentReviewContainer
def content_review(notebook_section: str):
return DatatopsContentReviewContainer(
"", # No text prompt
notebook_section,
{
"url": "https://pmyvdlilci.execute-api.us-east-1.amazonaws.com/klab",
"name": "neuromatch_cn",
"user_key": "y1x3mpx5",
},
).render()
feedback_prefix = "W3D1_T2"
# Imports
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from scipy.stats import gamma as gamma_distribution
from matplotlib.transforms import Affine2D
# @title Figure Settings
import logging
logging.getLogger('matplotlib.font_manager').disabled = True
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import ipywidgets as widgets
from IPython.display import clear_output
from ipywidgets import FloatSlider, Dropdown, interactive_output
from ipywidgets import interact, fixed, HBox, Layout, VBox, interactive, Label
plt.style.use("https://raw.githubusercontent.com/NeuromatchAcademy/course-content/main/nma.mplstyle")
import warnings
warnings.filterwarnings("ignore")
# @title Plotting Functions
def plot_mixture_prior(x, gaussian1, gaussian2, combined):
""" Plots a prior made of a mixture of gaussians
Args:
x (numpy array of floats): points at which the likelihood has been evaluated
gaussian1 (numpy array of floats): normalized probabilities for Gaussian 1 evaluated at each `x`
gaussian2 (numpy array of floats): normalized probabilities for Gaussian 2 evaluated at each `x`
posterior (numpy array of floats): normalized probabilities for the posterior evaluated at each `x`
Returns:
Nothing
"""
fig, ax = plt.subplots()
ax.plot(x, gaussian1, '--b', LineWidth=2, label='Gaussian 1')
ax.plot(x, gaussian2, '-.b', LineWidth=2, label='Gaussian 2')
ax.plot(x, combined, '-r', LineWidth=2, label='Gaussian Mixture')
ax.legend()
ax.set_ylabel('Probability')
ax.set_xlabel('Orientation (Degrees)')
def plot_gaussian(mu, sigma):
x = np.linspace(-7, 7, 1000, endpoint=True)
y = gaussian(x, mu, sigma)
plt.figure(figsize=(6, 4))
plt.plot(x, y, c='blue')
plt.fill_between(x, y, color='b', alpha=0.2)
plt.ylabel('$\mathcal{N}(x, \mu, \sigma^2)$')
plt.xlabel('x')
plt.yticks([])
plt.show()
def plot_losses(mu, sigma):
x = np.linspace(-2, 2, 400, endpoint=True)
y = gaussian(x, mu, sigma)
error = x - mu
mse_loss = (error)**2
abs_loss = np.abs(error)
zero_one_loss = (np.abs(error) >= 0.02).astype(np.float)
fig, (ax_gaus, ax_error) = plt.subplots(2, 1, figsize=(6, 8))
ax_gaus.plot(x, y, color='blue', label='true distribution')
ax_gaus.fill_between(x, y, color='blue', alpha=0.2)
ax_gaus.set_ylabel('$\\mathcal{N}(x, \\mu, \\sigma^2)$')
ax_gaus.set_xlabel('x')
ax_gaus.set_yticks([])
ax_gaus.legend(loc='upper right')
ax_error.plot(x, mse_loss, color='c', label='Mean Squared Error', linewidth=3)
ax_error.plot(x, abs_loss, color='m', label='Absolute Error', linewidth=3)
ax_error.plot(x, zero_one_loss, color='y', label='Zero-One Loss', linewidth=3)
ax_error.legend(loc='upper right')
ax_error.set_xlabel('$\\hat{\\mu}$')
ax_error.set_ylabel('Error')
plt.show()
def gaussian_mixture(mu1, mu2, sigma1, sigma2, factor):
assert 0.0 < factor < 1.0
x = np.linspace(-7.0, 7.0, 1000, endpoint=True)
y_1 = gaussian(x, mu1, sigma1)
y_2 = gaussian(x, mu2, sigma2)
mixture = y_1 * factor + y_2 * (1.0 - factor)
plt.figure(figsize=(8, 6))
plt.plot(x, y_1, c='deepskyblue', label='p(x)', linewidth=3.0)
plt.fill_between(x, y_1, color='deepskyblue', alpha=0.2)
plt.plot(x, y_2, c='aquamarine', label='q(x)', linewidth=3.0)
plt.fill_between(x, y_2, color='aquamarine', alpha=0.2)
plt.plot(x, mixture, c='b', label='$\pi \cdot p(x) + (1-\pi) \cdot q(x)$', linewidth=3.0)
plt.fill_between(x, mixture, color='b', alpha=0.2)
plt.yticks([])
plt.legend(loc="upper left")
plt.xlabel('x')
plt.show()
def plot_utility_mixture_dist(mu1, sigma1, mu2, sigma2, mu_g, sigma_g,
mu_loc, mu_dist, plot_utility_row=True):
x = np.linspace(-7, 7, 1000, endpoint=True)
prior = gaussian(x, mu1, sigma1)
likelihood = gaussian(x, mu2, sigma2)
mu_post, sigma_post = product_guassian(mu1, mu2, sigma1, sigma2)
posterior = gaussian(x, mu_post, sigma_post)
gain = gaussian(x, mu_g, sigma_g)/2
sigma_mix, factor = 1.0, 0.5
mu_mix1, mu_mix2 = mu_loc - mu_dist/2, mu_loc + mu_dist/2
gaus_mix1, gaus_mix2 = gaussian(x, mu_mix1, sigma_mix), gaussian(x, mu_mix2, sigma_mix)
loss = factor * gaus_mix1 + (1 - factor) * gaus_mix2
utility = np.multiply(posterior, gain) - np.multiply(posterior, loss)
if plot_utility_row:
plot_bayes_utility_rows(x, prior, likelihood, posterior, gain, loss, utility)
else:
plot_bayes_row(x, prior, likelihood, posterior)
return None
def plot_mvn2d(mu1, mu2, sigma1, sigma2, corr):
x, y = np.mgrid[-2:2:.02, -2:2:.02]
cov12 = corr * sigma1 * sigma2
z = mvn2d(x, y, mu1, mu2, sigma1, sigma2, cov12)
plt.figure(figsize=(6, 6))
plt.contourf(x, y, z, cmap='Reds')
plt.axis("off")
plt.show()
def plot_marginal(sigma1, sigma2, c_x, c_y, corr):
mu1, mu2 = 0.0, 0.0
cov12 = corr * sigma1 * sigma2
xx, yy = np.mgrid[-2:2:.02, -2:2:.02]
x, y = xx[:, 0], yy[0]
p_x = gaussian(x, mu1, sigma1)
p_y = gaussian(y, mu2, sigma2)
zz = mvn2d(xx, yy, mu1, mu2, sigma1, sigma2, cov12)
mu_x_y = mu1+cov12*(c_y-mu2)/sigma2**2
mu_y_x = mu2+cov12*(c_x-mu1)/sigma1**2
sigma_x_y = np.sqrt(sigma1**2 - cov12**2/sigma2**2)
sigma_y_x = np.sqrt(sigma2**2 - cov12**2/sigma1**2)
p_x_y = gaussian(x, mu_x_y, sigma_x_y)
p_y_x = gaussian(x, mu_y_x, sigma_y_x)
p_c_y = gaussian(mu_x_y-sigma_x_y, mu_x_y, sigma_x_y)
p_c_x = gaussian(mu_y_x-sigma_y_x, mu_y_x, sigma_y_x)
# definitions for the axes
left, width = 0.1, 0.65
bottom, height = 0.1, 0.65
spacing = 0.01
rect_z = [left, bottom, width, height]
rect_x = [left, bottom + height + spacing, width, 0.2]
rect_y = [left + width + spacing, bottom, 0.2, height]
# start with a square Figure
fig = plt.figure(figsize=(8, 8))
ax_z = fig.add_axes(rect_z)
ax_x = fig.add_axes(rect_x, sharex=ax_z)
ax_y = fig.add_axes(rect_y, sharey=ax_z)
ax_z.set_axis_off()
ax_x.set_axis_off()
ax_y.set_axis_off()
ax_x.set_xlim(np.min(x), np.max(x))
ax_y.set_ylim(np.min(y), np.max(y))
ax_z.contourf(xx, yy, zz, cmap='Greys')
ax_z.hlines(c_y, mu_x_y-sigma_x_y, mu_x_y+sigma_x_y, color='c', zorder=9, linewidth=3)
ax_z.vlines(c_x, mu_y_x-sigma_y_x, mu_y_x+sigma_y_x, color='m', zorder=9, linewidth=3)
ax_x.plot(x, p_x, label='$p(x)$', c = 'b', linewidth=3)
ax_x.plot(x, p_x_y, label='$p(x|y = C_y)$', c = 'c', linestyle='dashed', linewidth=3)
ax_x.hlines(p_c_y, mu_x_y-sigma_x_y, mu_x_y+sigma_x_y, color='c', linestyle='dashed', linewidth=3)
ax_y.plot(p_y, y, label='$p(y)$', c = 'r', linewidth=3)
ax_y.plot(p_y_x, y, label='$p(y|x = C_x)$', c = 'm', linestyle='dashed', linewidth=3)
ax_y.vlines(p_c_x, mu_y_x-sigma_y_x, mu_y_x+sigma_y_x, color='m', linestyle='dashed', linewidth=3)
ax_x.legend(loc="upper left", frameon=False)
ax_y.legend(loc="lower right", frameon=False)
plt.show()
def plot_bayes(mu1, mu2, sigma1, sigma2):
x = np.linspace(-7, 7, 1000, endpoint=True)
prior = gaussian(x, mu1, sigma1)
likelihood = gaussian(x, mu2, sigma2)
mu_post, sigma_post = product_guassian(mu1, mu2, sigma1, sigma2)
posterior = gaussian(x, mu_post, sigma_post)
plt.figure(figsize=(8, 6))
plt.plot(x, prior, c='b', label='prior')
plt.fill_between(x, prior, color='b', alpha=0.2)
plt.plot(x, likelihood, c='r', label='likelihood')
plt.fill_between(x, likelihood, color='r', alpha=0.2)
plt.plot(x, posterior, c='k', label='posterior')
plt.fill_between(x, posterior, color='k', alpha=0.2)
plt.yticks([])
plt.legend(loc="upper left")
plt.ylabel('$\mathcal{N}(x, \mu, \sigma^2)$')
plt.xlabel('x')
plt.show()
def plot_information(mu1, sigma1, mu2, sigma2):
x = np.linspace(-7, 7, 1000, endpoint=True)
mu3, sigma3 = product_guassian(mu1, mu2, sigma1, sigma2)
prior = gaussian(x, mu1, sigma1)
likelihood = gaussian(x, mu2, sigma2)
posterior = gaussian(x, mu3, sigma3)
plt.figure(figsize=(8, 6))
plt.plot(x, prior, c='b', label='Satellite')
plt.fill_between(x, prior, color='b', alpha=0.2)
plt.plot(x, likelihood, c='r', label='Space Mouse')
plt.fill_between(x, likelihood, color='r', alpha=0.2)
plt.plot(x, posterior, c='k', label='Center')
plt.fill_between(x, posterior, color='k', alpha=0.2)
plt.yticks([])
plt.legend(loc="upper left")
plt.ylabel('$\mathcal{N}(x, \mu, \sigma^2)$')
plt.xlabel('x')
plt.show()
def plot_information_global(mu3, sigma3, mu1, mu2):
x = np.linspace(-7, 7, 1000, endpoint=True)
sigma1, sigma2 = reverse_product(mu3, sigma3, mu1, mu2)
prior = gaussian(x, mu1, sigma1)
likelihood = gaussian(x, mu2, sigma2)
posterior = gaussian(x, mu3, sigma3)
plt.figure(figsize=(8, 6))
plt.plot(x, prior, c='b', label='Satellite')
plt.fill_between(x, prior, color='b', alpha=0.2)
plt.plot(x, likelihood, c='r', label='Space Mouse')
plt.fill_between(x, likelihood, color='r', alpha=0.2)
plt.plot(x, posterior, c='k', label='Center')
plt.fill_between(x, posterior, color='k', alpha=0.2)
plt.yticks([])
plt.legend(loc="upper left")
plt.ylabel('$\mathcal{N}(x, \mu, \sigma^2)$')
plt.xlabel('x')
plt.show()
def plot_loss_utility_gaussian(loss_f, mu, sigma, mu_true):
x = np.linspace(-7, 7, 1000, endpoint=True)
posterior = gaussian(x, mu, sigma)
y_label = "$p(x)$"
plot_loss_utility(x, posterior, loss_f, mu_true, y_label)
def plot_loss_utility_mixture(loss_f, mu1, mu2, sigma1, sigma2, factor, mu_true):
x = np.linspace(-7, 7, 1000, endpoint=True)
y_1 = gaussian(x, mu1, sigma1)
y_2 = gaussian(x, mu2, sigma2)
posterior = y_1 * factor + y_2 * (1.0 - factor)
y_label = "$\pi \cdot p(x) + (1-\pi) \cdot q(x)$"
plot_loss_utility(x, posterior, loss_f, mu_true, y_label)
def plot_loss_utility(x, posterior, loss_f, mu_true, y_label):
mean, median, mode = calc_mean_mode_median(x, posterior)
loss = calc_loss_func(loss_f, mu_true, x)
utility = calc_expected_loss(loss_f, posterior, x)
min_expected_loss = x[np.argmin(utility)]
plt.figure(figsize=(12, 8))
plt.subplot(2, 2, 1)
plt.title("Probability")
plt.plot(x, posterior, c='b')
plt.fill_between(x, posterior, color='b', alpha=0.2)
plt.yticks([])
plt.xlabel('x')
plt.ylabel(y_label)
plt.axvline(mean, ls='dashed', color='red', label='Mean')
plt.axvline(median, ls='dashdot', color='blue', label='Median')
plt.axvline(mode, ls='dotted', color='green', label='Mode')
plt.legend(loc="upper left")
plt.subplot(2, 2, 2)
plt.title(loss_f)
plt.plot(x, loss, c='c', label=loss_f)
plt.ylabel('loss')
plt.xlabel('x')
plt.subplot(2, 2, 3)
plt.title("Expected Loss")
plt.plot(x, utility, c='y', label='$\mathbb{E}[L]$')
plt.axvline(min_expected_loss, ls='dashed', color='red', label='$Min~ \mathbb{E}[Loss]$')
plt.legend(loc="lower right")
plt.xlabel('x')
plt.ylabel('$\mathbb{E}[L]$')
plt.show()
def plot_loss_utility_bayes(mu1, mu2, sigma1, sigma2, mu_true, loss_f):
x = np.linspace(-4, 4, 1000, endpoint=True)
prior = gaussian(x, mu1, sigma1)
likelihood = gaussian(x, mu2, sigma2)
mu_post, sigma_post = product_guassian(mu1, mu2, sigma1, sigma2)
posterior = gaussian(x, mu_post, sigma_post)
loss = calc_loss_func(loss_f, mu_true, x)
utility = - calc_expected_loss(loss_f, posterior, x)
plt.figure(figsize=(18, 5))
plt.subplot(1, 3, 1)
plt.title("Posterior distribution")
plt.plot(x, prior, c='b', label='prior')
plt.fill_between(x, prior, color='b', alpha=0.2)
plt.plot(x, likelihood, c='r', label='likelihood')
plt.fill_between(x, likelihood, color='r', alpha=0.2)
plt.plot(x, posterior, c='k', label='posterior')
plt.fill_between(x, posterior, color='k', alpha=0.2)
plt.yticks([])
plt.legend(loc="upper left")
plt.xlabel('x')
plt.subplot(1, 3, 2)
plt.title(loss_f)
plt.plot(x, loss, c='c')
plt.ylabel('loss')
plt.subplot(1, 3, 3)
plt.title("Expected utility")
plt.plot(x, utility, c='y', label='utility')
plt.legend(loc="upper left")
plt.show()
def plot_simple_utility_gaussian(mu, sigma, mu_g, mu_c, sigma_g, sigma_c):
x = np.linspace(-7, 7, 1000, endpoint=True)
posterior = gaussian(x, mu, sigma)
gain = gaussian(x, mu_g, sigma_g)
loss = gaussian(x, mu_c, sigma_c)
utility = np.multiply(posterior, gain) - np.multiply(posterior, loss)
plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1)
plt.title("Probability")
plt.plot(x, posterior, c='b', label='posterior')
plt.fill_between(x, posterior, color='b', alpha=0.2)
plt.yticks([])
plt.xlabel('x')
plt.subplot(1, 3, 2)
plt.title("utility function")
plt.plot(x, gain, c='m', label='gain')
plt.plot(x, -loss, c='c', label='loss')
plt.legend(loc="upper left")
plt.subplot(1, 3, 3)
plt.title("expected utility")
plt.plot(x, utility, c='y', label='utility')
plt.legend(loc="upper left")
plt.show()
def plot_utility_gaussian(mu1, mu2, sigma1, sigma2, mu_g, mu_c, sigma_g, sigma_c, plot_utility_row=True):
x = np.linspace(-7, 7, 1000, endpoint=True)
prior = gaussian(x, mu1, sigma1)
likelihood = gaussian(x, mu2, sigma2)
mu_post, sigma_post = product_guassian(mu1, mu2, sigma1, sigma2)
posterior = gaussian(x, mu_post, sigma_post)
if plot_utility_row:
gain = gaussian(x, mu_g, sigma_g)
loss = gaussian(x, mu_c, sigma_c)
utility = np.multiply(posterior, gain) - np.multiply(posterior, loss)
plot_bayes_utility_rows(x, prior, likelihood, posterior, gain, loss, utility)
else:
plot_bayes_row(x, prior, likelihood, posterior)
return None
def plot_utility_mixture(mu_m1, mu_m2, sigma_m1, sigma_m2, factor,
mu, sigma, mu_g, mu_c, sigma_g, sigma_c, plot_utility_row=True):
x = np.linspace(-7, 7, 1000, endpoint=True)
y_1 = gaussian(x, mu_m1, sigma_m1)
y_2 = gaussian(x, mu_m2, sigma_m2)
prior = y_1 * factor + y_2 * (1.0 - factor)
likelihood = gaussian(x, mu, sigma)
posterior = np.multiply(prior, likelihood)
posterior = posterior / (posterior.sum() * (x[1] - x[0]))
if plot_utility_row:
gain = gaussian(x, mu_g, sigma_g)
loss = gaussian(x, mu_c, sigma_c)
utility = np.multiply(posterior, gain) - np.multiply(posterior, loss)
plot_bayes_utility_rows(x, prior, likelihood, posterior, gain, loss, utility)
else:
plot_bayes_row(x, prior, likelihood, posterior)
return None
def plot_utility_uniform(mu, sigma, mu_g, mu_c, sigma_g, sigma_c, plot_utility_row=True):
x = np.linspace(-7, 7, 1000, endpoint=True)
prior = np.ones_like(x) / (x.max() - x.min())
likelihood = gaussian(x, mu, sigma)
posterior = likelihood
if plot_utility_row:
gain = gaussian(x, mu_g, sigma_g)
loss = gaussian(x, mu_c, sigma_c)
utility = np.multiply(posterior, gain) - np.multiply(posterior, loss)
plot_bayes_utility_rows(x, prior, likelihood, posterior, gain, loss, utility)
else:
plot_bayes_row(x, prior, likelihood, posterior)
return None
def plot_utility_gamma(alpha, beta, offset, mu, sigma, mu_g, mu_c, sigma_g, sigma_c, plot_utility_row=True):
x = np.linspace(-4, 10, 1000, endpoint=True)
prior = gamma_pdf(x-offset, alpha, beta)
likelihood = gaussian(x, mu, sigma)
posterior = np.multiply(prior, likelihood)
posterior = posterior / (posterior.sum() * (x[1] - x[0]))
if plot_utility_row:
gain = gaussian(x, mu_g, sigma_g)
loss = gaussian(x, mu_c, sigma_c)
utility = np.multiply(posterior, gain) - np.multiply(posterior, loss)
plot_bayes_utility_rows(x, prior, likelihood, posterior, gain, loss, utility)
else:
plot_bayes_row(x, prior, likelihood, posterior)
return None
def plot_bayes_row(x, prior, likelihood, posterior):
mean, median, mode = calc_mean_mode_median(x, posterior)
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.title("Prior and likelihood distribution")
plt.plot(x, prior, c='b', label='prior')
plt.fill_between(x, prior, color='b', alpha=0.2)
plt.plot(x, likelihood, c='r', label='likelihood')
plt.fill_between(x, likelihood, color='r', alpha=0.2)
plt.yticks([])
plt.legend(loc="upper left")
plt.xlabel('x')
plt.subplot(1, 2, 2)
plt.title("Posterior distribution")
plt.plot(x, posterior, c='k', label='posterior')
plt.fill_between(x, posterior, color='k', alpha=0.1)
plt.axvline(mean, ls='dashed', color='red', label='Mean')
plt.axvline(median, ls='dashdot', color='blue', label='Median')
plt.axvline(mode, ls='dotted', color='green', label='Mode')
plt.legend(loc="upper left")
plt.yticks([])
plt.xlabel('x')
plt.show()
def plot_bayes_utility_rows(x, prior, likelihood, posterior, gain, loss, utility):
mean, median, mode = calc_mean_mode_median(x, posterior)
max_utility = x[np.argmax(utility)]
plt.figure(figsize=(12, 8))
plt.subplot(2, 2, 1)
plt.title("Prior and likelihood distribution")
plt.plot(x, prior, c='b', label='prior')
plt.fill_between(x, prior, color='b', alpha=0.2)
plt.plot(x, likelihood, c='r', label='likelihood')
plt.fill_between(x, likelihood, color='r', alpha=0.2)
# plt.plot(x, posterior, c='k', label='posterior')
# plt.fill_between(x, posterior, color='k', alpha=0.2)
plt.yticks([])
plt.legend(loc="upper left")
# plt.ylabel('$f(x)$')
plt.xlabel('x')
plt.subplot(2, 2, 2)
plt.title("Posterior distribution")
plt.plot(x, posterior, c='k', label='posterior')
plt.fill_between(x, posterior, color='k', alpha=0.1)
plt.axvline(mean, ls='dashed', color='red', label='Mean')
plt.axvline(median, ls='dashdot', color='blue', label='Median')
plt.axvline(mode, ls='dotted', color='green', label='Mode')
plt.legend(loc="upper left")
plt.yticks([])
plt.xlabel('x')
plt.subplot(2, 2, 3)
plt.title("utility function")
plt.plot(x, gain, c='m', label='gain')
# plt.fill_between(x, gain, color='m', alpha=0.2)
plt.plot(x, -loss, c='c', label='loss')
# plt.fill_between(x, -loss, color='c', alpha=0.2)
plt.legend(loc="upper left")
plt.xlabel('x')
plt.subplot(2, 2, 4)
plt.title("expected utility")
plt.plot(x, utility, c='y', label='utility')
# plt.fill_between(x, utility, color='y', alpha=0.2)
plt.axvline(max_utility, ls='dashed', color='red', label='Max utility')
plt.legend(loc="upper left")
plt.xlabel('x')
plt.ylabel('utility')
plt.legend(loc="lower right")
plt.show()
def plot_bayes_loss_utility_gaussian(loss_f, mu_true, mu1, mu2, sigma1, sigma2):
x = np.linspace(-7, 7, 1000, endpoint=True)
prior = gaussian(x, mu1, sigma1)
likelihood = gaussian(x, mu2, sigma2)
mu_post, sigma_post = product_guassian(mu1, mu2, sigma1, sigma2)
posterior = gaussian(x, mu_post, sigma_post)
loss = calc_loss_func(loss_f, mu_true, x)
plot_bayes_loss_utility(x, prior, likelihood, posterior, loss, loss_f)
return None
def plot_bayes_loss_utility_uniform(loss_f, mu_true, mu, sigma):
x = np.linspace(-7, 7, 1000, endpoint=True)
prior = np.ones_like(x) / (x.max() - x.min())
likelihood = gaussian(x, mu, sigma)
posterior = likelihood
loss = calc_loss_func(loss_f, mu_true, x)
plot_bayes_loss_utility(x, prior, likelihood, posterior, loss, loss_f)
return None
def plot_bayes_loss_utility_gamma(loss_f, mu_true, alpha, beta, offset, mu, sigma):
x = np.linspace(-7, 7, 1000, endpoint=True)
prior = gamma_pdf(x-offset, alpha, beta)
likelihood = gaussian(x, mu, sigma)
posterior = np.multiply(prior, likelihood)
posterior = posterior / (posterior.sum() * (x[1] - x[0]))
loss = calc_loss_func(loss_f, mu_true, x)
plot_bayes_loss_utility(x, prior, likelihood, posterior, loss, loss_f)
return None
def plot_bayes_loss_utility_mixture(loss_f, mu_true, mu_m1, mu_m2, sigma_m1, sigma_m2, factor, mu, sigma):
x = np.linspace(-7, 7, 1000, endpoint=True)
y_1 = gaussian(x, mu_m1, sigma_m1)
y_2 = gaussian(x, mu_m2, sigma_m2)
prior = y_1 * factor + y_2 * (1.0 - factor)
likelihood = gaussian(x, mu, sigma)
posterior = np.multiply(prior, likelihood)
posterior = posterior / (posterior.sum() * (x[1] - x[0]))
loss = calc_loss_func(loss_f, mu_true, x)
plot_bayes_loss_utility(x, prior, likelihood, posterior, loss, loss_f)
return None
def plot_bayes_loss_utility(x, prior, likelihood, posterior, loss, loss_f):
mean, median, mode = calc_mean_mode_median(x, posterior)
expected_loss = calc_expected_loss(loss_f, posterior, x)
min_expected_loss = x[np.argmin(expected_loss)]
plt.figure(figsize=(12, 8))
plt.subplot(2, 2, 1)
plt.title("Prior and Likelihood")
plt.plot(x, prior, c='b', label='prior')
plt.fill_between(x, prior, color='b', alpha=0.2)
plt.plot(x, likelihood, c='r', label='likelihood')
plt.fill_between(x, likelihood, color='r', alpha=0.2)
plt.yticks([])
plt.legend(loc="upper left")
plt.xlabel('x')
plt.subplot(2, 2, 2)
plt.title("Posterior")
plt.plot(x, posterior, c='k', label='posterior')
plt.fill_between(x, posterior, color='k', alpha=0.1)
plt.axvline(mean, ls='dashed', color='red', label='Mean')
plt.axvline(median, ls='dashdot', color='blue', label='Median')
plt.axvline(mode, ls='dotted', color='green', label='Mode')
plt.legend(loc="upper left")
plt.yticks([])
plt.xlabel('x')
plt.subplot(2, 2, 3)
plt.title(loss_f)
plt.plot(x, loss, c='c', label=loss_f)
# plt.fill_between(x, loss, color='c', alpha=0.2)
plt.ylabel('loss')
plt.xlabel('x')
plt.subplot(2, 2, 4)
plt.title("expected loss")
plt.plot(x, expected_loss, c='y', label='$\mathbb{E}[L]$')
# plt.fill_between(x, expected_loss, color='y', alpha=0.2)
plt.axvline(min_expected_loss, ls='dashed', color='red', label='$Min~ \mathbb{E}[Loss]$')
plt.legend(loc="lower right")
plt.xlabel('x')
plt.ylabel('$\mathbb{E}[L]$')
plt.show()
global global_loss_plot_switcher
global_loss_plot_switcher = False
def loss_plot_switcher(what_to_plot):
global global_loss_plot_switcher
if global_loss_plot_switcher:
clear_output()
else:
global_loss_plot_switcher = True
if what_to_plot == "Gaussian":
loss_f_options = Dropdown(
options=["Mean Squared Error", "Absolute Error", "Zero-One Loss"],
value="Mean Squared Error", description="Loss: ")
mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description="µ_estimate", continuous_update=True)
sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_estimate", continuous_update=True)
mu_true_slider = FloatSlider(min=-3.0, max=3.0, step=0.01, value=0.0, description="µ_true", continuous_update=True)
widget_ui = HBox([VBox([loss_f_options, mu_true_slider]), VBox([mu_slider, sigma_slider])])
widget_out = interactive_output(plot_loss_utility_gaussian,
{'loss_f': loss_f_options,
'mu': mu_slider,
'sigma': sigma_slider,
'mu_true': mu_true_slider})
elif what_to_plot == "Mixture of Gaussians":
loss_f_options = Dropdown(
options=["Mean Squared Error", "Absolute Error", "Zero-One Loss"],
value="Mean Squared Error", description="Loss: ")
mu1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description="µ_est_p", continuous_update=True)
mu2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description="µ_est_q", continuous_update=True)
sigma1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_est_p", continuous_update=True)
sigma2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_est_q", continuous_update=True)
factor_slider = FloatSlider(min=0.0, max=1.0, step=0.01, value=0.5, description="π", continuous_update=True)
mu_true_slider = FloatSlider(min=-3.0, max=3.0, step=0.01, value=0.0, description="µ_true", continuous_update=True)
widget_ui = HBox([VBox([loss_f_options, factor_slider, mu_true_slider]),
VBox([mu1_slider, sigma1_slider]),
VBox([mu2_slider, sigma2_slider])])
widget_out = interactive_output(plot_loss_utility_mixture,
{'mu1': mu1_slider,
'mu2': mu2_slider,
'sigma1': sigma1_slider,
'sigma2': sigma2_slider,
'factor': factor_slider,
'mu_true': mu_true_slider,
'loss_f': loss_f_options})
display(widget_ui, widget_out)
global global_plot_prior_switcher
global_plot_prior_switcher = False
def plot_prior_switcher(what_to_plot):
global global_plot_prior_switcher
if global_plot_prior_switcher:
clear_output()
else:
global_plot_prior_switcher = True
if what_to_plot == "Gaussian":
mu1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description="µ_prior", continuous_update=True)
mu2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description="µ_likelihood", continuous_update=True)
sigma1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_prior", continuous_update=True)
sigma2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_likelihood", continuous_update=True)
widget_ui = HBox([VBox([mu1_slider, sigma1_slider]),
VBox([mu2_slider, sigma2_slider])])
widget_out = interactive_output(plot_utility_gaussian,
{'mu1': mu1_slider,
'mu2': mu2_slider,
'sigma1': sigma1_slider,
'sigma2': sigma2_slider,
'mu_g': fixed(1.0),
'mu_c': fixed(-1.0),
'sigma_g': fixed(0.5),
'sigma_c': fixed(value=0.5),
'plot_utility_row': fixed(False)})
elif what_to_plot == "Uniform":
mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description="µ_likelihood", continuous_update=True)
sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_likelihood", continuous_update=True)
widget_ui = VBox([mu_slider, sigma_slider])
widget_out = interactive_output(plot_utility_uniform,
{'mu': mu_slider,
'sigma': sigma_slider,
'mu_g': fixed(1.0),
'mu_c': fixed(-1.0),
'sigma_g': fixed(0.5),
'sigma_c': fixed(value=0.5),
'plot_utility_row': fixed(False)})
elif what_to_plot == "Gamma":
alpha_slider = FloatSlider(min=1.0, max=10.0, step=0.1, value=2.0, description="α_prior", continuous_update=True)
beta_slider = FloatSlider(min=0.5, max=2.0, step=0.01, value=1.0, description="β_prior", continuous_update=True)
# offset_slider = FloatSlider(min=-6.0, max=2.0, step=0.1, value=0.0, description="offset", continuous_update=True)
offset_slider = fixed(0.0)
mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description="µ_likelihood", continuous_update=True)
sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_likelihood", continuous_update=True)
gaus_label = Label(value="normal likelihood", layout=Layout(display="flex", justify_content="center"))
gamma_label = Label(value="gamma prior", layout=Layout(display="flex", justify_content="center"))
widget_ui = HBox([VBox([gamma_label, alpha_slider, beta_slider]),
VBox([gaus_label, mu_slider, sigma_slider])])
widget_out = interactive_output(plot_utility_gamma,
{'alpha': alpha_slider,
'beta': beta_slider,
'offset': offset_slider,
'mu': mu_slider,
'sigma': sigma_slider,
'mu_g': fixed(1.0),
'mu_c': fixed(-1.0),
'sigma_g': fixed(0.5),
'sigma_c': fixed(value=0.5),
'plot_utility_row': fixed(False)})
elif what_to_plot == "Mixture of Gaussians":
mu_m1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description="µ_mix_p", continuous_update=True)
mu_m2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description="µ_mix_q", continuous_update=True)
sigma_m1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_mix_p", continuous_update=True)
sigma_m2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_mix_q", continuous_update=True)
factor_slider = FloatSlider(min=0.0, max=1.0, step=0.01, value=0.5, description="π", continuous_update=True)
mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description="µ_likelihood", continuous_update=True)
sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_likelihood", continuous_update=True)
widget_ui = HBox([VBox([mu_m1_slider, sigma_m1_slider, factor_slider]),
VBox([mu_m2_slider, sigma_m2_slider]),
VBox([mu_slider, sigma_slider])])
widget_out = interactive_output(plot_utility_mixture,
{'mu_m1': mu_m1_slider,
'mu_m2': mu_m2_slider,
'sigma_m1': sigma_m1_slider,
'sigma_m2': sigma_m2_slider,
'factor': factor_slider,
'mu': mu_slider,
'sigma': sigma_slider,
'mu_g': fixed(1.0),
'mu_c': fixed(-1.0),
'sigma_g': fixed(0.5),
'sigma_c': fixed(value=0.5),
'plot_utility_row': fixed(False)})
display(widget_ui, widget_out)
global global_plot_bayes_loss_utility_switcher
global_plot_bayes_loss_utility_switcher = False
def plot_bayes_loss_utility_switcher(what_to_plot):
global global_plot_bayes_loss_utility_switcher
if global_plot_bayes_loss_utility_switcher:
clear_output()
else:
global_plot_bayes_loss_utility_switcher = True
if what_to_plot == "Gaussian":
loss_f_options = Dropdown(
options=["Mean Squared Error", "Absolute Error", "Zero-One Loss"],
value="Mean Squared Error", description="Loss: ")
mu1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description="µ_prior", continuous_update=True)
sigma1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_prior", continuous_update=True)
mu2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description="µ_likelihood", continuous_update=True)
sigma2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_likelihood", continuous_update=True)
mu_true_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description="µ_true", continuous_update=True)
widget_ui = HBox([VBox([loss_f_options, mu1_slider, sigma1_slider]),
VBox([mu_true_slider, mu2_slider, sigma2_slider])])
widget_out = interactive_output(plot_bayes_loss_utility_gaussian,
{'mu1': mu1_slider,
'mu2': mu2_slider,
'sigma1': sigma1_slider,
'sigma2': sigma2_slider,
'mu_true': mu_true_slider,
'loss_f': loss_f_options})
elif what_to_plot == "Uniform":
loss_f_options = Dropdown(
options=["Mean Squared Error", "Absolute Error", "Zero-One Loss"],
value="Mean Squared Error", description="Loss: ")
mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description="µ_likelihood", continuous_update=True)
sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_likelihood", continuous_update=True)
mu_true_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description="µ_true", continuous_update=True)
widget_ui = HBox([VBox([loss_f_options, mu_slider, sigma_slider]),
VBox([mu_true_slider])])
widget_out = interactive_output(plot_bayes_loss_utility_uniform,
{'mu': mu_slider,
'sigma': sigma_slider,
'mu_true': mu_true_slider,
'loss_f': loss_f_options})
elif what_to_plot == "Gamma":
loss_f_options = Dropdown(
options=["Mean Squared Error", "Absolute Error", "Zero-One Loss"],
value="Mean Squared Error", description="Loss: ")
alpha_slider = FloatSlider(min=1.0, max=10.0, step=0.1, value=2.0, description="α_prior", continuous_update=True)
beta_slider = FloatSlider(min=0.5, max=2.0, step=0.01, value=1.0, description="β_prior", continuous_update=True)
# offset_slider = FloatSlider(min=-6.0, max=2.0, step=0.1, value=0.0, description="offset", continuous_update=True)
offset_slider = fixed(0.0)
mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description="µ_likelihood", continuous_update=True)
sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_likelihood", continuous_update=True)
mu_true_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description="µ_true", continuous_update=True)
gaus_label = Label(value="normal likelihood", layout=Layout(display="flex", justify_content="center"))
gamma_label = Label(value="gamma prior", layout=Layout(display="flex", justify_content="center"))
widget_ui = HBox([VBox([loss_f_options, gamma_label, alpha_slider, beta_slider]),
VBox([mu_true_slider, gaus_label, mu_slider, sigma_slider])])
widget_out = interactive_output(plot_bayes_loss_utility_gamma,
{'alpha': alpha_slider,
'beta': beta_slider,
'offset': offset_slider,
'mu': mu_slider,
'sigma': sigma_slider,
'mu_true': mu_true_slider,
'loss_f': loss_f_options})
elif what_to_plot == "Mixture of Gaussians":
mu_m1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description="µ_mix_p", continuous_update=True)
mu_m2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description="µ_mix_q", continuous_update=True)
sigma_m1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_mix_p", continuous_update=True)
sigma_m2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_mix_q", continuous_update=True)
factor_slider = FloatSlider(min=0.0, max=1.0, step=0.01, value=0.5, description="π", continuous_update=True)
mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5, description="µ_likelihood", continuous_update=True)
sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_likelihood", continuous_update=True)
mu_true_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5, description="µ_true", continuous_update=True)
loss_f_options = Dropdown(
options=["Mean Squared Error", "Absolute Error", "Zero-One Loss"],
value="Mean Squared Error", description="Loss: ")
empty_label = Label(value=" ")
widget_ui = HBox([VBox([loss_f_options, mu_m1_slider, sigma_m1_slider]),
VBox([mu_true_slider, mu_m2_slider, sigma_m2_slider]),
VBox([empty_label, mu_slider, sigma_slider])])
widget_out = interactive_output(plot_bayes_loss_utility_mixture,
{'mu_m1': mu_m1_slider,
'mu_m2': mu_m2_slider,
'sigma_m1': sigma_m1_slider,
'sigma_m2': sigma_m2_slider,
'factor': factor_slider,
'mu': mu_slider,
'sigma': sigma_slider,
'mu_true': mu_true_slider,
'loss_f': loss_f_options})
display(widget_ui, widget_out)
# @title Helper Functions
def gaussian(x, μ, σ):
""" Compute Gaussian probability density function for given value of the
random variable, mean, and standard deviation
Args:
x (scalar): value of random variable
μ (scalar): mean of Gaussian
σ (scalar): standard deviation of Gaussian
Returns:
scalar: value of probability density function
"""
return np.exp(-((x - μ) / σ)**2 / 2) / np.sqrt(2 * np.pi * σ**2)
def gamma_pdf(x, α, β):
return gamma_distribution.pdf(x, a=α, scale=1/β)
def mvn2d(x, y, mu1, mu2, sigma1, sigma2, cov12):
mvn = multivariate_normal([mu1, mu2], [[sigma1**2, cov12], [cov12, sigma2**2]])
return mvn.pdf(np.dstack((x, y)))
def product_guassian(mu1, mu2, sigma1, sigma2):
J_1, J_2 = 1/sigma1**2, 1/sigma2**2
J_3 = J_1 + J_2
mu_prod = (J_1*mu1/J_3) + (J_2*mu2/J_3)
sigma_prod = np.sqrt(1/J_3)
return mu_prod, sigma_prod
def reverse_product(mu3, sigma3, mu1, mu2):
J_3 = 1/sigma3**2
J_1 = J_3 * (mu3 - mu2) / (mu1 - mu2)
J_2 = J_3 * (mu3 - mu1) / (mu2 - mu1)
sigma1, sigma2 = 1/np.sqrt(J_1), 1/np.sqrt(J_2)
return sigma1, sigma2
def calc_mean_mode_median(x, y):
"""
"""
pdf = y * (x[1] - x[0])
# Calc mode of an arbitrary function
mode = x[np.argmax(pdf)]
# Calc mean of an arbitrary function
mean = np.multiply(x, pdf).sum()
# Calc median of an arbitrary function
cdf = np.cumsum(pdf)
idx = np.argmin(np.abs(cdf - 0.5))
median = x[idx]
return mean, median, mode
def calc_expected_loss(loss_f, posterior, x):
dx = x[1] - x[0]
expected_loss = np.zeros_like(x)
for i in np.arange(x.shape[0]):
loss = calc_loss_func(loss_f, x[i], x) # or mse or zero_one_loss
expected_loss[i] = np.sum(loss * posterior) * dx
return expected_loss
セクション0: はじめに
# @title Video 1: Introduction
from ipywidgets import widgets
from IPython.display import YouTubeVideo
from IPython.display import IFrame
from IPython.display import display
class PlayVideo(IFrame):
def __init__(self, id, source, page=1, width=400, height=300, **kwargs):
self.id = id
if source == 'Bilibili':
src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'
elif source == 'Osf':
src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'
super(PlayVideo, self).__init__(src, width, height, **kwargs)
def display_videos(video_ids, W=400, H=300, fs=1):
tab_contents = []
for i, video_id in enumerate(video_ids):
out = widgets.Output()
with out:
if video_ids[i][0] == 'Youtube':
video = YouTubeVideo(id=video_ids[i][1], width=W,
height=H, fs=fs, rel=0)
print(f'Video available at https://youtube.com/watch?v={video.id}')
else:
video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,
height=H, fs=fs, autoplay=False)
if video_ids[i][0] == 'Bilibili':
print(f'Video available at https://www.bilibili.com/video/{video.id}')
elif video_ids[i][0] == 'Osf':
print(f'Video available at https://osf.io/{video.id}')
display(video)
tab_contents.append(out)
return tab_contents
video_ids = [('Youtube', '7FTCHEEfAaA'), ('Bilibili', 'BV18P4y147nL')]
tab_contents = display_videos(video_ids, W=854, H=480)
tabs = widgets.Tab()
tabs.children = tab_contents
for i in range(len(tab_contents)):
tabs.set_title(i, video_ids[i][0])
display(tabs)
# @title Submit your feedback
content_review(f"{feedback_prefix}_Introduction_Video")
セクション1: アストロキャット!
あなたは宇宙飛行士の猫、アストロキャットだとしましょう!ジェットパックを使って宇宙を航行し、ネズミを追いかけるのが目標です。
猫なので、親指がなく自分でジェットパックを操作できません。地球の管制センターが操作します。
管制センターがあなたを導くためには、あなたの位置を知る必要があります。彼らはあなたの位置を推測しようとしています。彼らはあなたが宇宙ネズミの近くにいるのが好きだという事前知識を持っています。また、隠れた状態であるあなたの位置の信頼性の低い素早い測定も得ています。
彼らはベイズの定理とベイズ意思決定を使ってあなたの位置を推定しようとします。これからこのチュートリアルで詳しく見ていきます。
# @title Video 2: Astrocat
from ipywidgets import widgets
from IPython.display import YouTubeVideo
from IPython.display import IFrame
from IPython.display import display
class PlayVideo(IFrame):
def __init__(self, id, source, page=1, width=400, height=300, **kwargs):
self.id = id
if source == 'Bilibili':
src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'
elif source == 'Osf':
src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'
super(PlayVideo, self).__init__(src, width, height, **kwargs)
def display_videos(video_ids, W=400, H=300, fs=1):
tab_contents = []
for i, video_id in enumerate(video_ids):
out = widgets.Output()
with out:
if video_ids[i][0] == 'Youtube':
video = YouTubeVideo(id=video_ids[i][1], width=W,
height=H, fs=fs, rel=0)
print(f'Video available at https://youtube.com/watch?v={video.id}')
else:
video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,
height=H, fs=fs, autoplay=False)
if video_ids[i][0] == 'Bilibili':
print(f'Video available at https://www.bilibili.com/video/{video.id}')
elif video_ids[i][0] == 'Osf':
print(f'Video available at https://osf.io/{video.id}')
display(video)
tab_contents.append(out)
return tab_contents
video_ids = [('Youtube', 'N2u-1joYo6E'), ('Bilibili', 'BV14q4y1H7Bj')]
tab_contents = display_videos(video_ids, W=854, H=480)
tabs = widgets.Tab()
tabs.children = tab_contents
for i in range(len(tab_contents)):
tabs.set_title(i, video_ids[i][0])
display(tabs)
# @title Submit your feedback
content_review(f"{feedback_prefix}_Astrocat_Video")
アストロキャットは宇宙にいて、1次元の位置を考えています。隠れ状態 は真の位置です。衛星は潜在的な損失を、宇宙ネズミは利益を表します。間接的な測定を使って、管制センターとしてあなたはアストロキャットの位置を推定したり、どこに送るのが良さそうか決めたりできます。
この例では、あなた自身を科学者として考え、アストロキャットがどこにいると信じるか、誤差を考慮して点推定(単一の位置推定)を選ぶ方法、衛星と宇宙ネズミの位置に関する不確実性をどう扱うかを決める問題として考えられます。実際、遠隔衛星を制御する科学者が直面する問題と同じです!また、これは脳が動作目標やテニスボールを打つ目標を決めるときに行うこととも考えられます。多くの古典的実験がこの枠組みを使って、人間の意思決定や動作の最適性を研究しています。詳細は参考文献を参照してください。
セクション2: アストロキャットの位置の確率分布
チュートリアル開始からここまでの推定所要時間: 5分
まずは管制センターが位置をどう推定すべきか考えます。まだ測定は考慮せず、一般的にどのように不確実性を表現するかを考えます。今扱うのは連続分布で、アストロキャットの位置は任意の実数値を取れます。前回のチュートリアルでは離散分布で、魚は片側かもう片側かでした。
では、アストロキャットがいる可能性のある無限の点それぞれの確率をどう表現するでしょうか?
ベイズアプローチはどんな確率分布にも使えます。多くの変数は複雑または未知(例えば経験的)な分布を使う必要がありますが、ここではガウス分布やその拡張を使います。
セクション2.1: ガウス分布
# @title Video 3: The Gaussian distribution
from ipywidgets import widgets
from IPython.display import YouTubeVideo
from IPython.display import IFrame
from IPython.display import display
class PlayVideo(IFrame):
def __init__(self, id, source, page=1, width=400, height=300, **kwargs):
self.id = id
if source == 'Bilibili':
src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'
elif source == 'Osf':
src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'
super(PlayVideo, self).__init__(src, width, height, **kwargs)
def display_videos(video_ids, W=400, H=300, fs=1):
tab_contents = []
for i, video_id in enumerate(video_ids):
out = widgets.Output()
with out:
if video_ids[i][0] == 'Youtube':
video = YouTubeVideo(id=video_ids[i][1], width=W,
height=H, fs=fs, rel=0)
print(f'Video available at https://youtube.com/watch?v={video.id}')
else:
video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,
height=H, fs=fs, autoplay=False)
if video_ids[i][0] == 'Bilibili':
print(f'Video available at https://www.bilibili.com/video/{video.id}')
elif video_ids[i][0] == 'Osf':
print(f'Video available at https://osf.io/{video.id}')
display(video)
tab_contents.append(out)
return tab_contents
video_ids = [('Youtube', 'w76wZWrDMNE'), ('Bilibili', 'BV1RX4y1F7G7')]
tab_contents = display_videos(video_ids, W=854, H=480)
tabs = widgets.Tab()
tabs.children = tab_contents
for i in range(len(tab_contents)):
tabs.set_title(i, video_ids[i][0])
display(tabs)
ビデオのテキスト要約はこちらをクリック
このチュートリアルで使う分布の一つはガウス分布(正規分布とも呼ばれます)です。
これは特別でよく使われる分布で、いくつか理由があります。統計学で最も重要な定理の一つ、中心極限定理の焦点でもあります。この定理は、ある変数のサンプルを多数足し合わせると、その和は元の変数の分布に関わらず正規分布に近づくことを示します。これは今は詳しく扱いませんが、ボーナスのリンクを参照してください。さらに、ガウス分布は数学的に扱いやすく、重要な問題に対して簡単な閉形式解を提供します。後で見るように、ガウス分布は混合ガウス分布として拡張でき、他の分布をよく近似できます。要するに、ガウス分布は連続分布で最も重要なものの一つです。
ガウス分布は2つのパラメータを持ちます。平均 は中心の位置を決めます。広がりは標準偏差 またはその二乗の分散 で制御されます。標準偏差と分散は混同しやすいので注意してください。
# @title Submit your feedback
content_review(f"{feedback_prefix}_The_Gaussian_distribution_Video")
ガウス分布の変数 に対する式は次の通りです:
この例では、 はアストロキャットの1次元の位置です。 はNormal(ガウス)分布を表す標準的な記法です。例えば、 は平均0、分散1のガウス分布を意味します。この式の形は直感的ではありませんが、平均と標準偏差の値が確率分布にどう影響するかを見ていきます。
ここではコードでガウス分布を実装しませんが、必要なら前提知識の復習 W0D5 T1$ を参照してください。
# @markdown Execute this cell to enable the function `gaussian`
def gaussian(x, mu, sigma):
return np.exp(-(x - mu)**2 / (2*sigma**2)) / np.sqrt(2 * np.pi * sigma**2)
インタラクティブデモ 2.1: ガウス分布のパラメータを探る
ガウス分布のパラメータが分布にどう影響するかを探ります。以下のデモで平均 と標準偏差 を変えてみてください。
以下について議論しましょう:
- を増やすとどうなる? を増やすとどうなる?
- ある事象が で起こる確率を知りたい場合、同じ確率を生む異なる と の組み合わせは2つ見つけられますか?
- で同じ確率を生むガウス分布はいくつ作れますか?
# @markdown Execute this cell to enable the widget
widget = interact(plot_gaussian,
mu=FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.0,
continuous_update=False, description="µ"),
sigma=FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5,
continuous_update=False,description="σ"))
# @title Submit your feedback
content_review(f"{feedback_prefix}_Exploring_Gaussian_parameters_Interactive_Demo_and_Discussion")
セクション2.2: ガウス分布の掛け算
# @title Video 4: Multiplying Gaussians
from ipywidgets import widgets
from IPython.display import YouTubeVideo
from IPython.display import IFrame
from IPython.display import display
class PlayVideo(IFrame):
def __init__(self, id, source, page=1, width=400, height=300, **kwargs):
self.id = id
if source == 'Bilibili':
src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'
elif source == 'Osf':
src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'
super(PlayVideo, self).__init__(src, width, height, **kwargs)
def display_videos(video_ids, W=400, H=300, fs=1):
tab_contents = []
for i, video_id in enumerate(video_ids):
out = widgets.Output()
with out:
if video_ids[i][0] == 'Youtube':
video = YouTubeVideo(id=video_ids[i][1], width=W,
height=H, fs=fs, rel=0)
print(f'Video available at https://youtube.com/watch?v={video.id}')
else:
video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,
height=H, fs=fs, autoplay=False)
if video_ids[i][0] == 'Bilibili':
print(f'Video available at https://www.bilibili.com/video/{video.id}')
elif video_ids[i][0] == 'Osf':
print(f'Video available at https://osf.io/{video.id}')
display(video)
tab_contents.append(out)
return tab_contents
video_ids = [('Youtube', 'tmi3JwbvmpY'), ('Bilibili', 'BV1pM4y1T78x')]
tab_contents = display_videos(video_ids, W=854, H=480)
tabs = widgets.Tab()
tabs.children = tab_contents
for i in range(len(tab_contents)):
tabs.set_title(i, video_ids[i][0])
display(tabs)
このビデオでは2つのガウス分布の掛け算を扱います。
ビデオのテキスト要約はこちらをクリック
ガウス分布の掛け算は、ランダム変数の掛け算ではなく、分布自体の掛け算です。平均 , 、標準偏差 , の2つのガウス分布を掛けると、またガウス分布になります。掛け算の結果のガウス分布の平均 と標準偏差 は次のようになります:
\begin{align}
&= a \
&= \
a &= \frac{\sigma_{1}^{-2}}{\sigma_{1}^{-2} + \sigma_{2}^{-2}}
\end{align}
これは少し複雑に見えますが、ガウス分布の情報量は分散の逆数 で表されます。つまり、ガウス分布の掛け算では、結果の平均は元の平均の情報量に比例した重み付き平均で、結果の情報量は元の2つの情報量の和になります。次のデモで詳しく見ていきます。
# @title Submit your feedback
content_review(f"{feedback_prefix}_Multiplying_Gaussians_Video")
インタラクティブデモ 2.2: ガウス分布の掛け算
2つのガウス分布の掛け算を実装しました。以下のウィジェットを使って、2つのガウス分布の情報量と組み合わせを考えます。ここでは衛星と宇宙ネズミの間の中間位置を見つけたいと想像してください。これは2つの位置の中心(平均)です。不確実性があるので、最もありそうな場所を考える際に不確実性を重み付けする必要があります。
このデモでは、 と は衛星位置の分布の平均と標準偏差、 と は宇宙ネズミ位置の分布の平均と標準偏差、 と は掛け算によって得られた中心位置の分布の平均と標準偏差です。
質問:
- のとき、 に関してどれくらいの不確実性(情報量)がありますか?
- のとき、 の推定はどうなりますか?(この場合、 は最大11までですが十分大きいです)
- の場合、 の推定はどう違いますか?最初の例と何が変わりましたか?
- に設定し、 が約2になるように を変えてみてください。同じ を生む はいくつありますか?
- 続けて、 に設定した場合、 にするにはどの を変える必要がありますか?
- のとき、平均に関してどれくらいの情報がありますか?
# @markdown Execute this cell to enable the widget
mu1_slider = FloatSlider(min=-5.0, max=-0.51, step=0.01, value=-2.0,
description="µ_1", continuous_update=True)
mu2_slider = FloatSlider(min=0.5, max=5.01, step=0.01, value=2.0,
description="µ_2",continuous_update=True)
sigma1_slider = FloatSlider(min=0.1, max=11.01, step=0.01, value=1.0,
description="σ_1", continuous_update=True)
sigma2_slider = FloatSlider(min=0.1, max=11.01, step=0.01, value=1.0,
description="σ_2", continuous_update=True)
distro_1_label = Label(value="Satellite",
layout=Layout(display="flex", justify_content="center"))
distro_2_label = Label(value="Space Mouse",
layout=Layout(display="flex", justify_content="center"))
widget_ui = HBox([VBox([distro_1_label, mu1_slider, sigma1_slider]),
VBox([distro_2_label, mu2_slider, sigma2_slider])])
widget_out = interactive_output(plot_information, {'mu1': mu1_slider,
'mu2': mu2_slider,
'sigma1': sigma1_slider,
'sigma2': sigma2_slider})
display(widget_ui, widget_out)
# @title Submit your feedback
content_review(f"{feedback_prefix}_Multiplying_Gaussians_Interactive_Demo_and_Discussion")
これらの概念をシステム神経科学で直接使うことを考えると、感覚空間の位置を表す2つのニューロンの応答を組み合わせ(平均化)したときにどれだけ情報が得られるか(受容野がどれだけ共有されているか)を知りたい場合があります。これはガウス分布の掛け算に相当します!
セクション2.3: ガウス混合分布
# @title Video 5: Mixtures of Gaussians
from ipywidgets import widgets
from IPython.display import YouTubeVideo
from IPython.display import IFrame
from IPython.display import display
class PlayVideo(IFrame):
def __init__(self, id, source, page=1, width=400, height=300, **kwargs):
self.id = id
if source == 'Bilibili':
src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'
elif source == 'Osf':
src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'
super(PlayVideo, self).__init__(src, width, height, **kwargs)
def display_videos(video_ids, W=400, H=300, fs=1):
tab_contents = []
for i, video_id in enumerate(video_ids):
out = widgets.Output()
with out:
if video_ids[i][0] == 'Youtube':
video = YouTubeVideo(id=video_ids[i][1], width=W,
height=H, fs=fs, rel=0)
print(f'Video available at https://youtube.com/watch?v={video.id}')
else:
video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,
height=H, fs=fs, autoplay=False)
if video_ids[i][0] == 'Bilibili':
print(f'Video available at https://www.bilibili.com/video/{video.id}')
elif video_ids[i][0] == 'Osf':
print(f'Video available at https://osf.io/{video.id}')
display(video)
tab_contents.append(out)
return tab_contents
video_ids = [('Youtube', '2mv-5qFYV5o'), ('Bilibili', 'BV17B4y1K7En')]
tab_contents = display_videos(video_ids, W=854, H=480)
tabs = widgets.Tab()
tabs.children = tab_contents
for i in range(len(tab_contents)):
tabs.set_title(i, video_ids[i][0])
display(tabs)
ビデオのテキスト要約はこちらをクリック
もし連続分布が単一の山でうまく表せない場合はどうでしょう?例えば、アストロキャットがよくある場所が2か所ある場合、ガウス分布ではうまく表現できません。多峰性分布が必要です。幸い、ガウス混合分布としてガウス分布を拡張できます。
ガウス混合分布では、複数の標準ガウス分布を重み付きで足し合わせ(正規化して積分が1になるように)、より複雑な分布を作ります。各ガウス分布は通常通り平均と標準偏差で表されます。混合分布の追加パラメータは各ガウスの重み です。以下のデモでガウス混合分布と標準ガウス成分の関係を理解しましょう。導出は扱いませんが、ボーナス課題として挑戦できます。
混合分布はベイズモデリングでよく使われる重要なツールです。
# @title Submit your feedback
content_review(f"{feedback_prefix}_Mixtures_of_Gaussians_Video")
インタラクティブデモ 2.3: ガウス混合分布を探る
2つのガウス分布の混合を調べます。重みパラメータ が1つあり、1つのガウス分布の重みを決めます。もう一方は です。
以下のウィジェットで各ガウスのパラメータと混合重み () を操作し、混合ガウス分布の挙動を理解しましょう。
以下の質問について議論してください:
- 重み を増やすと混合分布(濃い青)はどう変わる?
- 混合分布の2つの山をもっと離すには?
- 混合分布を単一の山(ガウスのように)にできますか?
- 混合分布で単一の丸い山や2つの分離した山以外の形は作れますか?
# @markdown Execute this cell to enable the widget
mu1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-1.0,
description="µ_p", continuous_update=True)
mu2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=1.0,
description="µ_q", continuous_update=True)
sigma1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5,
description="σ_p", continuous_update=True)
sigma2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5,
description="σ_q", continuous_update=True)
factor_slider = FloatSlider(min=0.1, max=0.9, step=0.01, value=0.5,
description="π", continuous_update=True)
distro_1_label = Label(value="p(x)",
layout=Layout(display="flex", justify_content="center"))
distro_2_label = Label(value="q(x)",
layout=Layout(display="flex", justify_content="center"))
mixture_label = Label(value="mixing coefficient",
layout=Layout(display="flex", justify_content="center"))
widget_ui = HBox([VBox([distro_1_label, mu1_slider, sigma1_slider]),
VBox([distro_2_label, mu2_slider, sigma2_slider]),
VBox([mixture_label, factor_slider])])
widget_out = interactive_output(gaussian_mixture, {'mu1': mu1_slider,
'mu2': mu2_slider,
'sigma1': sigma1_slider,
'sigma2': sigma2_slider,
'factor': factor_slider})
display(widget_ui, widget_out)
# @title Submit your feedback
content_review(f"{feedback_prefix}_Exploring_Gaussian_mixtures_Interactive_Demo_and_Discussion")
セクション3: 効用
チュートリアル開始からここまでの推定所要時間: 25分
# @title Video 6: Utility Loss Estimators
from ipywidgets import widgets
from IPython.display import YouTubeVideo
from IPython.display import IFrame
from IPython.display import display
class PlayVideo(IFrame):
def __init__(self, id, source, page=1, width=400, height=300, **kwargs):
self.id = id
if source == 'Bilibili':
src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'
elif source == 'Osf':
src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'
super(PlayVideo, self).__init__(src, width, height, **kwargs)
def display_videos(video_ids, W=400, H=300, fs=1):
tab_contents = []
for i, video_id in enumerate(video_ids):
out = widgets.Output()
with out:
if video_ids[i][0] == 'Youtube':
video = YouTubeVideo(id=video_ids[i][1], width=W,
height=H, fs=fs, rel=0)
print(f'Video available at https://youtube.com/watch?v={video.id}')
else:
video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,
height=H, fs=fs, autoplay=False)
if video_ids[i][0] == 'Bilibili':
print(f'Video available at https://www.bilibili.com/video/{video.id}')
elif video_ids[i][0] == 'Osf':
print(f'Video available at https://osf.io/{video.id}')
display(video)
tab_contents.append(out)
return tab_contents
video_ids = [('Youtube', 'M_pIY0xv_xw'), ('Bilibili', 'BV1zP4y147bK')]
tab_contents = display_videos(video_ids, W=854, H=480)
tabs = widgets.Tab()
tabs.children = tab_contents
for i in range(len(tab_contents)):
tabs.set_title(i, video_ids[i][0])
display(tabs)
ビデオのテキスト要約はこちらをクリック
アストロキャットの位置を知りたいです。管制センターに座標を提供するとき、確率分布全体を示すわけにはいきません。単一の座標を推定して示しますが、その推定は簡単に最もありそうな値を取るだけではありません。推定した位置と実際の位置が違う場合の良し悪しを考えます。
効用は、ある行動を取ったときに隠れ状態の値に対して得られる利益(負なら損失)を表します。連続例では、状態のすべての連続値に対して効用を定義できる関数が必要です。ここでの行動はアストロキャットの位置の推定値です。
損失関数は、真の値 に対して推定値 を選んだときの「コスト」(誤差の罰則)を決めます。損失関数は効用関数の負の値に相当します。
# @title Submit your feedback
content_review(f"{feedback_prefix}_Utility_loss_estimators_Video")
セクション3.1: 標準的な損失関数
多くの損失関数がありますが、ここでは3つに注目します。平均二乗誤差(真値と推定値の差の二乗)、絶対誤差(差の絶対値)、ゼロ・ワン損失(推定値が真値と完全に一致しなければ損失1)です。以下の式で表されます:
\begin{eqnarray}
&=& ( \
&=& \
&=& \begin{cases}
0,& \
1, & \textrm{それ以外}
\end{cases}
\end{eqnarray}
これらの損失関数が期待効用にどう影響するかを探ります。
次のセルで calc_loss_func 関数の実装を見てください。
# @markdown Execute this cell to enable the function `calc_loss_func`
def calc_loss_func(loss_f, mu_true, x):
error = x - mu_true
if loss_f == "Mean Squared Error":
loss = (error)**2
elif loss_f == "Absolute Error":
loss = np.abs(error)
elif loss_f == "Zero-One Loss":
loss = (np.abs(error) >= 0.03).astype(float)
return loss
インタラクティブデモ 3: 異なる分布で損失を探る
損失関数が確率分布とどう相互作用して期待効用や行動に影響するか見てみましょう。
以下のウィジェットを操作し、次のことを議論してください:
- ガウス分布の場合、3つの損失関数で期待効用のピークはx軸上で変わりますか?このピークは選択する行動(推定位置)を示します。つまり、損失関数の違いで行動は変わりますか?
- 2つの山を持つガウス混合分布の場合、期待損失のピークはx軸上で変わりますか?
- 平均、最頻値、中央値がすべて異なる(等しくない)ガウス混合分布のパラメータを見つけてください。この分布で、3つの損失関数それぞれの期待効用のピークは確率分布の平均・中央値・最頻値とどう対応しますか?
- 2つの山が同じ高さのガウス混合分布では、モードはいくつありますか?
# @markdown Execute this cell to enable the widget
widget = interact(loss_plot_switcher,
what_to_plot = Dropdown(
options=["Gaussian", "Mixture of Gaussians"],
value="Gaussian", description="Distribution: "))
# @title Submit your feedback
content_review(f"{feedback_prefix}_Exploring_loss_with_different_distributions_Interactive_Demo_and_Discussion")
アストロキャットの座標を推定するのは確率分布だけでは簡単ではないことがわかります。推定には効用/損失の概念と特定の損失関数が必要です。
対称分布では平均、中央値、最頻値は同じですが、ガンマ分布や指数分布のような歪みのある分布では異なります。以下で事前分布としてさらに多くの分布を探れます。
セクション3.2: より複雑な損失関数
先ほどの損失関数は比較的単純でよく使われますが、現実は複雑です。ここではアストロキャットは宇宙ネズミの近くにいることを望み、衛星からは離れていたいので、これを反映した複雑な損失関数が必要です。
損失関数ではなく、アストロキャットにとって重要な利益と損失を考慮した一般化効用関数を考えます。
不確実性に応じて、アストロキャットが「良い」空間の近くにいると推定し、「悪い」空間から離れていると推定したいと考えます。これらの効用はガウス分布の利益・損失領域としてモデル化し、幅は宇宙ネズミと衛星の位置の不確実性から来ると仮定します。
次のインタラクティブデモでこれを探りましょう。
インタラクティブデモ 3.2: 複雑な猫のコスト
損失関数を使って推定量や誤差に基づく期待損失を決める方法を探ったので、今度は一般化効用関数を使った場合の推定値の変化を見ます。
質問:
- アストロキャットの を変えると期待効用はどう変わる?
- 期待効用がどこでも正確にゼロになることはある?
- アストロキャットの周囲の領域で期待効用がゼロで、他の場所で正負に分かれることはある?
- アストロキャットの位置の不確実性が増すと期待効用はどうなる?
# @markdown Execute this cell to enable the widget
mu_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.0, description="µ", continuous_update=True)
mu_g_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-1.0, description="µ_gain", continuous_update=True)
mu_c_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=1.0, description="µ_cost", continuous_update=True)
sigma_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ", continuous_update=True)
sigma_g_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_gain", continuous_update=True)
sigma_c_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5, description="σ_cost", continuous_update=True)
distro_label = Label(value="probability", layout=Layout(display="flex", justify_content="center"))
gain_label = Label(value="gain", layout=Layout(display="flex", justify_content="center"))
loss_label = Label(value="loss", layout=Layout(display="flex", justify_content="center"))
widget_ui = HBox([VBox([distro_label, mu_slider, sigma_slider]),
VBox([gain_label, mu_g_slider, sigma_g_slider]),
VBox([loss_label, mu_c_slider, sigma_c_slider])])
widget_out = interactive_output(plot_simple_utility_gaussian,
{'mu': mu_slider,
'sigma': sigma_slider,
'mu_g': mu_g_slider,
'mu_c': mu_c_slider,
'sigma_g': sigma_g_slider,
'sigma_c': sigma_c_slider})
display(widget_ui, widget_out)
# @title Submit your feedback
content_review(f"{feedback_prefix}_Complicated_cat_costs_Interactive_Demo_and_Discussion")
セクション4: 相関と周辺化
チュートリアル開始からここまでの推定所要時間: 40分
このセクションでは2次元ガウス分布を扱います。これは2つのガウス確率変数の結合分布です。
# @title Video 7: Correlation and marginalization
from ipywidgets import widgets
from IPython.display import YouTubeVideo
from IPython.display import IFrame
from IPython.display import display
class PlayVideo(IFrame):
def __init__(self, id, source, page=1, width=400, height=300, **kwargs):
self.id = id
if source == 'Bilibili':
src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'
elif source == 'Osf':
src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'
super(PlayVideo, self).__init__(src, width, height, **kwargs)
def display_videos(video_ids, W=400, H=300, fs=1):
tab_contents = []
for i, video_id in enumerate(video_ids):
out = widgets.Output()
with out:
if video_ids[i][0] == 'Youtube':
video = YouTubeVideo(id=video_ids[i][1], width=W,
height=H, fs=fs, rel=0)
print(f'Video available at https://youtube.com/watch?v={video.id}')
else:
video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,
height=H, fs=fs, autoplay=False)
if video_ids[i][0] == 'Bilibili':
print(f'Video available at https://www.bilibili.com/video/{video.id}')
elif video_ids[i][0] == 'Osf':
print(f'Video available at https://osf.io/{video.id}')
display(video)
tab_contents.append(out)
return tab_contents
video_ids = [('Youtube', 'KbtmWhGlBYg'), ('Bilibili', 'BV1Eg411M72d')]
tab_contents = display_videos(video_ids, W=854, H=480)
tabs = widgets.Tab()
tabs.children = tab_contents
for i in range(len(tab_contents)):
tabs.set_title(i, video_ids[i][0])
display(tabs)
# @title Submit your feedback
content_review(f"{feedback_prefix}_Correlation_and_marginalization_Video")
セクション4.1: 相関
関連部分のビデオテキスト要約はこちらをクリック
2変数の2次元ガウスで変数が独立なら、片方を見てももう片方について何もわかりません。しかし相関(共分散)がある場合はどうでしょう?
平均 , 、標準偏差 , の2つのガウスの共分散は:
ここで は期待値を示します。つまり、 の値から平均を引いたものと の値から平均を引いたものの積の期待値です。
相関は共分散を正規化したもので、-1(完全な逆相関)から1(完全な相関)までの値を取ります。
これは重要な概念で、2つの隠れ状態(または確率変数)だけでなく、次元のガウス確率変数ベクトルにも拡張されます。計算論的神経科学で広く使われています。
インタラクティブデモ 4.1: 共変する2Dガウス
この2Dガウス(2つのガウスの結合分布)を探りましょう。
以下のウィジェットで次の質問を考えてください:
- これらの変数が関心のある隠れ状態を表す場合、片方を観測するともう片方について何がわかる?相関にどう依存する?
- 平均、分散、相関を変えると分布の形はどう変わる?
- 片方の隠れ状態分布だけを取り出したい場合はどうする?(ヒント: チュートリアル1を思い出してください)
# @markdown Execute the cell to enable the widget
mu1_slider = FloatSlider(min=-1.0, max=1.0, step=0.01, value=0.0,
description="µ_x", continuous_update=True)
mu2_slider = FloatSlider(min=-1.0, max=1.0, step=0.01, value=0.0,
description="µ_y", continuous_update=True)
sigma1_slider = FloatSlider(min=0.1, max=1.5, step=0.01, value=0.5,
description="σ_x", continuous_update=True)
sigma2_slider = FloatSlider(min=0.1, max=1.5, step=0.01, value=0.5,
description="σ_y", continuous_update=True)
corr_slider = FloatSlider(min=-0.99, max=0.99, step=0.01, value=0.0,
description="ρ", continuous_update=True)
distro1_label = Label(value="p(x)",
layout=Layout(display="flex", justify_content="center"))
distro2_label = Label(value="p(y)",
layout=Layout(display="flex", justify_content="center"))
corr_label = Label(value="correlation",
layout=Layout(display="flex", justify_content="center"))
widget_ui = HBox([VBox([distro1_label, mu1_slider, sigma1_slider]),
VBox([corr_label, corr_slider]),
VBox([distro2_label, mu2_slider, sigma2_slider])])
widget_out = interactive_output(plot_mvn2d, {'mu1': mu1_slider,
'mu2': mu2_slider,
'sigma1': sigma1_slider,
'sigma2': sigma2_slider,
'corr': corr_slider})
display(widget_ui, widget_out)
# @title Submit your feedback
content_review(f"{feedback_prefix}_Covarying_2D_Gaussian_Interactive_Demo_and_Discussion")
セクション4.2: 周辺化と情報量
関連部分のビデオテキスト要約はこちらをクリック
チュートリアル1で、ある変数の確率を求めるには他の変数を平均化(周辺化)する必要があると学びました。相関のあるガウス分布でも同様です。例えば、2つの変数がアストロキャットの空間上の位置(2次元)を表すとします。XかYの不確実性を得たい場合は周辺化します。
しかし、Xの測定値があってYの不確実性を知りたい場合は周辺化ではなく条件付き確率 を計算します。次のインタラクティブデモでこの関係を探ります。
また、不確実性は情報量の逆数として考えられます。結合情報量は相関によって決まります。ベイズアプローチでは、事前分布と尤度の相互情報量として考えることも重要です。
インタラクティブデモ 4.2: 2Dガウスの周辺化
以下のウィジェットで次の質問を考えてください:
- いつ周辺分布と条件付き確率分布は同じになる?なぜ?
- が大きいとき、両方の変数を見ることでどれだけ情報が増える?
- が0に近いが2つの変数の分散が大きく異なる場合、条件付き確率は周辺分布とどう違う? が変わると?
# @markdown Execute this cell to enable the widget
sigma1_slider = FloatSlider(min=0.1, max=1.1, step=0.01, value=0.5,
description="σ_x", continuous_update=True)
sigma2_slider = FloatSlider(min=0.1, max=1.1, step=0.01, value=0.5,
description="σ_y", continuous_update=True)
c_x_slider = FloatSlider(min=-1.0, max=1.0, step=0.01, value=0.0,
description="Cx", continuous_update=True)
c_y_slider = FloatSlider(min=-1.0, max=1.0, step=0.01, value=0.0,
description="Cy", continuous_update=True)
corr_slider = FloatSlider(min=-1.0, max=1.0, step=0.01, value=0.0,
description="ρ", continuous_update=True)
distro1_label = Label(value="x",
layout=Layout(display="flex", justify_content="center"))
distro2_label = Label(value="y",
layout=Layout(display="flex", justify_content="center"))
corr_label = Label(value="correlation",
layout=Layout(display="flex", justify_content="center"))
widget_ui = HBox([VBox([distro1_label, sigma1_slider, c_x_slider]),
VBox([corr_label, corr_slider]),
VBox([distro2_label, sigma2_slider, c_y_slider])])
widget_out = interactive_output(plot_marginal, {'sigma1': sigma1_slider,
'sigma2': sigma2_slider,
'c_x': c_x_slider,
'c_y': c_y_slider,
'corr': corr_slider})
display(widget_ui, widget_out)
# @title Submit your feedback
content_review(f"{feedback_prefix}_Marginalization_and_information_Interactive_Demo_and_Discussion")
セクション5: 連続分布に対するベイズの定理
チュートリアル開始からここまでの推定所要時間: 55分
# @title Video 8: Posterior Beliefs
from ipywidgets import widgets
from IPython.display import YouTubeVideo
from IPython.display import IFrame
from IPython.display import display
class PlayVideo(IFrame):
def __init__(self, id, source, page=1, width=400, height=300, **kwargs):
self.id = id
if source == 'Bilibili':
src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'
elif source == 'Osf':
src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'
super(PlayVideo, self).__init__(src, width, height, **kwargs)
def display_videos(video_ids, W=400, H=300, fs=1):
tab_contents = []
for i, video_id in enumerate(video_ids):
out = widgets.Output()
with out:
if video_ids[i][0] == 'Youtube':
video = YouTubeVideo(id=video_ids[i][1], width=W,
height=H, fs=fs, rel=0)
print(f'Video available at https://youtube.com/watch?v={video.id}')
else:
video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,
height=H, fs=fs, autoplay=False)
if video_ids[i][0] == 'Bilibili':
print(f'Video available at https://www.bilibili.com/video/{video.id}')
elif video_ids[i][0] == 'Osf':
print(f'Video available at https://osf.io/{video.id}')
display(video)
tab_contents.append(out)
return tab_contents
video_ids = [('Youtube', '84Fqve5i28U'), ('Bilibili', 'BV1xU4y1n7SX')]
tab_contents = display_videos(video_ids, W=854, H=480)
tabs = widgets.Tab()
tabs.children = tab_contents
for i in range(len(tab_contents)):
tabs.set_title(i, video_ids[i][0])
display(tabs)
# @title Submit your feedback
content_review(f"{feedback_prefix}_Posterior_beliefs_Video")
連続の場合、事後分布の形が事前分布とどう異なるかを考えられます。ガウスの場合が基本ですが、非対称な事前分布(または尤度)や事後分布では、ベイズの定理適用時に平均、中央値、最頻値がどう影響を受けるか見られます。
セクション5.1: ガウスの例
ベイズの定理は2つの情報源を結合する方法を示します。事前分布(例: 管制センターのアストロキャット位置のノイズのある期待)と尤度(例: 測定後のアストロキャットのノイズのある表現)を組み合わせて、両方の情報を考慮した事後分布(信念分布)を得ます。ベイズの定理は次の通りです:
事前と尤度がガウス分布の場合、ベイズの定理は次の形になります:
\begin{align}
&= \
&= \
&= \
&\propto \mathcal{N}(\mu_{likelihood},\sigma_{likelihood}^2) \times \mathcal{N}(\mu_{prior},\sigma_{prior}^2)
\end{align}
事後分布のパラメータは、セクション2.2でのガウス分布の掛け算と同様に得られます。
インタラクティブデモ 5.1: ガウスのベイズ
以下のインタラクティブデモで次の質問を考えましょう:
- ガウスの事後分布では情報がどう結合されているように見えますか?(ヒント: 事前演習を思い出してください)
- ここでの事後分布と前の演習での2つのガウスの平均を表すガウスの違いは?
- 事前分布と尤度の情報の相対的な重み付けはどう考えればよい?
# @markdown Execute this cell to enable the widget
mu1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5,
description="µ_prior", continuous_update=True)
mu2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5,
description="µ_likelihood", continuous_update=True)
sigma1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5,
description="σ_prior", continuous_update=True)
sigma2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5,
description="σ_likelihood", continuous_update=True)
distro1_label = Label(value="prior distribution",
layout=Layout(display="flex", justify_content="center"))
distro2_label = Label(value="likelihood distribution",
layout=Layout(display="flex", justify_content="center"))
widget_ui = HBox([VBox([distro1_label, mu1_slider, sigma1_slider]),
VBox([distro2_label, mu2_slider, sigma2_slider])])
widget_out = interactive_output(plot_bayes, {'mu1': mu1_slider,
'mu2': mu2_slider,
'sigma1': sigma1_slider,
'sigma2': sigma2_slider})
display(widget_ui, widget_out)
セクション5.2: 事前分布の探求
インタラクティブデモ 5.2: 事前分布の探求
もしアストロキャットの位置の事前分布がガウスでなかったらどうなるでしょう?ベイズの定理は同じですが、解析解はより複雑か不可能かもしれません。異なる事前分布のときの事後分布の挙動を見てみましょう。
以下の質問を考えてください:
- 非ガウスの事前分布を使うと事後分布がガウスに見えないのはなぜ?
- 平坦な事前分布とは何を意味する?
- ガンマ分布の事前分布は他とどう違う?
- 尤度がガウス以外になることは想像できますか?
# @markdown Execute this cell to enable the widget
widget = interact(plot_prior_switcher,
what_to_plot = Dropdown(
options=["Gaussian", "Mixture of Gaussians",
"Uniform", "Gamma"],
value="Gaussian", description="Prior: "))
# @title Submit your feedback
content_review(f"{feedback_prefix}_Prior_exploration_Interactive_Demo_and_Discussion")
セクション6: ベイズ意思決定
チュートリアル開始からここまでの推定所要時間: 1時間15分
# @title Video 9: Bayesian Decisions
from ipywidgets import widgets
from IPython.display import YouTubeVideo
from IPython.display import IFrame
from IPython.display import display
class PlayVideo(IFrame):
def __init__(self, id, source, page=1, width=400, height=300, **kwargs):
self.id = id
if source == 'Bilibili':
src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'
elif source == 'Osf':
src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'
super(PlayVideo, self).__init__(src, width, height, **kwargs)
def display_videos(video_ids, W=400, H=300, fs=1):
tab_contents = []
for i, video_id in enumerate(video_ids):
out = widgets.Output()
with out:
if video_ids[i][0] == 'Youtube':
video = YouTubeVideo(id=video_ids[i][1], width=W,
height=H, fs=fs, rel=0)
print(f'Video available at https://youtube.com/watch?v={video.id}')
else:
video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,
height=H, fs=fs, autoplay=False)
if video_ids[i][0] == 'Bilibili':
print(f'Video available at https://www.bilibili.com/video/{video.id}')
elif video_ids[i][0] == 'Osf':
print(f'Video available at https://osf.io/{video.id}')
display(video)
tab_contents.append(out)
return tab_contents
video_ids = [('Youtube', 'lqh41KXtwL0'), ('Bilibili', 'BV1Q54y1E7eE')]
tab_contents = display_videos(video_ids, W=854, H=480)
tabs = widgets.Tab()
tabs.children = tab_contents
for i in range(len(tab_contents)):
tabs.set_title(i, video_ids[i][0])
display(tabs)
# @title Submit your feedback
content_review(f"{feedback_prefix}_Bayesian_Decisions_Video")
連続次元でのベイズ意思決定は、チュートリアル1の二値の場合と同じです。違いは期待効用の計算が積分で行われ、確率分布がすべて連続であることです。
セクション6.1: 事後分布に対するベイズ推定
事後分布がガウス以外の形を取ることがあると理解したので、損失関数を再検討します。事後分布は多様な形を取れます。
インタラクティブデモ 6.1: さまざまな事前分布での標準損失関数
質問:
- 双峰性の事前分布がある場合、異なる損失関数は学習内容にどう影響する?
- 損失関数が事後分布の形に対して異なる振る舞いをするのはなぜ?いつ異なる期待損失を生む?
- ガウス混合分布の場合、期待損失がガウスの場合と異なる状況を説明してください。
# @markdown Execute this cell to enable the widget
widget = interact(plot_bayes_loss_utility_switcher,
what_to_plot = Dropdown(
options=["Gaussian", "Mixture of Gaussians",
"Uniform", "Gamma"],
value="Gaussian", description="Prior: "))
# @title Submit your feedback
content_review(f"{feedback_prefix}_Standard_loss_functions_with_various_priors_Interactive_Demo_and_Discussion")
セクション6.2: ベイズ意思決定
最後に、これまで学んだことをすべて組み合わせます!
新しいアストロキャットの位置の測定値を受け取ったと想像してください。どこにいるか決めて、どれだけ移動させるか指示する必要があります。ただし、衛星と宇宙ネズミの位置も考慮します。衛星方向の誤差は宇宙ネズミ方向より悪いので、セクション3.2の複雑な効用関数を使います。
インタラクティブデモ 6.2: さまざまな事前分布での複雑な猫のコスト
質問:
- 弱い事前分布と尤度の場合、推定にどれだけ効用関数を頼っていますか?
- 分散の小さい良い測定(尤度)があればどれだけ助かりますか?
- 意思決定で最も重要な要素は何ですか?
# @markdown Execute this cell to enable the widget
mu1_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=-0.5,
description="µ_prior", continuous_update=True)
mu2_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.5,
description="µ_likelihood", continuous_update=True)
sigma1_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5,
description="σ_prior", continuous_update=True)
sigma2_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5,
description="σ_likelihood", continuous_update=True)
mu_g_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.0,
description="µ_gain", continuous_update=True)
sigma_g_slider = FloatSlider(min=0.1, max=2.0, step=0.01, value=0.5,
description="σ_gain", continuous_update=True)
dist_label = Label(value="loss distance: µ1_c - µ2_c",
layout=Layout(display="flex", justify_content="center"))
loc_label = Label(value="loss center: (µ1_c + µ2_c) / 2",
layout=Layout(display="flex", justify_content="center"))
mu_dist_slider = FloatSlider(min=0.0, max=8.0, step=0.01, value=4.0,
description="distance", continuous_update=True)
mu_loc_slider = FloatSlider(min=-4.0, max=4.0, step=0.01, value=0.0,
description="center", continuous_update=True)
widget_ui = HBox([VBox([mu1_slider, sigma1_slider, mu2_slider, sigma2_slider]),
VBox([dist_label, mu_dist_slider, loc_label, mu_loc_slider]),
VBox([mu_g_slider, sigma_g_slider])])
widget_out = interactive_output(plot_utility_mixture_dist,
{'mu1': mu1_slider,
'sigma1': sigma1_slider,
'mu2': mu2_slider,
'sigma2': sigma2_slider,
'mu_g': mu_g_slider,
'sigma_g': sigma_g_slider,
'mu_dist': mu_dist_slider,
'mu_loc': mu_loc_slider,
'plot_utility_row': fixed(True)})
display(widget_ui, widget_out)
# @title Submit your feedback
content_review(f"{feedback_prefix}_Complicated_cat_costs_with_various_priors_Interactive_Demo_and_Discussion")
まとめ
チュートリアルの推定所要時間: 1時間35分
このチュートリアルでは、アストロキャットの位置を見つけて選ぶ文脈でベイズの定理とベイズアプローチの探求を拡張しました。
具体的には以下を扱いました:
-
ガウス分布とその性質
-
尤度は測定値が隠れ状態のもとでの確率であること
-
ガウス分布間の情報共有(PDFの掛け算や2次元分布を通じて)
-
事前分布と尤度がどのように相互作用して事後分布(測定後の隠れ状態の確率)を作るかは共分散に依存すること
-
効用は行動と状態の組み合わせから得られる利益であり、行動の期待効用はその状態の確率で重み付けされた効用の総和であり、期待効用が最大の行動を選べること