In [8]:

    

import numpy as np
import test
import matplotlib.pyplot as plt


    

In [5]:

    

pima = np.genfromtxt('pima-indians-diabetes.data', delimiter=',')


    

In [6]:

    

# Load data
X = np.concatenate((np.ones((pima.shape[0],1)),pima[:,0:8]), axis=1)
Y = pima[:,8]

Xs = (X - np.mean(X, axis=0))/np.concatenate((np.ones(1),np.std(X[:,1:], axis=0)))
n, p = X.shape

nsample = 1
nbatch = 768
M = np.identity(p)
C = 0 * np.identity(p)
eps = 0.1
m = 10
V = 0 * np.identity(p)
theta = np.zeros(p)


    

In [9]:

    

phi=1
test.U(np.ones(p)*.1,Y,X)


    

    
Out[9]:





16642.14181768343

In [ ]:

    

# HMC - Unscaled
nsample = 10000
m = 20
eps = .0001
#theta = np.zeros(p)
theta = beta_true_unscale.copy()
phi = 0.01

np.random.seed(2)
samples = np.zeros((nsample, p))
u = np.zeros(nsample)
for i in range(nsample):
    theta = hmc(Y, X, gradU, M, eps, m, theta, C, V)
    samples[i] = theta
    u[i] = U(theta, Y, X)
    
np.mean(samples, axis=0) - beta_true_unscale