<jupyter_start><jupyter_text>Plot  a Bunch of Samples from the prior distribution<jupyter_code>for i in range(10):
    plt.plot(x_samples,mvn.rvs())<jupyter_output><empty_output><jupyter_text>$f∗|X∗, X,f ∼ \mathcal{N}(K(X∗, X)K(X, X)^{−1}f,K(X∗, X∗) − K(X∗, X)K(X, X)^{−1}K(X, X∗)).$$p(x|m, Σ) = (2π)−D/2|Σ|−1/2exp^{−\frac{1}{2}(x − m)'Σ^{−1}(x − m)}$<jupyter_code>#D = np.array([[3,.1],[4,.1],[5,.1],[40,2],[80,-4]])
xs = np.linspace(0.01,10,num=60)
ys = np.sin(.5*np.pi*xs)*3 + np.sin(2*np.pi*xs) + (1/np.exp(xs/2))*6*np.sin(4*np.pi*xs) + xs//2
#ys = (1/np.exp(xs/2))*4*np.sin(2*np.pi*xs) + xs/2
#ys = np.hstack([np.zeros(50),np.ones(50)*2]) - 1
plt.plot(xs,ys)

D = np.vstack([xs,ys+np.random.normal(scale=1,size=xs.shape)]).T

plt.plot(xs,D[:,1])

n_samps = 220
x_samples = np.linspace(0,10,num=n_samps)
plt.plot(1/np.exp(xs))
#define parameters for the GP regression
k = 2
l = 1e-1
K = np.zeros([D.shape[0],D.shape[0]])


for idx1,samp1 in enumerate(D[:,0]):
    
    for idx2,samp2 in enumerate(D[:,0]):
        
        K[idx1,idx2] = kernel(samp1,samp2)  #+= 1e-6
        #if idx1==idx2:
K += 2*np.eye(D.shape[0])



Kst =  np.zeros([x_samples.shape[0],D.shape[0]])


for idx1,samp1 in enumerate(x_samples):
    
    for idx2,samp2 in enumerate(D[:,0]):
        
        Kst[idx1,idx2] = kernel(samp1,samp2)  #+= 1e-6
        #if idx1==idx2:


Kstst =  np.zeros([x_samples.shape[0],x_samples.shape[0]])


for idx1,samp1 in enumerate(x_samples):
    
    for idx2,samp2 in enumerate(x_samples):
        
        Kstst[idx1,idx2] = kernel(samp1,samp2)  #+= 1e-6
        #if idx1==idx2:



        
GP_mean = Kst.dot(np.linalg.inv(K)).dot(D[:,1])
GP_cov = (Kstst - Kst.dot(np.linalg.inv(K)).dot(Kst.T)) + np.eye(GP_mean.shape[0])*1e-6
GP_posterior = multivariate_normal(GP_mean,GP_cov)<jupyter_output><empty_output><jupyter_text>This is success enjoy the beautiful plot of samples from the posterior<jupyter_code>plt.figure(figsize=(12,6))

ax = plt.gca()

#for i in range(4):
#    ax.plot(x_samples,GP_posterior.rvs())
    
#ax.plot(x_samples,GP_posterior.rvs(),label='Posterior Samples')


ax.plot(x_samples,GP_mean,linewidth=2,color='k',label='Posterior Mean')



ax.fill_between(x_samples, GP_mean-np.sqrt(np.diag(GP_cov)), GP_mean+np.sqrt(np.diag(GP_cov)), facecolor=[.3]*3,alpha=.4, interpolate=True,label='CI')
ax.plot([0], [0], color=[.3]*3,alpha=.4,linewidth=12, label='$\pm$2SD')


ax.plot(D[:,0],ys,color=[.9,.4,.9],linestyle='--',linewidth=4,label='True Function')
ax.scatter(D[:,0],D[:,1],s=64,color=[.8,.3,.2])

ax.legend(loc='best')
plt.xlabel(r'$x$')
plt.ylabel(r'$f(x)$')
plt.xlim(xmin=0,xmax=6.3)<jupyter_output><empty_output><jupyter_text>Samples From Posterior Distribution of Functions<jupyter_code>for idx1,samp1 in enumerate(x_samples):
    
    for idx2,samp2 in enumerate(x_samples):
        
        
        SIGMA[idx1,idx2] = kernel(samp1,samp2)  #+= 1e-6
        #if idx1==idx2:
        #    SIGMA[idx1,idx2]
            
SIGMA += np.eye(SIGMA.shape[0])*1e-6

mu = np.zeros(x_samples.shape)
GP_prior = multivariate_normal(mu,cov=SIGMA)
plt.figure(figsize=(16,8))
plt.suptitle('Samples from Distributions of GP \n\n')

plt.subplot(1,2,1)
plt.title('Samples From Prior')
ax = plt.gca()
for i in range(3):
    ax.plot(x_samples,mvn.rvs())
    
ax.plot(x_samples,mvn.rvs(),label='Samples')
plt.xlabel(r'$x$')
plt.ylabel(r'$f(x)$')


plt.subplot(1,2,2)
ax = plt.gca()
plt.title('\nSamples From Posterior')
for i in range(3):
    ax.plot(x_samples,GP_posterior.rvs())
    
ax.plot(x_samples,GP_posterior.rvs(),label='Samples')
plt.xlabel(r'$x$')<jupyter_output><empty_output><jupyter_text>Repeat Procedure in 2D<jupyter_code>n_samps = 100
x_pos = np.linspace(1,5,num=n_samps)
y_pos = np.linspace(3,5,num=n_samps)

kernel = lambda x1,x2: k*np.exp((-.5)*((x1-x2)**2)/l**2)
SIGMA = np.zeros([n_samps,n_samps])
#define parameters for the GP regression
k = 3
l = .5
for idx1,samp1 in enumerate(x_samples):
    
    for idx2,samp2 in enumerate(x_samples):
        
        
        SIGMA[idx1,idx2] = kernel(samp1,samp2)  #+= 1e-6
        #if idx1==idx2:
        #    SIGMA[idx1,idx2]
            
SIGMA += np.eye(SIGMA.shape[0])*1e-6<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.21327919751612132, 'fraction_numerical': 0.0489610699785049, 'mean_word_length': 4.4319066147859925, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '374153', 'license': 'MIT', 'n_tokens_mistral': 2031, 'n_tokens_neox': 1929, 'n_words': 254, 'path': 'GP_regression_first_test.ipynb', 'repo_name': 'yves-weissenberger/GPR', 'size': '4187'}