import numpy as np

n_sample = 1000
true_coef = np.array([2.0, 3.0])
u0 = np.array([0.0, 0.0])
A0 = 0.1 * np.eye(2)
sample_X = np.random.normal(size=(n_sample, 2))
sample_y = np.dot(sample_X, true_coef.reshape(-1, 1)).reshape(-1,) + np.random.normal(size=n_sample)

tmp1 = np.linalg.inv(np.dot(sample_X.T, sample_X))
tmp2 = np.dot(sample_X.T, sample_y.reshape(-1, 1))
beta = np.dot(tmp1, tmp2)
tmp1 = np.linalg.inv(np.dot(sample_X.T, sample_X) + A0)
tmp2 = np.dot(np.dot(sample_X.T, sample_X), beta) + np.dot(A0, u0.reshape(-1, 1))
estimation_coef = np.dot(tmp1, tmp2)
print(estimation_coef)