def sample_w(y, X, dia_alphas, beta):
    d = X.shape[1]
    S = np.linalg.inv(dia_alphas + beta*X.T*X)
    print(S)
    print(S*(beta*X.T*y))
    m = np.array(S*(beta*X.T*y)).reshape(d,d)
    print(m.shape)
    return multivariate_normal.rvs(mean=m, cov=S)