add bbvi

Author: saurabhmathur96

bbvi.py

@@ -1 +1,79 @@
 # posterior inference by black-box variational inference
+
+
+import autograd.numpy as np
+from autograd import grad
+import autograd.scipy as scipy
+import autograd.scipy.stats.norm as norm
+
+from autograd.misc.optimizers import adam, sgd
+
+def sigmoid(x):
+  return 0.5*(np.tanh(x/2)+1)
+
+def predict(w, x):
+  return sigmoid(np.dot(x, w))
+
+def log_sigmoid(x):
+  a = np.array([np.zeros_like(x), -x])
+  return -scipy.special.logsumexp(a, axis=0)
+
+def log_joint(w, x, y, alpha=0.1):
+  log_prior = alpha*np.sum(w**2)
+
+  score = np.dot(x,  w)
+  logp0 = log_sigmoid(score)
+
+  logp1 =  -score + logp0
+  log_likelihood = np.sum(y*logp0  + (1-y)*logp1, axis=0)
+
+  return log_likelihood - log_prior
+
+
+
+
+
+def bbvi(x, y, params, log_joint, T=25):
+  ''' Perform black-box variational inference
+    with q = diagonal gaussian
+  '''
+  def objective(params, i=0):
+    D = len(params) // 2
+    mu, log_sigma2 =  params[:D], params[D:]
+
+    entropy = 0.5*D*(1 + np.log(2*np.pi)) + np.sum(log_sigma2)
+    
+    # samples = np.random.randn(T, D) * np.exp(0.5*log_sigma2) + mu
+    logp = 0
+    sample = mu + np.random.randn(T, D) * np.exp(0.5*log_sigma2) 
+
+    for t in range(T):
+      logp += log_joint(sample[t], x, y)
+    return -(logp/T + entropy)
+
+  gradient = grad(objective)
+  return objective, gradient
+
+
+
+
+x = np.array([[0.52, 1.12,  0.77],
+               [0.88, -1.08, 0.15],
+               [0.52, 0.06, -1.30],
+               [0.74, -2.49, 1.39]])
+
+y = np.array([True, True, False, True])
+
+x = np.hstack([np.ones(( len(x),1)), x])
+
+params = np.zeros(4+4) 
+objective, gradient = bbvi(x, y, params, log_joint, T=100)
+print (objective(params))
+print (params)
+print (predict(params[:4], x))
+
+params = adam(gradient, params, step_size=0.01, num_iters=500)
+
+print (objective(params), objective(params))
+print (np.exp(params[4:]))
+print (predict(params[:4], x))

laplace.py

@@ -1,4 +1,5 @@
 # posterior inference by laplace approximation
+from autograd.misc.optimizers import adam, sgd
 from autograd import grad
 import autograd.numpy as np
 import autograd.scipy as scipy
@@ -12,23 +13,25 @@ def predict(w, x):
 
 def log_sigmoid(x):
   a = np.array([np.zeros_like(x), -x])
-  return -scipy.special.logsumexp(a)
+  return -scipy.special.logsumexp(a, axis=0)
 
 def nll_loss(w, x, y, alpha=None):
   score = np.dot(x, w)
 
-  logp1 = log_sigmoid(score)
-  logp0 = -score + logp1
-  loss = y*logp0  + (1-y)*logp1
-  reg = alpha*np.dot(w, w) if alpha else 0
-  return sum(loss) + reg
+  logp0 = log_sigmoid(score)
+
+  logp1 = -score+logp0
+  
+  loss = -np.sum(y*logp0  + (1-y)*logp1)
+  reg = alpha*np.sum(w**2) if alpha else 0
+  return  loss + reg
 
 def compute_precision(x, y, w, alpha):
   d = np.size(x, 1)
   y_hat = predict(w, x)
   R = np.diag(y_hat*(1 - y_hat))
-  precision = alpha * np.eye(d) + x.T.dot(R).dot(x)
-  return precision + 1e-9*np.eye(d)
+  precision = 1e-9*np.eye(d) + alpha * np.eye(d) + x.T.dot(R).dot(x)
+  return precision 
 
 def predict_mc(mu, sigma, x, T=100):
   ps = []
@@ -44,30 +47,34 @@ def predict_var(mu, sigmainv, x):
   kappa = np.sqrt(1 + sigma2_a*np.pi*.125)
   return sigmoid(mu_a/kappa) 
 
+
 x = np.array([[0.52, 1.12,  0.77],
                [0.88, -1.08, 0.15],
                [0.52, 0.06, -1.30],
                [0.74, -2.49, 1.39],
                [0.52, 1.12,  0.77]])
+
 y = np.array([True, True, False, True, False])
 
 
 x = np.hstack([np.ones(( len(x),1)), x])
-training_loss = lambda w: nll_loss(w, x, y, alpha=1)
+training_loss = lambda w, i: nll_loss(w, x, y, alpha=0.1)
 g = grad(training_loss)
-w = np.array([1, 0, 0, 0], dtype=np.float)
-print("Initial loss:", training_loss(w))
-for i in range(100):
-    w -= g(w) * 0.01
-print("Trained loss:", training_loss(w))
+w = np.array([1, 1, 1, 1], dtype=np.float)
+print("Initial loss:", training_loss(w, 0))
+#for i in range(100):
+#    w -= g(w) * 0.01
+
+w = sgd(g, w)
+print("Trained loss:", training_loss(w, 0))
 
 pred = predict(w, x) > 0.5
 
 print (y.astype(int))
-print (predict(w, x) )
+print ('ml', predict(w, x) )
 
 
-sigmainv = compute_precision(x,y,w,alpha=1)
+sigmainv = compute_precision(x,y,w,alpha=0.1)
 
-print (predict_var(w, sigmainv, x))
-print (predict_mc(w, np.linalg.inv(sigmainv), x))
+print ('var', predict_var(w, sigmainv, x))
+print ('mc', predict_mc(w, np.linalg.inv(sigmainv), x))