In [4]:

    

import pandas as pd
import pymc3 as pm
import matplotlib.pyplot as plt
import numpy as np
import theano.tensor as t
from scipy.stats import mode
def tinvlogit(x):
    return t.exp(x) / (1 + t.exp(x))

df=pd.read_csv('/Users/arrigal001/Desktop/thads2013n.txt',sep=',') 
df=df[df['BURDEN']>0] 
df=df[df['AGE1']>0] 
df['OWN']=[1 if obj=='2' else 0 for obj in df['OWNRENT']]


    

    




/anaconda3/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters

In [9]:

    

with pm.Model() as model: 
    # Define priors
    intercept = pm.Normal('Intercept', 0, sd=10)
    x_coeff = pm.Normal('x', 0, sd=10)
    price_coef = pm.Normal('price', 0, sd=10)
    
    # Define likelihood
    likelihood = pm.Bernoulli('y', 
                              pm.math.sigmoid(intercept+x_coeff*df['BEDRMS']+price_coef*df['COSTMED']),
                              observed=df['OWN'])
    WTP=pm.Deterministic('WTP',-x_coeff/price_coef)

    # Inference!
#     trace = pm.sample(3000)
    advi = pm.ADVI()
    approx = advi.fit(20000)
    
# pm.traceplot(trace)


    

    




Average Loss = 2.5601e+06: 100%|██████████| 20000/20000 [00:30<00:00, 648.07it/s]
Finished [100%]: Average Loss = 2.5601e+06

In [11]:

    

plt.plot(approx.hist)


    

    
Out[11]:





[<matplotlib.lines.Line2D at 0x1c123e7ac8>]






    

In [43]:

    

import seaborn as sns
ax = sns.kdeplot(trace['Intercept'], label='NUTS');
sns.kdeplot(approx.sample(10000)['Intercept'], label='ADVI');


    

In [15]:

    

pm.traceplot(trace)


    

    
Out[15]:





array([[<matplotlib.axes._subplots.AxesSubplot object at 0x1c1a2183c8>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1c1257c9b0>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x1c12f01198>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1c12f29e48>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x1c12f539e8>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1c12f7c518>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x1c12fa6080>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1c12fbc518>]],
      dtype=object)






    

In [41]:

    

trace_advi = approx.sample(draws=5000)
pm.traceplot(trace_advi)


    

    
Out[41]:





array([[<matplotlib.axes._subplots.AxesSubplot object at 0x1c1b26a710>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1c1b24c240>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x1c1b249eb8>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1c12669b70>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x1c1264e6a0>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1c12664cf8>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x1c1b19b978>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1c126034e0>]],
      dtype=object)






    

In [1]:

    

import pandas as pd
import pymc3 as pm
import matplotlib.pyplot as plt
import numpy as np

def logistic(x, b, noise=None):
    L = x.T.dot(b)
    if noise is not None:
        L = L+noise
    return 1/(1+np.exp(-L))

x1 = np.linspace(-10., 10, 10000)
x2 = np.linspace(0., 20, 10000)
bias = np.ones(len(x1))
X = np.vstack([x1,x2,bias]) # Add intercept
B =  [-10., 2., 1.] # Sigmoid params for X + intercept

# Noisy mean
pnoisy = logistic(X, B, noise=np.random.normal(loc=0., scale=0., size=len(x1)))
# dichotomize pnoisy -- sample 0/1 with probability pnoisy
# y = np.random.binomial(1., pnoisy)
y = [1 if x > 0.5 else 0 for x in pnoisy]


    

    




/anaconda3/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters

In [2]:

    

with pm.Model() as model: 
    # Define priors
    intercept = pm.Normal('Intercept', 0, sd=10)
    x1_coef = pm.Normal('x1', 0, sd=10)
    x2_coef = pm.Normal('x2', 0, sd=10)
    
# Define likelihood
    likelihood = pm.Bernoulli('y', 
                              pm.math.sigmoid(intercept+x1_coef*X[0]+x2_coef*X[1]),
                              observed=y)

#     trace = pm.sample(3000)
#     advi = pm.ADVI()
#     approx = advi.fit(150000)
    approx = pm.fit(60000, method='advi')


    

    




Average Loss = 140.89: 100%|██████████| 60000/60000 [00:54<00:00, 1106.56it/s]
Finished [100%]: Average Loss = 140.89

In [3]:

    

trace_advi = approx.sample(draws=900)
pm.traceplot(trace_advi)


    

    
Out[3]:





array([[<matplotlib.axes._subplots.AxesSubplot object at 0x1c152819e8>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1c15745e10>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x1c1563c2b0>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1c15f3d470>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x1c1500ccf8>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1c157e6f60>]],
      dtype=object)






    

In [147]:

    

pm.traceplot(trace)


    

    
Out[147]:





array([[<matplotlib.axes._subplots.AxesSubplot object at 0x1c24d423c8>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1c24d72320>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x1c24da4e48>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1c24dc9390>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x1c24df7e80>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1c24e2ba90>]],
      dtype=object)






    

In [158]:

    

sns.distplot(trace['Intercept'])
sns.distplot(trace_advi['Intercept'])


    

    
Out[158]:





<matplotlib.axes._subplots.AxesSubplot at 0x1c243e0da0>






    

In [7]:

    




    

    
Out[7]:





0

In [ ]: