<jupyter_start><jupyter_text>Data Science 2: Advanced Topics in Data Science  Homework 3 - Bayesian Analysis**Harvard University****Spring 2020****Instructors**: <NAME>, <NAME>, & <NAME>   Homework 3 is due March 5th<jupyter_code>#PLEASE RUN THIS CELL
import requests
from IPython.core.display import HTML
styles = requests.get("https://raw.githubusercontent.com/Harvard-IACS/2018-CS109A/master/content/styles/cs109.css").text
HTML(styles)<jupyter_output><empty_output><jupyter_text>INSTRUCTIONS- To submit your assignment, please follow the instructions on Canvas.- Please restart the kernel and run the entire notebook again before you submit.  Please use the libraries below:<jupyter_code>import pandas as pd
import numpy as np
import pymc3 as pm
import matplotlib.pyplot as plt
%matplotlib inline 
from datetime import datetime as dt
from scipy.special import expit<jupyter_output><empty_output><jupyter_text>Contraceptive Usage by Bangladeshi Women In the second part of this assignment, you are provided with data sets `dataset_2_train.csv` and `dataset_2_test.csv`containing details of contraceptive usage by 1934 Bangladeshi women. There are four attributes for each woman, along with a label indicating if she uses contraceptives. The attributes include* district: identifying code for the district the woman lives in* urban: type of region of residence* living.children: number of living children* age-mean: age of the woman (in years, centered around mean)The women are grouped into 60 districts. The task is to build a classification model that can predict if a given woman uses contraceptives.  Part A: Bayesian Logistic Regression The first model we will fit to the contraceptives data is a varying-intercept logistic regression model, where the intercept varies by district.Prior distribution:$\beta_{0j} \sim N(\mu_0,\sigma_0^2),$ with $\mu_0 \sim N(0,10000)$ and $\frac{1}{\sigma^2_0} \sim \text{Gamma}(0.1,0.1)$.$\beta_1 \sim N(0,10000)$, $\beta_2\sim N(0,10000)$, $\beta_3 \sim N(0,10000)$Model for data:$Y_{ij} \sim \text{Bernoulli}(p_{ij})$$\text{logit}\:p_{ij} = \beta_{0j} + \beta_1 \times \text{urban} + \beta_2 \times \text{living.children} + \beta_3 \times \text{age-mean}$,where $Y_{ij}$ is 1 if woman $i$ in district $j$ uses contraceptives, and 0 otherwise, and where $i \in \{1,...,N\}$ and $j \in \{1,...,J\}$. $N$ is the number of observations in the data, and $J$ is the number of districts. The above notation assumes $N(\mu, \sigma^2)$ is a Normal distribution with mean $\mu$ and variance $\sigma^2$.  **A1** As a preprocessing step, map the district number 61 to the number 54 so that the districts are in order. Also, re-name district 60 to be district 0 so that the districts are numbered 0 through 59. **We use this numbering throughout the homework, and will grade using these district numbers**<jupyter_code>#your code here
train_df = pd.read_csv("../data/dataset_2_train.csv");
test_df = pd.read_csv("../data/dataset_2_test.csv");
train_df.info()
train_df.describe()
train_df['district'].sort_values().unique()
train_df['district'].replace({61:54, 60:0}, inplace=True)
train_df.rename(columns = {"living.children":"living_children"}, inplace = True) 
test_df['district'].replace({61:54, 60:0}, inplace=True)
test_df.rename(columns = {"living.children":"living_children"}, inplace = True) 
train_df['district'].sort_values().unique()
train_df.head()<jupyter_output><empty_output><jupyter_text>**A2** First, we'll verify that pymc3 can recover the hidden parameter values. To do so, we'll hard-code known values and simulate data from the model, then using that data we'll check if pymc3 can get back the parameters we hard-coded. If it does, we'll have hope that it can get the hidden parameter values that generated the real data.Use the following hard-coded values:  $\mu_0 = 2$$\sigma^2_0 = 1$$\beta_{0j} \sim N(\mu_0,\sigma_0^2)$ for all sixty districts$\beta_1 = 4$$\beta_2 = -3$$\beta_3 = -2$**Verify that these parameters give an overall average of 46 to 48 percent of subjects using contraceptives**<jupyter_code>#your code here
np.random.seed(109) #use this seed tto generate the same results

print('The urban mean is', train_df['urban'].mean())
print('The living_children mean is', train_df['living_children'].mean())
print('The age_mean mean is ', train_df['age_mean'].mean())
print('The current contraceptive_use mean is', train_df['contraceptive_use'].mean())
print('The calculated contraceptive_use mean is', expit(2 + 4*.34 + -3*1 + -2*.24))

alpha, sigma = 2, 1
beta = [4, -3, -2]

colNames = ('district', 'beta0', 'urban', 'living_children', 'age_mean', 'contraceptive_use')
sim_df = pd.DataFrame(columns = colNames)

for i in np.unique(test_df['district']):
    beta = np.random.normal(alpha, sigma , 1)
    #u_mean = np.mean(test_df['urban'][test_df['district'] == i]) # replaced with calculated mean above
    #l_mean = np.mean(test_df['living_children'][test_df['district'] == i]) # replaced with calculated mean above
    l_sd = np.std(test_df['living_children'][test_df['district'] == i])
    #a_mean = np.mean(test_df['age_mean'][test_df['district'] == i]) # replaced with calculated mean above
    a_sd = np.std(test_df['age_mean'][test_df['district'] == i])
    for j in test_df[test_df['district'] == i].iterrows():
        temp_df = pd.DataFrame(columns = colNames)
        temp_df['district'] = [i]
        temp_df['beta0'] = beta
        temp_df['urban'] = np.random.binomial(1, .34, size=1)
        temp_df['living_children'] = np.random.normal(1.21, l_sd, size=1) # adjusted to generated .47 contraceptive_use mean
        temp_df['age_mean'] = np.random.normal(.24, a_sd, size=1)
        sim_df = sim_df.append(temp_df)

sim_df['district'] = sim_df['district'].astype(int)
sim_df['urban'] = sim_df['urban'].astype(int)

sim_df['contraceptive_use'] = np.round(expit(np.asfarray(sim_df['beta0'] + 4*sim_df['urban'] + -3*sim_df['living_children'] + -2*sim_df['age_mean'])))

sim_df['contraceptive_use'] = sim_df['contraceptive_use'].astype(int)

print('The new contraceptive_use mean is',np.mean(sim_df['contraceptive_use']))

plt.hist(sim_df['contraceptive_use'])
plt.show()<jupyter_output>The urban mean is 0.30713547052740436
The living_children mean is 2.591520165460186
The age_mean mean is  -0.31062595656670683
The current contraceptive_use mean is 0.38986556359875907
The calculated contraceptive_use mean is 0.4700359482354283
The new contraceptive_use mean is 0.4705274043433299<jupyter_text>**A3** Fit the varying-intercept model specified above to your simulated data. Reasonably good convergence may require adjustments to the number of tuning samples and the target acceptance rate as suggested in pymc3's warnings.Note: It is wise to take advantage of a model's built-in `summary` method to inspect it after fitting<jupyter_code>#your code here

with pm.Model() as contraceptive_model:
    mu = pm.Normal('mu', mu=0, tau=10000)
    sigma = pm.Gamma('sigma', alpha=.001, beta=.001)
    
    beta0 = pm.Normal('beta0', mu=mu, sigma=sigma, shape=len(np.unique(sim_df['district'])))[sim_df['district']]
    beta1 = pm.Normal('beta1', mu=0, tau=10000)
    beta2 = pm.Normal('beta2', mu=0, tau=10000)
    beta3 = pm.Normal('beta3', mu=0, tau=10000)
    
    sP_i = pm.Deterministic('$sP_i$', pm.math.invlogit(beta0 + beta1*sim_df['urban'] + beta2*sim_df['living_children'] + beta3*sim_df['age_mean']))
    
    obs = pm.Bernoulli('obs', p=sP_i, shape=len(sim_df), observed=sim_df['contraceptive_use'])
    
    trace = pm.sample(1500, tune=1500, target_accept=0.95)
    
pm.summary(trace, var_names=['mu', 'sigma', 'beta0', 'beta1', 'beta2', 'beta3'])<jupyter_output>/usr/share/anaconda3/lib/python3.7/site-packages/theano/tensor/subtensor.py:2197: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  rval = inputs[0].__getitem__(inputs[1:])
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta3, beta2, beta1, beta0, sigma, mu]
Sampling 2 chains, 0 divergences: 100%|██████████| 6000/6000 [03:16<00:00, 30.55draws/s] 
/usr/share/anaconda3/lib/python3.7/site-packages/theano/tensor/subtensor.py:2197: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different re[...]<jupyter_text>**A4** Plot the trace plots of the MCMC sampler for the parameters $\mu_0, \sigma_0, \beta_1, \beta_2, \beta_3$. Based on these and the R-hat values, does it look like the samplers converged?<jupyter_code>#your code here
pm.traceplot(trace, var_names=('mu', 'sigma', 'beta1', 'beta2', 'beta3'))<jupyter_output>/usr/share/anaconda3/lib/python3.7/site-packages/theano/tensor/subtensor.py:2197: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  rval = inputs[0].__getitem__(inputs[1:])<jupyter_text>your response hereThe Bayseian logistic regression model we created to fit our simulated data was NOT able to capture the hard-coded values within the parameters, $\mu_0, \sigma_0, \beta_1, \beta_2, \beta_3$, distributions of our simulated data. In addition, based on the R-hat values, larger than 1.05, the samplers do not seem to have converged.  **A5** Plot histograms of the posterior distributions for the parameters $\beta_{0,9}$, $\beta_{0,19}$, $\beta_{0,29}$, ..., $\beta_{0,59}$. Are the actual parameters that you generated contained within these posterior distributions?<jupyter_code>#your code here
#pm.plot_posterior(trace, var_names=('beta0'))
pm.plot_posterior(trace, var_names='beta0', coords={'beta0_dim_0': [9, 19, 29, 39, 49, 59]}, kind='hist')

pm.plot_forest(trace,
               kind='forestplot',
               var_names=['beta0'],
               combined=True,
               ridgeplot_overlap=60,
               figsize=(20, 5))<jupyter_output><empty_output><jupyter_text>your response hereThe posterior distributions for the parameters $\beta_{0,9}$, $\beta_{0,19}$, $\beta_{0,29}$, ..., $\beta_{0,59}$ do not contain the actual parameters we generated. The means are very close to 0.  **A6** We now fit our model to the actual data. Fit the varying-intercept model to the real training data.<jupyter_code>#your code here

with contraceptive_model:
    
    tP_i = pm.Deterministic('$tP_i$', pm.math.invlogit(beta0 + beta1*train_df['urban'] + beta2*train_df['living_children'] + beta3*train_df['age_mean']))
    
    obs_cu = pm.Bernoulli('obs_cu', p=tP_i, shape=len(train_df), observed=train_df['contraceptive_use'])

    trace = pm.sample(1500, tune=1500, target_accept=0.95)

    print(trace.varnames)
    print('DONE')<jupyter_output>Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
/usr/share/anaconda3/lib/python3.7/site-packages/theano/tensor/subtensor.py:2197: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  rval = inputs[0].__getitem__(inputs[1:])
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta3, beta2, beta1, beta0, sigma, mu]
Sampling 2 chains, 0 divergences: 100%|██████████| 6000/6000 [01:45<00:00, 56.66draws/s] 
The estimated number of effective samples is smaller than 200 for some parameters.<jupyter_text>**A7** Check the convergence by examining the trace plots and R-hats, as you did with the simulated data. What do you observe?<jupyter_code>#your code here
pm.traceplot(trace, var_names=('mu', 'sigma', 'beta1', 'beta2', 'beta3'))<jupyter_output><empty_output><jupyter_text>your response here  **A8** Based on the posterior means, which district has the highest base rate of contraceptive usage (independent of other factors like urban population)? Which district has the lowest?<jupyter_code>#your code here
beta0_sum = pm.summary(trace, var_names=['beta0'])
beta0_sum[(beta0_sum['mean']==beta0_sum['mean'].min()) | (beta0_sum['mean']==beta0_sum['mean'].max())]

pm.plots.plot_posterior(trace, var_names = 'beta0', coords={'beta0_dim_0': [27, 15]})

pm.plot_forest(trace,
               kind='forestplot',
               var_names=['beta0'],
               combined=True,
               ridgeplot_overlap=60,
               figsize=(20, 5))<jupyter_output><empty_output><jupyter_text>your response here  **A9** What are the posterior means of $\mu_0$ and $\sigma^2_0$? Do these values offer any evidence in support of or against the varying-intercept model, compared to a model with a single intercept value for all districts or a model with an intercept of 0 for all districts?<jupyter_code>#your code here
pm.plot_posterior(trace, var_names=('mu', 'sigma'))<jupyter_output><empty_output><jupyter_text>your response here  Part B: Varying-Coefficients ModelThe next model we will fit to the contraceptives data is a varying-coefficients logistic regression model, where the coefficients on living.children, age-mean, and urban vary by district.$\beta_{0j} \sim N(\mu_0,\sigma_0^2),$ with $\mu_0 \sim N(0,10000)$ and $\frac{1}{\sigma_0^2} \sim \text{Gamma}(0.1,0.1)$.$\beta_{1j} \sim N(0,\sigma_1^2)$, with $\frac{1}{\sigma^2_1} \sim \text{Gamma}(0.1,0.1)$$\beta_{2j}\sim N(0,\sigma_2^2)$, with $\frac{1}{\sigma^2_2}\sim \text{Gamma}(0.1,0.1)$$\beta_{3j} \sim N(0,\sigma_3^2)$, with $\frac{1}{\sigma^2_3} \sim \text{Gamma}(0.1,0.1)$Model for data:$Y_{ij} \sim \text{Bernoulli}(p_{ij})$$\text{logit}\:p_{ij} = \beta_{0j} + \beta_{1j} \times \text{urban} + \beta_{2j} \times \text{living.children} + \beta_{3j} \times \text{age-mean}$,  **B1** Fit the model to the real data. Plot posterior means for each of the coefficients to the predictors. What do you conclude from these graphs?<jupyter_code>#your code here
with pm.Model() as contraceptive_model:
    mu = pm.Normal('mu', mu=0, tau=10000)
    
    sigma0 = pm.Gamma('sigma0', alpha=.001, beta=.001)
    sigma1 = pm.Gamma('sigma1', alpha=.001, beta=.001)
    sigma2 = pm.Gamma('sigma2', alpha=.001, beta=.001)
    sigma3 = pm.Gamma('sigma3', alpha=.001, beta=.001)
    
    beta0 = pm.Normal('beta0', mu=mu, sigma=sigma0, shape=len(np.unique(train_df['district'])))[train_df['district']]
    beta1 = pm.Normal('beta1', mu=0, sigma=sigma1, shape=len(np.unique(train_df['district'])))[train_df['district']]
    beta2 = pm.Normal('beta2', mu=0, sigma=sigma2, shape=len(np.unique(train_df['district'])))[train_df['district']]
    beta3 = pm.Normal('beta3', mu=0, sigma=sigma3, shape=len(np.unique(train_df['district'])))[train_df['district']]
    
    p_i = pm.Deterministic('$P_i$', pm.math.invlogit(beta0 + beta1*train_df['urban'] + beta2*train_df['living_children'] + beta3*train_df['age_mean']))
    
    obs = pm.Bernoulli('obs', p=p_i, shape=len(train_df), observed=train_df['contraceptive_use'])
    
    trace = pm.sample(2000, tune=2000, target_accept=0.95)
    
    print(trace.varnames)    
    print('DONE')

pm.summary(trace, var_names=['mu', 'sigma0', 'sigma1', 'sigma2', 'sigma3'])

pm.traceplot(trace, var_names=('mu', 'sigma0', 'sigma1', 'sigma2', 'sigma3'))

pm.plot_forest(trace,
               kind='forestplot',
               var_names=['mu', 'sigma0', 'sigma1', 'sigma2', 'sigma3'],
               combined=True,
               ridgeplot_overlap=5,
               figsize=(20, 2))
beta0_sum = pm.summary(trace, var_names=['beta0'])
beta0_sum[(beta0_sum['mean']==beta0_sum['mean'].min()) | (beta0_sum['mean']==beta0_sum['mean'].max())]

pm.plot_posterior(trace, var_names='beta0', coords={'beta0_dim_0': [11, 56]}, kind='hist')

pm.plot_forest(trace,
               kind='forestplot',
               var_names=['beta0'],
               combined=True,
               ridgeplot_overlap=60,
               figsize=(20, 5))

beta1_sum = pm.summary(trace, var_names=['beta1'])
beta1_sum[(beta1_sum['mean']==beta1_sum['mean'].min()) | (beta1_sum['mean']==beta1_sum['mean'].max())]

pm.plot_posterior(trace, var_names='beta1', coords={'beta1_dim_0': [52, 14]}, kind='hist')

pm.plot_forest(trace,
               kind='forestplot',
               var_names=['beta1'],
               combined=True,
               ridgeplot_overlap=60,
               figsize=(20, 5))

beta2_sum = pm.summary(trace, var_names=['beta2'])
beta2_sum[(beta2_sum['mean']==beta2_sum['mean'].min()) | (beta2_sum['mean']==beta2_sum['mean'].max())]

pm.plot_posterior(trace, var_names='beta2', coords={'beta2_dim_0': [54, 30]}, kind='hist')

pm.plot_forest(trace,
               kind='forestplot',
               var_names=['beta2'],
               combined=True,
               ridgeplot_overlap=60,
               figsize=(20, 5))

beta3_sum = pm.summary(trace, var_names=['beta3'])
beta3_sum[(beta3_sum['mean']==beta3_sum['mean'].min()) | (beta3_sum['mean']==beta3_sum['mean'].max())]

pm.plot_posterior(trace, var_names='beta3', coords={'beta3_dim_0': [54, 25]}, kind='hist')

pm.plot_forest(trace,
               kind='forestplot',
               var_names=['beta3'],
               combined=True,
               ridgeplot_overlap=60,
               figsize=(20, 5))<jupyter_output><empty_output><jupyter_text>your response here  **B2** Use all of the information you've gleaned thus far to build a final Bayesian logistic regression classifier on the train set. Then, use your model to make predictions on the test set. Report your model's classification percentage.<jupyter_code>#your code here
with pm.Model() as contraceptive_model:
    district = pm.Data('district', train_df['district'])
    urban = pm.Data('urban', train_df['urban'])
    living_children = pm.Data('living_children', train_df['living_children'])
    age_mean = pm.Data('age_mean', train_df['age_mean'])
    contraceptive_use = pm.Data('contraceptive_use', train_df['contraceptive_use'])
    
    mu = pm.Normal('mu', mu=0, tau=10000)
    
    sigma0 = pm.Gamma('sigma0', alpha=.001, beta=.001)
    #sigma1 = pm.Gamma('sigma1', alpha=.001, beta=.001)
    #sigma2 = pm.Gamma('sigma2', alpha=.001, beta=.001)
    #sigma3 = pm.Gamma('sigma3', alpha=.001, beta=.001)
    
    beta0 = pm.Normal('beta0', mu=mu, sigma=sigma0, shape=len(np.unique(train_df['district'])))[train_df['district']]
    beta1 = pm.Normal('beta1', mu=0, tau=10000)#, sigma=sigma1, shape=len(np.unique(train_df['district'])))[train_df['district']]
    beta2 = pm.Normal('beta2', mu=0, tau=10000)#, sigma=sigma2, shape=len(np.unique(train_df['district'])))[train_df['district']]
    beta3 = pm.Normal('beta3', mu=0, tau=10000)#, sigma=sigma3, shape=len(np.unique(train_df['district'])))[train_df['district']]
    
    #p_i = pm.Deterministic('$P_i$', pm.math.invlogit(beta0 + beta1*train_input['urban'] + beta2*train_input['living_children'] + beta3*train_input['age_mean']))
    p_i = pm.Deterministic('$P_i$', pm.math.invlogit(beta0 + beta1*urban + beta2*living_children + beta3*age_mean))
    
    obs = pm.Bernoulli('obs', p=p_i, shape=len(train_df['contraceptive_use']), observed=contraceptive_use)
    
    trace = pm.sample(1500, tune=1500, target_accept=0.95)
    
    print(trace.varnames)  
    print('DONE')

with contraceptive_model:
    # Switch out the observations and use `sample_posterior_predictive` to predict
    district.set_value(test_df['district'])
    urban.set_value(test_df['urban'])
    living_children.set_value(test_df['living_children'])
    age_mean.set_value(test_df['age_mean'])
    contraceptive_use.set_value(test_df['contraceptive_use'])
    
    post_pred = pm.sample_posterior_predictive(trace, samples=len(test_df['contraceptive_use']))
    
    print('DONE')

prediction = post_pred['obs'][:,0]
print('The true contraceptive use mean is', test_df['contraceptive_use'].mean())
print('The predicted contraceptive use mean is', prediction.mean())

plt.hist(test_df['contraceptive_use'])
plt.title('Histogram of Contraceptive Use (test)')
plt.show()

plt.hist(prediction)
plt.title('Histogram of Contraceptive Use (prediction)')
plt.show()

# https://www.geeksforgeeks.org/confusion-matrix-machine-learning/
from sklearn.metrics import confusion_matrix, accuracy_score, classification_report

c_matrix = confusion_matrix(test_df['contraceptive_use'], prediction)

print('Confusion Matrix :')
print(c_matrix)
print('Accuracy Score :',accuracy_score(test_df['contraceptive_use'], prediction))
print('Report : ')
print(classification_report(test_df['contraceptive_use'], prediction))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15001672160909654, 'fraction_numerical': 0.040179637857722994, 'mean_word_length': 5.040981240981241, 'pattern_counts': {'":': 2, '<': 54, '<?xml version=': 0, '>': 52, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '71918', 'license': 'MIT', 'n_tokens_mistral': 7632, 'n_tokens_neox': 6956, 'n_words': 1774, 'path': 'content/homeworks/hw03/notebook/cs109b_hw3_edited.ipynb', 'repo_name': 'medeirosf/2020-CS109B', 'size': '20955'}