前回のチュートリアルでは、データがノイズを含む線形関係から生成されていると仮定し、平均二乗誤差を最小化することでモデルパラメータを推定する効果的な方法を見つけました。

その際、ノイズは単なる邪魔者として扱いましたが、もしノイズをモデルに直接組み込んだらどうなるでしょうか？

線形モデルを思い出してください:

\begin{align}
y = $\theta x + \epsilon$.
\end{align}

ノイズ成分 $\epsilon$ はしばしばガウス分布（正規分布）からのランダム変数としてモデル化されます。

ガウス分布はその確率密度関数（pdf）で表されます。
\begin{align}
\mathcal{N}(x; \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{1}{2\sigma^2}(x-\mu)^2} \end{align}

この分布は平均 $\mu$ と分散 $\sigma^2$ の2つのパラメータに依存します。ノイズ信号は平均0、分散1のガウス「ホワイトノイズ」とみなすことが多いです：

\begin{align}
$\epsilon \sim \mathcal{N}(0, 1)$.
\end{align}

インタラクティブデモ 1.1: ガウス分布エクスプローラー

以下のエクスプローラーウィジェットを使って、$\mu$ と $\sigma$ のパラメータを変化させるとサンプルの位置や形状がどう変わるかを見てみましょう。

1. $\mu$ を変えるとpdfにどんな影響がありますか？
2. $\sigma$ を変えるとpdfにどんな影響がありますか？

# @markdown Make sure you execute this cell to enable the widget!

@widgets.interact(mu=widgets.FloatSlider(0.0, min=-2.0, max=2.0),
                  sigma=widgets.FloatSlider(1.0, min=0.5, max=2.0))
def plot_normal_dist(mu=0, sigma=1):

  # Generate pdf & samples from normal distribution with mu/sigma
  rv = stats.norm(mu, sigma)
  x = np.linspace(-5, 5, 100)
  y = rv.pdf(x)
  samples = rv.rvs(1000)

  # Plot
  fig, ax = plt.subplots()
  ax.hist(samples, 20, density=True, color='g', histtype='stepfilled',
          alpha=0.8, label='histogram')
  ax.plot(x, y, color='orange', linewidth=3, label='pdf')
  ax.vlines(mu, 0, rv.pdf(mu), color='y', linewidth=3, label='$\mu$')
  ax.vlines([mu-sigma, mu+sigma], 0, rv.pdf([mu-sigma, mu+sigma]), colors='red',
            color='b', linewidth=3, label='$\sigma$')
  ax.set(xlabel='x', ylabel='probability density', xlim=[-5, 5], ylim=[0, 1.0])
  ax.legend()
  plt.show()

解答を見る$

# @title Submit your feedback
content_review(f"{feedback_prefix}_Gaussian_Distribution_Explorer_Interactive_Demo_and_Discussion")

セクション1.2: 確率モデル

チュートリアル開始からここまでの推定所要時間: 11分

ノイズ成分 $\epsilon$ をランダム変数としてモデル化したので、これを元の線形モデルにどう組み込むか考えましょう。ノイズは平均0、分散1のガウス分布に従うとします：$\epsilon \sim \mathcal{N}(0, 1)$。すると、$y$ も平均 $\mu = \theta x$、分散 $\sigma^2 = 1$ のガウス分布からのランダム変数として扱えます：

これは、パラメータ $\theta$ と入力 $x$ が与えられたときに、$y$ が観測される確率が

であることを意味します。

この節および次の節では、単一のデータ点（$x$ と $y$ のペア）に注目するため、添字 $i$ は省略しています（つまり $x_i$ ではなく $x$ を使います）。

さて、真の基礎モデルが $\theta = 1.2$ である元のサンプルデータセットを再度見てみましょう。

# @markdown Execute this cell to generate some simulated data

# setting a fixed seed to our random number generator ensures we will always
# get the same psuedorandom number sequence

np.random.seed(121)
theta = 1.2
n_samples = 30
x = 10 * np.random.rand(n_samples) # sample from a uniform distribution over [0,10)
noise = np.random.randn(n_samples) # sample from a standard normal distribution
y = theta * x + noise

今回は、$p(y|x,\theta=1.2)$ の密度をプロットし、異なる $x$ の値で $p(y)$ がどう変わるかを見てみます。

# @markdown Execute this cell to visualize p(y|x, theta=1.2)

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10, 5))

# Invokes helper function to generate density image plots from data and parameters
im = plot_density_image(x, y, 1.2, ax=ax1)
plt.colorbar(im, ax=ax1)
ax1.axvline(8, color='k')
ax1.set(title=r'p(y | x, $\theta$=1.2)')

# Plot pdf for given x
ylim = ax1.get_ylim()
yy = np.linspace(ylim[0], ylim[1], 50)
ax2.plot(yy, stats.norm(theta * 8, 1).pdf(yy), color='orange', linewidth=2)
ax2.set(
    title=r'p(y|x=8, $\theta$=1.2)',
    xlabel='y',
    ylabel='probability density')
plt.show()

セクション1.3: 尤度推定

チュートリアル開始からここまでの推定所要時間: 15分

確率モデルができたので、次はデータに合う良い $\theta$ の推定値を見つける課題に戻ります。確率の不確実性を考慮して、推定値 $\hat{\theta}$ がデータに適合する尤度を考えます。尤度関数 $\mathcal{L}(\theta)$ は、その $\theta$ でパラメータ化された確率密度関数に等しいです：

コーディング演習 1.3: 尤度関数

この演習では、$\sigma = 1$ の線形モデルに対する尤度関数 $\mathcal{L}(\theta|x, y)$ を実装します。

この関数を実装したら、推定値 $\hat{\theta}$ が与えられた観測値を生成した確率を計算できます。データセットのサンプルの一つで試してみましょう。

ヒント: 指数関数と平方根にはそれぞれ np.exp と np.sqrt を使いましょう。

def likelihood(theta_hat, x, y):
  """The likelihood function for a linear model with noise sampled from a
    Gaussian distribution with zero mean and unit variance.

  Args:
    theta_hat (float): An estimate of the slope parameter.
    x (ndarray): An array of shape (samples,) that contains the input values.
    y (ndarray): An array of shape (samples,) that contains the corresponding
      measurement values to the inputs.

  Returns:
    ndarray: the likelihood values for the theta_hat estimate
  """
  sigma = 1
  ##############################################################################
  ## TODO for students: implement the likelihood function
  # Fill out function and remove
  raise NotImplementedError("Student exercise: implement the likelihood function")
  ##############################################################################

  # Compute Gaussian likelihood
  pdf = ...

  return pdf


print(likelihood(1.0, x[1], y[1]))

サンプル出力

0.11344443599846923

解答を見る$

# @title Submit your feedback
content_review(f"{feedback_prefix}_Likelihood_function_Exercise")

$\mathcal{L}(\theta=1.0|x=2.1,y=3.7) \approx 0.11$ となるはずです。ここまでは良いですが、これが他の推定値より良いことをどう示せるでしょうか？

データ点の集合を扱う場合、我々はそれらの結合確率、すなわち全てのデータ点がパラメータで説明される尤度を考えます。ノイズが各出力に独立に影響すると仮定したので、尤度は分解可能で、次のように書けます：

ここで、$N$ 個のデータ点 $\mathbf{x} = [x_1,...,x_N]$ と $\mathbf{y} = [y_1,...,y_N]$ があります。

実際には、この積は数値的に不安定になることがあります。小さな値を掛け合わせるとアンダーフローが起きるためです。この問題は尤度の対数を取ることで回避できます。対数は積を和に変換するためです：

likelihood メソッドの出力の対数をデータセット全体に適用して和を取ることで、異なる $\hat{\theta}$ の比較がしやすくなります。また、異なる分布密度をデータセット上でプロットして、定性的に比較することもできます。

# @markdown Execute this cell to visualize different distribution densities
theta_hats = [0.5, 1.0, 2.2]
fig, axes = plt.subplots(ncols=3, figsize=(15, 5))
for theta_hat, ax in zip(theta_hats, axes):
  ll = np.sum(np.log(likelihood(theta_hat, x, y)))  # log likelihood
  im = plot_density_image(x, y, theta_hat, ax=ax)
  ax.scatter(x, y)
  ax.set(title=fr'$\hat{{\theta}}$ = {theta_hat}, log likelihood: {ll:.2f}')
plt.colorbar(im, ax=ax)
plt.show()

対数尤度の計算を使うと、$\mathcal{L}(\theta=1.0) > \mathcal{L}(\theta=0.5) > \mathcal{L}(\theta=2.2)$ となります。

これは良いニュースです。これで尤度に基づいて推定器を比較する方法ができました。しかし、MSEアプローチと同様に、最良の推定器を解析的に求めたいです。この場合は尤度を最大化する推定器を見つけます。

セクション1.4: 最尤推定量の求め方

チュートリアル開始からここまでの推定所要時間: 23分

# Compute theta_hat_MLE
theta_hat_mle = (x @ y) / (x @ x)

# @markdown Execute this cell to visualize density with theta_hat_mle

# Plot the resulting distribution density
fig, ax = plt.subplots()
ll = np.sum(np.log(likelihood(theta_hat_mle, x, y))) # log likelihood
im = plot_density_image(x, y, theta_hat_mle, ax=ax)
plt.colorbar(im, ax=ax)
ax.scatter(x, y)
ax.set(title=fr'$\hat{{\theta}}$ = {theta_hat_mle:.2f}, log likelihood: {ll:.2f}')
plt.show()