linreg_elbo

def linreg_elbo(X, t, max_iter=100, tol=1e-4, display=False):
    delta = 1e-8
    n, d = X.shape
    alpha = beta = 1.
    xi = np.ones(n)
    mu = np.ones(d)
    sigma = np.identity(d)
    L = [compute_elbo(X, t, alpha, beta, xi, mu, sigma)]
    if display:
        print('%9s %15s' % ('iteration', 'F'))
        print('%9d %15.11f' % (0, L[-1]))

    for iteration in xrange(1, max_iter + 1):
        xi_old = xi
        xi = np.array([sigma.dot(x.reshape((-1, 1))).dot(x.reshape(1, -1)).trace() + mu.dot(x) ** 2 - 2 * mu.dot(x) * t_n + t_n ** 2
              for x, t_n in zip(X, t)])
        xi = np.sqrt(xi)
        xi[xi < delta] = delta

        alpha = 2 * n / (np.sum(xi_old) + np.sum(xi ** 2 / xi_old))
        beta = d / (sigma.trace() + mu.dot(mu))

        sigma = inv(alpha * ((X).T / xi).dot(X) + beta * np.identity(d))
        mu = alpha * sigma.dot(np.sum(X.T * t / xi, axis=1))
        L.append(compute_elbo(X, t, alpha, beta, xi, mu, sigma))
        if display:
            print('%9d %15.11f' % (iteration, L[-1]))
        if abs(L[-1] - L[-2]) <= tol:
            break

    return {'alpha': alpha, 'beta': beta, 'mu': mu, 'sigma': sigma, 'L': L}