<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 $$
\mathcal{L}(r, a) = \log \mathcal{N}(y_r; W_r a, q)  \prod_{i\neq r} \mathcal{N}(y_i; W_i a, s) 
$$
$$
\log\mathcal{N}(y_r; W_r a, q) = -\frac{1}{2}\log 2\pi q - \frac{1}{2} \frac{1}{q} (y_r - W_r a)^2
$$
$$
\log\mathcal{N}(y_i; W_i a, s) = -\frac{1}{2}\log 2\pi s - \frac{1}{2} \frac{1}{s} (y_i - W_i a)^2
$$
$$
\mathcal{L}(r, a) = -\frac{1}{2}\log 2\pi q - \frac{1}{2} \frac{1}{q} (y_r - W_r a)^2 + \sum_{i\neq r} -\frac{1}{2}\log 2\pi s - \frac{1}{2} \frac{1}{s} (y_i - W_i a)^2 
$$
$$
\mathcal{L}(r, a) = -\frac{1}{2}\log 2\pi q - \sum_{i\neq r} \frac{1}{2}\log 2\pi s - \frac{1}{2} \frac{1}{q} (y_r - W_r a)^2 - \frac{1}{2} \frac{1}{s} \sum_{i\neq r}  (y_i - W_i a)^2 
$$
$$
D_r = \diag(s,s,\dots, q, \dots, s)^{-1}
$$
$$
\mathcal{L}(r, a) = -\frac{1}{2}\log 2\pi q - \frac{N-1}{2}\log 2\pi s - \frac{1}{2} (y - Wa)^\top D_r (y - Wa) 
$$
$$
\mathcal{L}(r, a) =^+ - \frac{1}{2} (y - Wa)^\top D_r (y - Wa) 
$$
\begin{eqnarray}
\mathcal{L}(r, a) =^+ - \frac{1}{2} y D_r y^\top  + y^\top D_r W a - \frac{1}{2} a^\top W^\top D_r Wa
\end{eqnarray}
\begin{eqnarray}
\frac{\partial}{\partial a}\mathcal{L}(r, a) & = &    W^\top D_r y  - W^\top D_r W a = 0 \
 W^\top D_r y  & = & W^\top D_r W a = 0 \
 (W^\top D_r W)^{-1} W^\top D_r y  & = &  a_r^*  \
\end{eqnarray}
To use standard Least Square solver, we substitute
\begin{eqnarray}
V_r^\top \equiv W^\top D_r^{1/2} \
V_r \equiv  D_r^{1/2} W \
\end{eqnarray}
\begin{eqnarray}
(V_r^\top V_r)^{-1} V_r^\top D_r^{1/2} y  & = &  a_r^* 
\end{eqnarray}
In Matlab/Octave this is least square with
\begin{eqnarray}
 a_r^*  = V_r \backslash D_r^{1/2} y
\end{eqnarray}
Step1:  Todo
Step2:  SUppose we are given a data set $(y_i, x_i)$ for $i=1\dots N$
Assume we have a basis regression model (for example a polynomial basis where $f_k(x) = x^k$) and wish to fit
$y_i = \sum_k A_{ik} w_k + \epsilon_i$
for all $i = 1 \dots N$ where 
$
A_{ik} = f_k(x_i)
$
Assume the prior
$
w \sim \mathcal{N}(w; 0, P)
$
Derive an expression for $p(y_{\text{new}}| x_{\text{new}}, y_{1
<ASSISTANT_TASK:>
Python Code:
import scipy.linalg as la
LL = np.zeros(N)
for rr in range(N):
    ss = s*np.ones(N)
    ss[rr] = q
    D_r = np.diag(1/ss)
    V_r = np.dot(np.sqrt(D_r), W)
    b = y/np.sqrt(ss)
    a_r,re,ra, cond = la.lstsq(V_r, b)
    e = (y-np.dot(W, a_r))/np.sqrt(ss)
    LL[rr] = -0.5*np.dot(e.T, e)
    print(LL[rr])
#plt.plot(x, y, 'o')
#plt.plot(x, np.dot(W, a_r),'-')
#plt.plot(e)
plt.plot(LL)
plt.show()
Explanation: $$
\mathcal{L}(r, a) = \log \mathcal{N}(y_r; W_r a, q)  \prod_{i\neq r} \mathcal{N}(y_i; W_i a, s) 
$$
$$
\log\mathcal{N}(y_r; W_r a, q) = -\frac{1}{2}\log 2\pi q - \frac{1}{2} \frac{1}{q} (y_r - W_r a)^2
$$
$$
\log\mathcal{N}(y_i; W_i a, s) = -\frac{1}{2}\log 2\pi s - \frac{1}{2} \frac{1}{s} (y_i - W_i a)^2
$$
$$
\mathcal{L}(r, a) = -\frac{1}{2}\log 2\pi q - \frac{1}{2} \frac{1}{q} (y_r - W_r a)^2 + \sum_{i\neq r} -\frac{1}{2}\log 2\pi s - \frac{1}{2} \frac{1}{s} (y_i - W_i a)^2 
$$
$$
\mathcal{L}(r, a) = -\frac{1}{2}\log 2\pi q - \sum_{i\neq r} \frac{1}{2}\log 2\pi s - \frac{1}{2} \frac{1}{q} (y_r - W_r a)^2 - \frac{1}{2} \frac{1}{s} \sum_{i\neq r}  (y_i - W_i a)^2 
$$
$$
D_r = \diag(s,s,\dots, q, \dots, s)^{-1}
$$
$$
\mathcal{L}(r, a) = -\frac{1}{2}\log 2\pi q - \frac{N-1}{2}\log 2\pi s - \frac{1}{2} (y - Wa)^\top D_r (y - Wa) 
$$
$$
\mathcal{L}(r, a) =^+ - \frac{1}{2} (y - Wa)^\top D_r (y - Wa) 
$$
\begin{eqnarray}
\mathcal{L}(r, a) =^+ - \frac{1}{2} y D_r y^\top  + y^\top D_r W a - \frac{1}{2} a^\top W^\top D_r Wa
\end{eqnarray}
\begin{eqnarray}
\frac{\partial}{\partial a}\mathcal{L}(r, a) & = &    W^\top D_r y  - W^\top D_r W a = 0 \
 W^\top D_r y  & = & W^\top D_r W a = 0 \
 (W^\top D_r W)^{-1} W^\top D_r y  & = &  a_r^*  \
\end{eqnarray}
To use standard Least Square solver, we substitute
\begin{eqnarray}
V_r^\top \equiv W^\top D_r^{1/2} \
V_r \equiv  D_r^{1/2} W \
\end{eqnarray}
\begin{eqnarray}
(V_r^\top V_r)^{-1} V_r^\top D_r^{1/2} y  & = &  a_r^* 
\end{eqnarray}
In Matlab/Octave this is least square with
\begin{eqnarray}
 a_r^*  = V_r \backslash D_r^{1/2} y
\end{eqnarray}
End of explanation
import numpy as np
import scipy as sc
import scipy.linalg as la
def cond_Gauss(Sigma, mu, idx1, idx2, x2):
    Sigma11 = Sigma[idx1, idx1].reshape((len(idx1),len(idx1)))
    Sigma12 = Sigma[idx1, idx2].reshape((len(idx1),len(idx2)))
    Sigma22 = Sigma[idx2, idx2].reshape((len(idx2),len(idx2)))
#    print(Sigma11)
#    print(Sigma12)
#    print(Sigma22)
    
    
    mu1 = mu[idx1]
    mu2 = mu[idx2]
    G = np.dot(Sigma12, la.inv(Sigma22))
    cond_Sig_1 =  Sigma11 - np.dot(G, Sigma12.T)
    cond_mu_1 = mu1 + np.dot(G, (x2-mu2))
    return cond_mu_1, cond_Sig_1
mu = np.array([0,0])
#P = np.array([2])
#A = np.array([1])
idx1 = [0]
idx2 = [1]
x2 = 5
P = np.array(3).reshape((len(idx1), len(idx1)))
A = np.array(-1).reshape((len(idx2), len(idx1)))
rho = np.array(0)
#Sigma = np.array([[P, A*P],[P*A, A*P*A + rho ]])
I = np.eye(len(idx2))
Sigma = np.concatenate((np.concatenate((P,np.dot(P, A.T)),axis=1), np.concatenate((np.dot(A, P),np.dot(np.dot(A, P), A.T ) + rho*I ),axis=1)))
print(Sigma)
#print(mu)
cond_mu_1, cond_Sig_1 = cond_Gauss(Sigma, mu, idx1, idx2, x2)
print('E[x_1|x_2 = {}] = '.format(x2) , cond_mu_1)
print(cond_Sig_1)
Explanation: Todo: Evaluate the likelihood for all polynomial orders $K=1 \dots 8$
$p(x_1, x_2) = \mathcal{N}(\mu, \Sigma)$
$\mu = \left(\begin{array}{c} \mu_{1} \
 \mu_{2} \end{array} \right)$
$\Sigma = \left(\begin{array}{cc} \Sigma_{11} & \Sigma_{12} \
\Sigma_{12}^\top & \Sigma_{22} \end{array} \right)$
$
p(x_1 | x_2) = \mathcal{N}(\mu_1 + \Sigma_{12} \Sigma_{22}^{-1} (x_2 -\mu_2), \Sigma_{11} - \Sigma_{12} \Sigma_{22}^{-1}\Sigma_{12}^\top)
$
End of explanation
# Use this code to generate a dataset
N = 30
K = 4
s = 0.1
q = 10*s
x = 2*np.random.randn(N)
e = np.sqrt(s) * np.random.randn(N)
# Create the vandermonde matrix
A = x.reshape((N,1))**np.arange(K).reshape(1,K)
w = np.array([0,-1,0.5,0])
y = np.dot(A, w) + e
plt.plot(x, y, 'o')
#plt.plot(e)
plt.show()
# Sig = [P, A.T; A A*A.T+rho*I]
N1 = 3
N2 = 7
P = np.random.randn(N1,N1)
A = np.random.randn(N2,N1)
#Sig11 = np.mat(P)
#Sig12 = np.mat(A.T)
#Sig21 = np.mat(A)
#Sig22 = Sig21*Sig12
Sig11 = np.mat(P)
Sig12 = np.mat(A.T)
Sig21 = np.mat(A)
Sig22 = Sig21*Sig12
print(Sig11.shape)
print(Sig12.shape)
print(Sig21.shape)
print(Sig22.shape)
W = np.bmat([[Sig11, Sig12],[Sig21, Sig22]])
Sig22.shape
3500*1.18*12
%matplotlib inline
import matplotlib.pylab as plt
import numpy as np
x = np.array([3.7, 2.3, 6.9, 7.5])
N = len(x)
lam = np.arange(0.05,10,0.01)
ll = -N*np.log(lam) - np.sum(x)/lam
plt.plot(lam, np.exp(ll))
plt.plot(np.mean(x), 0, 'ok')
plt.show()
xx = np.arange(0, 10, 0.01)
lam = 1000
p = 1/lam*np.exp(-xx/lam)
plt.plot(xx, p)
plt.plot(x, np.zeros((N)), 'ok')
plt.ylim((0,1))
plt.show()
1-(5./6.)**4
1-18/37
import numpy as np
N = 7
A = np.diag(np.ones(7))
ep = 0.5
a = 1
idx = [1, 2, 3, 4, 5, 6, 0]
A = ep*A + (1-ep)*A[:,idx]
C = np.array([[a, 1-a, 1-a, a, a, 1-a, 1-a],[1-a, a, a, 1-a, 1-a, a, a]])
p = np.ones((1,N))/N
print(A)
y = [1, 1, 0, 0, 0]
print(p)
p = C[y[0] , :]*p
print(p/np.sum(p, axis=1))
Explanation: SUppose we are given a data set $(y_i, x_i)$ for $i=1\dots N$
Assume we have a basis regression model (for example a polynomial basis where $f_k(x) = x^k$) and wish to fit
$y_i = \sum_k A_{ik} w_k + \epsilon_i$
for all $i = 1 \dots N$ where 
$
A_{ik} = f_k(x_i)
$
Assume the prior
$
w \sim \mathcal{N}(w; 0, P)
$
Derive an expression for $p(y_{\text{new}}| x_{\text{new}}, y_{1:N}, x_{1:N})$ and implement a program that plots the mean and corresponding errorbars (from standard deviation of $p(y_{\text{new}}| x_{\text{new}}, y_{1:N}, x_{1:N})$) by choosing $y_{\text{new}}$ on a regular grid. 
Note that $y_{\text{new}} = \sum f_k(x_{\text{new}}) w_k + \epsilon$
End of explanation
<END_TASK>