(Your student-id is the number shown on your student card.)
E.g. if you work with 3 people, the notebook should be named: 12301230_3434343_1238938934_lab1.ipynb.
This will be parsed by a regexp, so please double check your filename.
Before you turn this problem in, please make sure everything runs correctly. First, restart the kernel (in the menubar, select Kernel$\rightarrow$Restart) and then run all cells (in the menubar, select Cell$\rightarrow$Run All).
Make sure you fill in any place that says YOUR CODE HERE or "YOUR ANSWER HERE", as well as your names and email adresses below.
In [1]:
NAME = "Laura Ruis"
NAME2 = "Fredie Haver"
NAME3 = "Lukás Jelínek"
EMAIL = "lauraruis@live.nl"
EMAIL2 = "frediehaver@hotmail.com"
EMAIL3 = "lukas.jelinek1@gmail.com"
Notes on implementation:
In [2]:
%pylab inline
plt.rcParams["figure.figsize"] = [20,10]
Write a method gen_cosine(N) that generates toy data like in fig 1.2 of Bishop's book. The method should have a parameter $N$, and should return $N$-dimensional vectors $\bx$ and $\bt$, where $\bx$ contains evenly spaced values from 0 to (including) 2$\pi$, and the elements $t_i$ of $\bt$ are distributed according to:
where $x_i$ is the $i$-th elements of $\bf{x}$, the mean $\mu_i = \cos(x_i)$ and the standard deviation $\sigma = 0.2$.
In [3]:
def gen_cosine(n):
x = np.array(np.linspace(0, 2 * np.pi, n))
t = np.random.normal(np.cos(x), 0.2)
return x, t
In [4]:
def gen_real_cosine(n):
real_x = np.array(np.linspace(0, 2 * np.pi, n))
real_y = np.cos(real_x)
return real_x,real_y
In [5]:
### Test your function
np.random.seed(5)
N = 10
x, t = gen_cosine(N)
assert x.shape == (N,), "the shape of x is incorrect"
assert t.shape == (N,), "the shape of t is incorrect"
Write a method fit_polynomial(x, t, M) that finds the maximum-likelihood solution of an unregularized $M$-th order polynomial for some dataset x. The error function to minimize w.r.t. $\bw$ is:
$E(\bw) = \frac{1}{2} (\bPhi\bw - \bt)^T(\bPhi\bw - \bt)$
where $\bPhi$ is the feature matrix (or design matrix) as explained in Bishop's book at section 3.1.1, $\bt$ is the vector of target values. Your method should return a vector $\bw$ with the maximum-likelihood parameter estimates, as well as the feature matrix $\bPhi$.
In [6]:
def designmatrix(x, M): # it is highly recommended to write a helper function that computes Phi
design_matrix = []
for i in range(M+1):
design_matrix.append([data ** i for data in x])
design_matrix = np.matrix(design_matrix).transpose()
return design_matrix
In [7]:
def LSE(phi, t):
phi_squared_inv = np.linalg.inv(np.matmul(phi.transpose(), phi))
mp_pseudo_inv = np.matmul(phi_squared_inv, phi.transpose())
return np.reshape(np.array(np.matmul(mp_pseudo_inv, t)), phi.shape[1])
def fit_polynomial(x, t, M):
Phi = designmatrix(x, M)
w_ml = LSE(Phi, t)
return w_ml, Phi
In [8]:
### Test your function
N = 10
x = np.square((np.linspace(-1, 1, N)))
t = 0.5*x + 1.5
m = 2
w, Phi = fit_polynomial(x,t,m)
assert w.shape == (m+1,), "The shape of w is incorrect"
assert Phi.shape == (N, m+1), "The shape of Phi is incorrect"
Sample a dataset with $N=10$, and fit four polynomials with $M \in (0, 2, 4, 8)$.
For each value of $M$, plot the prediction function, along with the data and the original cosine function. The resulting figure should look similar to fig 1.4 of the Bishop's book. Note that you can use matplotlib's plt.pyplot(.) functionality for creating grids of figures.
In [9]:
N = 10
x,t = gen_cosine(N)
real_x , real_y = gen_real_cosine(100)
w_ml = []
Phi = []
M = [0, 2, 4, 8]
for m in M:
w,phi = fit_polynomial(x,t,m)
w_ml.append(w)
Phi.append(phi)
predictions = []
for i in range(len(M)):
predictions.append(np.squeeze(np.asarray(np.matmul(Phi[i],w_ml[i].transpose()))))
f, ((ax1, ax2), (ax4, ax8)) = plt.subplots(2, 2, sharex='col', sharey='row')
ax1.scatter(x,t)
ax1.plot(real_x,real_y)
ax1.plot(x,predictions[0])
ax1.set_title("M = "+ str(M[0]))
ax2.scatter(x,t)
ax2.plot(real_x,real_y)
ax2.plot(x,predictions[1])
ax2.set_title("M = "+ str(M[1]))
ax4.scatter(x,t)
ax4.plot(real_x,real_y)
ax4.plot(x,predictions[2])
ax4.set_title("M = "+ str(M[2]))
ax8.scatter(x,t)
ax8.plot(real_x,real_y)
ax8.plot(x,predictions[3])
ax8.set_title("M = "+ str(M[3]))
Out[9]:
Write a method fit_polynomial_reg(x, t, M, lamb) that fits a regularized $M$-th order polynomial to the periodic data, as discussed in the lectures, where lamb is the regularization term lambda. (Note that 'lambda' cannot be used as a variable name in Python since it has a special meaning). The error function to minimize w.r.t. $\bw$:
$E(\bw) = \frac{1}{2} (\bPhi\bw - \bt)^T(\bPhi\bw - \bt) + \frac{\lambda}{2} \mathbf{w}^T \mathbf{w}$
For background, see section 3.1.4 of Bishop's book.
The function should return $\bw$ and $\bPhi$.
In [10]:
def LSEReg(phi,t,lamb):
lambid = lamb*np.identity(phi.shape[1])
phi_squared = np.matmul(phi.transpose(), phi)
lamphi_inv = np.linalg.inv(numpy.add(lambid,phi_squared))
mp_pseudo_inv = np.matmul(lamphi_inv, phi.transpose())
return np.reshape(np.array(np.matmul(mp_pseudo_inv, t)), phi.shape[1])
def fit_polynomial_reg(x, t, m, lamb):
Phi = designmatrix(x, m)
w = LSEReg(Phi, t, lamb)
return w, Phi
In [11]:
### Test your function
N = 10
x = np.square((np.linspace(-1, 1, N)))
t = 0.5*x + 1.5
m = 2
lamb = 0.1
w, Phi = fit_polynomial_reg(x,t,m, lamb)
assert w.shape == (m+1,), "The shape of w is incorrect"
assert Phi.shape == (N, m+1), "The shape of w is incorrect"
Use cross-validation to find a good choice of $M$ and $\lambda$, given a dataset of $N=10$ datapoints generated with gen_cosine(20). You should write a function that tries (loops over) a reasonable range of choices of $M$ and $\lambda$, and returns the choice with the best cross-validation error. In this case you use $K=5$ folds.
You can let $M \in (0, 1, ..., 10)$, and let $\lambda \in (e^{-10}, e^{-9}, ..., e^{0})$.
a) (5 points) First of all, write a method pred_error(x_train, x_valid, t_train, t_valid, M, lamb) that compares the prediction of your method fit_polynomial_reg for a given set of parameters $M$ and $\lambda$ to t_valid. It should return the prediction error for a single fold.
In [12]:
def pred_error(x_train, x_valid, t_train, t_valid, M, reg):
w, Phi = fit_polynomial_reg(x_train,t_train,M,reg)
# w, Phi = fit_polynomial(x_train,t_train,M)
valid_Phi = designmatrix(x_valid, M)
t_trained = np.matmul(valid_Phi, w.transpose())
err =0.5*np.matmul((t_valid - t_trained),(t_valid - t_trained).transpose())
return err
In [13]:
### Test your function
N = 10
x = np.linspace(-1, 1, N)
t = 0.5*np.square(x) + 1.5
M = 2
reg = 0.1
pred_err = pred_error(x[:-2], x[-2:], t[:-2], t[-2:], M, reg)
assert pred_err < 0.01, "pred_err is too big"
b) (10 points) Now write a method find_best_m_and_lamb(x, t) that finds the best values for $M$ and $\lambda$. The method should return the best $M$ and $\lambda$. To get you started, here is a method you can use to generate indices of cross-validation folds.
In [14]:
def kfold_indices(N, k):
all_indices = np.arange(N,dtype=int)
np.random.shuffle(all_indices)
idx = [int(i) for i in np.floor(np.linspace(0,N,k+1))]
train_folds = []
valid_folds = []
for fold in range(k):
valid_indices = all_indices[idx[fold]:idx[fold+1]]
valid_folds.append(valid_indices)
train_folds.append(np.setdiff1d(all_indices, valid_indices))
return train_folds, valid_folds
In [15]:
def find_best_m_and_lamb(x, t):
folds = 5
train_idx, valid_idx = kfold_indices(len(x),folds)
M_best = 2
lamb_best = 1
err_best = 1000000000
for m in range(11):
for exponent in range(0,11):
total_err = 0
for i in range(folds):
x_train = np.take(x, train_idx[i])
t_train = np.take(t, train_idx[i])
x_valid = np.take(x, valid_idx[i])
t_valid = np.take(t, valid_idx[i])
total_err =total_err + pred_error(x_train, x_valid, t_train, t_valid, m, math.exp(-exponent))
# print(m,exponent,err)
if total_err < err_best:
err_best = total_err
M_best = m
lamb_best = math.exp(-exponent)
return M_best, lamb_best
In [16]:
### If you want you can write your own test here
In [17]:
np.random.seed(9)
N=10
x, t = gen_cosine(N)
print(find_best_m_and_lamb(x,t))
For some dataset with $N = 10$, plot the model with the optimal $M$ and $\lambda$ according to the cross-validation error, using the method you just wrote. In addition, the plot should show the dataset itself and the function that we try to approximate. Let the plot make clear which $M$ and $\lambda$ were found.
In [18]:
N=20
x, t = gen_cosine(N)
real_x,real_y = gen_real_cosine(100)
M,lamb = find_best_m_and_lamb(x,t)
w, Phi = fit_polynomial(x,t,M)
res = np.matmul(Phi,w.transpose())
res = np.asarray(res).reshape(-1)
plt.plot(real_x,real_y)
plt.scatter(x,t)
plt.plot(x,res)
plt.figtext(0.6, 0.8, "M = "+ str(M), size='xx-large')
plt.figtext(0.6, 0.77, "Lambda = exp("+ str(math.log(lamb))+" )", size='xx-large')
plt.show()
In [19]:
def gen_cosine2(n):
x = np.array(np.random.uniform(0, 2 * np.pi, n))
t = np.random.normal(np.cos(x), 0.2)
return x, t
In [20]:
### Test your function
np.random.seed(5)
N = 10
x, t = gen_cosine2(N)
assert x.shape == (N,), "the shape of x is incorrect"
assert t.shape == (N,), "the shape of t is incorrect"
You're going to implement a Bayesian linear regression model, and fit it to the periodic data. Your regression model has a zero-mean isotropic Gaussian prior over the parameters, governed by a single (scalar) precision parameter $\alpha$, i.e.:
$$p(\bw \;|\; \alpha) = \mathcal{N}(\bw \;|\; 0, \alpha^{-1} \bI)$$The covariance and mean of the posterior are given by:
$$\bS_N= \left( \alpha \bI + \beta \bPhi^T \bPhi \right)^{-1} $$$$\bm_N = \beta\; \bS_N \bPhi^T \bt$$where $\alpha$ is the precision of the predictive distribution, and $\beta$ is the noise precision. See MLPR chapter 3.3 for background.
Write a method fit_polynomial_bayes(x, t, M, alpha, beta) that returns the mean $\bm_N$ and covariance $\bS_N$ of the posterior for a $M$-th order polynomial. In addition it should return the design matrix $\bPhi$. The arguments x, t and M have the same meaning as in question 1.2.
In [21]:
def bayesian_LR(phi, t, alpha, beta, M):
phi_squared = np.matmul(phi.transpose(), phi)
Alpha = np.multiply(alpha, np.identity(M + 1))
cov = np.linalg.inv(np.add(Alpha, np.multiply(beta, phi_squared)))
phi_t = np.matmul(phi.transpose(), t).transpose()
mean = np.multiply(beta, np.matmul(cov, phi_t))
mean = np.reshape(np.array(mean), M+1)
return cov, mean
def fit_polynomial_bayes(x, t, M, alpha, beta):
Phi = designmatrix(x, M)
S, m = bayesian_LR(Phi, t, alpha, beta, M)
return m, S, Phi
In [22]:
### Test your function
N = 10
x = np.linspace(-1, 1, N)
t = 0.5*np.square(x) + 1.5
M = 2
alpha = 0.5
beta = 25
m, S, Phi = fit_polynomial_bayes(x, t, M, alpha, beta)
assert m.shape == (M+1,), "the shape of m is incorrect"
assert S.shape == (M+1, M+1), "the shape of S is incorrect"
assert Phi.shape == (N, M+1), "the shape of Phi is incorrect"
The predictive distribution of Bayesian linear regression is:
$$ p(t \;|\; \bx, \bt, \alpha, \beta) = \mathcal{N}(t \;|\; \bm_N^T \phi(\bx), \sigma_N^2(\bx))$$$$ \sigma_N^2 = \frac{1}{\beta} + \phi(\bx)^T \bS_N \phi(\bx) $$where $\phi(\bx)$ are the computed features for a new datapoint $\bx$, and $t$ is the predicted variable for datapoint $\bx$.
Write a function that predict_polynomial_bayes(x, m, S, beta) that returns the predictive mean, variance and design matrix $\bPhi$ given a new datapoint x, posterior mean m, posterior variance S and a choice of model variance beta.
In [23]:
def predict_polynomial_bayes(x, m, S, beta):
M = len(m)
sigma = []
mean = []
for i in x:
phi = designmatrix([i], M - 1)
var = 1/beta + np.matmul(np.reshape(phi.transpose(),(1,M)), np.reshape(np.matmul(S, phi.transpose()),(M,1)))
mean.append(np.matmul(np.reshape(m, (1,M)), phi.transpose()))
sigma.append(var)
mean = np.reshape(np.array(mean), (N,))
sigma = np.reshape(np.array(sigma), (N,))
Phi = designmatrix(x, M - 1)
return mean, sigma, Phi
In [24]:
### Test your function
np.random.seed(5)
N = 10
x = np.linspace(-1, 1, N)
m = np.empty(3)
S = np.empty((3, 3))
beta = 25
mean, sigma, Phi = predict_polynomial_bayes(x, m, S, beta)
assert mean.shape == (N,), "the shape of mean is incorrect"
assert sigma.shape == (N,), "the shape of sigma is incorrect"
assert Phi.shape == (N, m.shape[0]), "the shape of Phi is incorrect"
a) (5 points) Generate 10 datapoints with gen_cosine2(10). Compute the posterior mean and covariance for a Bayesian polynomial regression model with $M=4$, $\alpha=\frac{1}{2}$ and $\beta=\frac{1}{0.2^2}$.
Plot the Bayesian predictive distribution, where you plot (for $x$ between 0 and $2 \pi$) $t$'s predictive mean and a 1-sigma predictive variance using plt.fill_between(..., alpha=0.1) (the alpha argument induces transparency).
Include the datapoints in your plot.
In [25]:
np.random.seed(101)
real_x, real_y = gen_real_cosine(100)
N = 10
M = 4
alpha = 0.5
x, t = gen_cosine2(N)
perm = x.argsort()
beta = 25
m, S, Phi = fit_polynomial_bayes(x, t, M, alpha, beta)
plt.plot(real_x, real_y)
plt.scatter(x[perm], t[perm])
mean_pt, var, Phi = predict_polynomial_bayes(x, m, S, beta)
plt.plot(x[perm], mean_pt[perm])
plt.fill_between(x[perm],mean_pt[perm] - var[perm],mean_pt[perm] + var[perm], alpha=0.1)
plt.show()
b) (5 points) For a second plot, draw 100 samples from the parameters' posterior distribution. Each of these samples is a certain choice of parameters for 4-th order polynomial regression. Display each of these 100 polynomials.
In [26]:
M = 4
alpha = 0.5
beta = 25
real_x, real_y = gen_real_cosine(100)
perm = x.argsort()
plt.plot(real_x, real_y)
for i in range(1, 100):
m, S, Phi = fit_polynomial_bayes(x, t, M, alpha, i)
mean, var, phi = predict_polynomial_bayes(x, m, S, i)
plt.plot(x[perm], mean[perm])
plt.show()
Since we know the distribution that generated the data has a standard deviation of 0.2 (variance of 0.2^2), the best choice for the precision parameter is the inverse variance.
b) (5 points) What problems do we face when it comes to choosing basis functions in linear models?
On the one hand we want to have a good fit to the training data, on the other hand we do not want to overfit the data, which allows for a greater prediction error. The problem we face is dimensionality, how many basis function do we need to best approximate the underlying structure of the data.
In [ ]: