<jupyter_start><jupyter_text><jupyter_code>!pip install pystan
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pystan
import seaborn
from scipy.stats import mstats
%matplotlib inline<jupyter_output><empty_output><jupyter_text>ガウシアンプロセスサンプラー下記の式に従ってサンプリングを行うガウシアンプロセスサンプラー。 共分散関数$${K(x|α,ρ)_{i,j} =α^2\exp\Big(−\frac{1}{2ρ^2} \sum_{d=1}^D(x_{i,d} −x_{j,d})^2\Big)}$$　 において、$a = 1$, $\rho=1$ を適用し、　 $${K(x|1,1)_{i,j} =\exp\Big(−\frac{1}{2} \sum_{d=1}^D(x_{i,d} −x_{j,d})^2\Big)}$$　 とした共分散関数を利用します。ここで $i, j$ はそれぞれデータ数方向のindexであり、$d$ はベクトル型データの次元方向のindexです。 平均関数平均関数は $ 0$ ベクトルを利用します。 ガウシアンプロセスサンプラー下記のガウシアンプロセスサンプラーを用います。 　$$\begin{align}f(x) & =\mathcal {GP} ({0}, K(x | \alpha, \rho)) \\y &\sim \mathcal N(f(x), 0.3)\end{align}$$　ちなみに実際にサンプリングするときは有限の数のデータ $x$ を利用するので $\mathcal f(x) = {GP}$ を二番目の式に代入してしまって、確率分布の計算を考えた場合には　 $$y \sim \mathcal N({0}, K(x | \alpha, \rho) + 0.3I)$$　 となることろ利用します。$I$ は単位行列で、元々の二番目の式の分だけ分散が大きくなる（共分散は変わらない）という雰囲気だけ…。<jupyter_code>stan_model = """
data {
  int<lower=1> N;
  real x[N];
}

transformed data {
  matrix[N, N] K = cov_exp_quad(x, 1.0, 1.0);
  vector[N] mu = rep_vector(0, N);
  for (n in 1:N) 
    K[n, n] = K[n, n] + 0.3;
}

parameters {
  vector[N] y;
}

model {
  y ~ multi_normal(mu, K);
}
"""
sm = pystan.StanModel(model_code=stan_model)
x_data = np.linspace(-5, 5, 50)
N_data = x_data.shape[0]
y_sample = sm.sampling(iter=500, warmup=100, chains=3, seed=1992, 
                       data={'N': N_data, 'x': x_data})
y_data = y_sample.extract("y")["y"]<jupyter_output><empty_output><jupyter_text>サンプルを可視化合計で1200個のサンプルが得られました。適当に何個かを選んで見てみましょう。<jupyter_code>y_data.shape

plt.figure(figsize=(16, 8))
for i in range(4):
    plt.subplot(2, 2, i+1)
    plt.scatter(x_data, y_data[i, :])<jupyter_output><empty_output><jupyter_text>ガウシアンプロセス回帰上記のガウシアンプロセスサンプラーでは勝手に共分散関数 $K(x|\alpha, \rho)$ を勝手に決めて、サンプリングを行いました。実際のデータの生成プロセスは、まさかこんな再生核ヒルベルト空間を介して出てきているわけがありませんので、これは単にガウシアンプロセスを見るための例題です。もう一度式を見ましょう。　  $$\begin{align}f(x) & =\mathcal {GP} ({0}, K(x | \alpha, \rho)) \\y &\sim \mathcal N(f(x), 1)\end{align}$$ 　      重要なのは、何らかの生成プロセスを通して観測したデータ $(x, y)^N$ に対して、ガウシアンプロセスで回帰を行う場合は、共分散関数 $\alpha, \rho$ を上手く調整ししてやることになります（なぜならガウシアンプロセスサンプラーが発生するデータは、これらのパラメータに依存しているため）。  より一般的には、観測時のノイズ もパラメータ化して  　   $$\begin{align}f(x) & =\mathcal {GP} ({0}, K(x | \alpha, \rho)) \\y &\sim \mathcal N(f(x), \sigma^2)\end{align}$$　  とします。更に事前分布としてパラメータには下記のような分布を敷いておきましょう。特に何も仮定したくなければ無情報事前分布を敷いておけばいいです。  　  $$ \begin{align} \rho &\sim \mathcal {InvGamma} (1, 1) \\ \alpha & \sim \mathcal N(0, 1) \\ \rho & \sim \mathcal N(0, 1) \\ \end{align} $$<jupyter_code>stan_model_fit = """
data {
  int<lower=1> N;
  real x[N];
  vector[N] y;
}

transformed data {
  vector[N] mu = rep_vector(0, N);
}

parameters {
  real<lower=0> rho;
  real<lower=0> alpha;
  real<lower=0> sigma;
}

model {
  matrix[N, N] L_K;
  matrix[N, N] K = cov_exp_quad(x, alpha, rho);
  real sq_sigma = square(sigma);

  // diagonal elements
  for (n in 1:N)
    K[n, n] = K[n, n] + sq_sigma;
  
  L_K = cholesky_decompose(K);
  
  rho ~ inv_gamma(5, 5);
  alpha ~ normal(0, 1);
  sigma ~ normal(0, 1);

  y ~ multi_normal_cholesky(mu, L_K);
}
"""<jupyter_output><empty_output><jupyter_text>補足事項気づいた方もいるかもしれませんが、`cholesky_decompose`という文言が入っています。これは、共分散関数（有限データに対して分散共分散行列）は対称行列ですので、実際 $N \times N$ の行列を保持する意味はありません。兎にも角にも、メモリや計算量を節約しているという認識でいいでしょう。`multi_normal_cholesky`の部分ではコレスキー分解された共分散関数を渡すという約束です（こいつも計算量を節約したガウス分布と思って構わない）。  パラメータの事後分布表示<jupyter_code>sm_fit = pystan.StanModel(model_code=stan_model_fit)
fit = sm_fit.sampling(iter=2000, warmup=500, chains=3, seed=1992, 
                           data={'N': N_data, 'x': x_data, 'y': y_data[0, :] })
rho = fit.extract('rho')['rho']
alpha = fit.extract('alpha')['alpha']
sigma = fit.extract('sigma')['sigma']<jupyter_output><empty_output><jupyter_text>$\rho$ の事後分布データを生成したガウシアンプロセスサンプラーは $1.0$ を設定していたのでした。<jupyter_code>plt.hist(rho, bins=20)<jupyter_output><empty_output><jupyter_text>$\alpha$ の事後分布データを生成したガウシアンプロセスサンプラーは $1.0$ を設定していたのでした。<jupyter_code>plt.hist(alpha, bins=20)<jupyter_output><empty_output><jupyter_text>$\sigma$ の事後分布データを生成したガウシアンプロセスサンプラーは $0.3$ を設定していたのでした。少し高めに出ているように見えますね。<jupyter_code>plt.hist(sigma, bins=20)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1791009649329254, 'fraction_numerical': 0.024476347375853143, 'mean_word_length': 5.921824104234528, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '569947', 'license': 'MIT', 'n_tokens_mistral': 2631, 'n_tokens_neox': 2509, 'n_words': 289, 'path': 'gaussian_process_fit.ipynb', 'repo_name': 'hellocybernetics/practice_pystan', 'size': '4249'}