def sample_beta(y, X, a0, b0, w):
    n = X.shape[1]
    a1 = a0 + n/2.0
    sum1 = 0
    for i, x in enumerate(X):
        sum1 += (y[i] + np.dot(w,x))**2
    b1 = b0 + sum1
    return gamma.rvs(a=a1, scale=1/b1)