In [1]:

    

%load_ext autoreload
%autoreload 2
from autograd import grad
import autograd.numpy as np
import autograd.numpy.random as npr
import autograd.scipy.special as sp

from gamma_def import *
from gamma_def_rejection import *
from gamma_def_generalized import *
from gamma_def_score import *


import matplotlib.pyplot as plt
%matplotlib inline


    

In [3]:

    

# Initialize
x = np.loadtxt('data/nips_train.csv')

N = x.shape[0]
D = x.shape[1]
alphaz = 0.1
K = np.array([100,40,15])

params = np.load('results/Nips_Eta0.05_B2_corrTrue_K1_100_K2_40_K3_15_params_R.npy')


    

In [ ]:

    

%timeit grep_gradient(params[:,0],params[:,1],x,K,alphaz)
%timeit reparam_gradient(params[:,0],params[:,1],x,K,alphaz,corr=True,B=4)


    

In [9]:

    

# Estimate variance at initial step
np.random.seed(0)
sigma = 0.1
params[:,0] = 0.5+sigma*np.random.normal(size=N*np.sum(K)+D*K[0]+K[1]*K[0]+K[2]*K[1])
params[:,1] = sigma*np.random.normal(size=N*np.sum(K)+D*K[0]+K[1]*K[0]+K[2]*K[1])
params = np.log(1.+np.exp(params))
params[params[:,0]<1.,0] = 1.

S = 10
grep = np.zeros((S,params.shape[0],params.shape[1]))
rej1 = np.zeros((S,params.shape[0],params.shape[1]))
rej2 = np.zeros((S,params.shape[0],params.shape[1]))
for s in range(S):
    grep[s,:,:] = grep_gradient(params[:,0],params[:,1],x,K,alphaz)
    rej1[s,:,:] = reparam_gradient(params[:,0],params[:,1],x,K,alphaz,corr=True,B=0)
    rej2[s,:,:] = reparam_gradient(params[:,0],params[:,1],x,K,alphaz,corr=True,B=4)


    

In [ ]:

    

print np.min(np.var(rej1,axis=0)),np.median(np.var(rej1,axis=0)),np.max(np.var(rej1,axis=0))
print np.min(np.var(rej2,axis=0)),np.median(np.var(rej2,axis=0)),np.max(np.var(rej2,axis=0))
print np.min(np.var(grep,axis=0)),np.median(np.var(grep,axis=0)),np.max(np.var(grep,axis=0))


    

In [ ]:

    

# Estimate variance after 2600 iterations
np.random.seed(0)
params = np.load('results/Nips_Eta0.05_B2_corrTrue_K1_100_K2_40_K3_15_params_R.npy')
S = 10
grep = np.zeros((S,params.shape[0],params.shape[1]))
rej1 = np.zeros((S,params[ind,:].shape[0],params[ind,:].shape[1]))
rej2 = np.zeros((S,params.shape[0],params.shape[1]))
for s in range(S):
    grep[s,:,:] = grep_gradient(params[:,0],params[:,1],x,K,alphaz)
    rej1[s,:,:] = reparam_gradient(params[:,0],params[:,1],x,K,alphaz,corr=True,B=1.)
    rej2[s,:,:] = reparam_gradient(params[:,0],params[:,1],x,K,alphaz,corr=True,B=4)


    

In [ ]:

    

print np.min(np.var(rej1,axis=0)),np.median(np.var(rej1,axis=0)),np.max(np.var(rej1,axis=0))
print np.min(np.var(rej2,axis=0)),np.median(np.var(rej2,axis=0)),np.max(np.var(rej2,axis=0))
print np.min(np.var(grep,axis=0)),np.median(np.var(grep,axis=0)),np.max(np.var(grep,axis=0))


    

In [ ]: