<jupyter_start><jupyter_text>Sampling from prior<jupyter_code>n_samples = 100
n_funcs = 8

X = np.expand_dims(np.linspace(-4,4, n_samples), 1)
Σ = exp_quadratic(X, X)

ys = np.random.multivariate_normal(mean=np.zeros(n_samples), cov=Σ, size=n_funcs)
for i in range(n_funcs):
    plt.plot(X, ys[i], linestyle='-', marker='o', markersize=3)
plt.xlabel('$x$', fontsize=13)
plt.ylabel('$y = f(x)$', fontsize=13)
plt.title((
    f'{n_funcs} different function realizations at {n_samples} points\n'
    'sampled from a Gaussian process with exponentiated quadratic kernel'))
plt.xlim([-4, 4])
plt.show()
exponentiated_quadratic = exp_quadratic
A = np.array([[1,-2j],[2j,5]])
A, A.shape
L = np.linalg.cholesky(A)
np.dot(L, L.T.conj())
A = [[1,-2j],[2j,5]] # what happens if A is only array_like?
np.linalg.cholesky(A) # an ndarray object is returned
np.linalg.cholesky(np.matrix(A))
def GP(X1, y1, X2, kernel_func):
    cov11 = kernel_func(X1, X1)
    cov12 = kernel_func(X1, X2)
    solved = scipy.linalg.solve(cov11, cov12, assume_a='pos').T
    
    mu2 = solved @ y1
    cov22 = kernel_func(X2, X2)
    cov2 = cov22 - (solved @ cov12)
    return mu2, cov2<jupyter_output><empty_output><jupyter_text>TODO: Convert K-1 by using choleskyhttps://jaketae.github.io/study/gaussian-process/<jupyter_code>def GP2(X1, y1, X2, kernel_func):
    K11 = kernel_func(X1, X1)
    K12 = kernel_func(X1, X2)
    K22 = kernel_func(X2, X2)
    
    #L = np.linalg.cholesky(K11)    
    mu2 = K12.T.dot(np.linalg.inv(K11)).dot(y1)
    cov2 = K22 - K12.T.dot(np.linalg.inv(K11)).dot(K12)
    return mu2, cov2
ny = 10 # Number of functions
domain = (-6, 6)
domain[0]+2, domain[1]-2, (n1, 1)
X1 = np.random.uniform(domain[0]+2, domain[1]-2, size=(n1, 1))
X1.shape
%%prun
f_sin = lambda x: (np.sin(x)).flatten()

n1 = 40 # Train points
n2 = 75 # Test points
ny = 5 # Number of functions
domain = (-6, 6)

X1 = np.random.uniform(domain[0]+2, domain[1]-2, size=(n1, 1))
y1 = f_sin(X1)

X2 = np.linspace(domain[0], domain[1], n2).reshape(-1, 1)

# mu2, cov2 = GP(X1, y1, X2, exp_quadratic)
mu2, cov2 = GP2(X1, y1, X2, exp_quadratic)
sigma2 = np.sqrt(np.diag(cov2))

y2 = np.random.multivariate_normal(mean=mu2, cov=cov2, size=ny)
# Plot the postior distribution and some samples
fig, (ax1, ax2) = plt.subplots(
    nrows=2, ncols=1, figsize=(6, 6))

# Plot the distribution of the function (mean, covariance)
ax1.plot(X2, f_sin(X2), 'b--', label='$sin(x)$')
ax1.fill_between(X2.flat, mu2-2*sigma2, mu2+2*sigma2, color='red', 
                 alpha=0.15, label='$2 \sigma_{2|1}$')
ax1.plot(X2, mu2, 'r-', lw=2, label='$\mu_{2|1}$')
ax1.plot(X1, y1, 'ko', linewidth=2, label='$(x_1, y_1)$')

# Plot some samples from this function
ax2.plot(X2, y2.T, '-')
ax2.set_xlabel('$x$', fontsize=13)
ax2.set_ylabel('$y$', fontsize=13)
ax2.set_title('5 different function realizations from posterior')
ax1.axis([domain[0], domain[1], -3, 3])
ax2.set_xlim([-6, 6])
plt.tight_layout()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18651459241865145, 'fraction_numerical': 0.06642066420664207, 'mean_word_length': 4.5427509293680295, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 1, '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': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '37920', 'license': 'MIT', 'n_tokens_mistral': 1397, 'n_tokens_neox': 1336, 'n_words': 273, 'path': 'experiments/10_gp_from_scratch.ipynb', 'repo_name': 'danhphan/gp_experiments', 'size': '2981'}