<jupyter_start><jupyter_text>Introduction to Probabilistic Programming Probabilistic programming is a method for Bayesian statistical analysis. Bayesian methods, in contrast to Frequentist statistics, do not believe that the data represents the entire truth about the world. Instead, Bayesian methods are model based meaning there is a prior probability which is updated with more information. This is captured in the famous Bayes' Rule.$$P(Y | X) = \frac{P(X | Y) * P (Y)}{P(X)}$$Where * $P(Y | X)$: posterior probability of Y conditioned on the value of X - what we want to figure out, * $(P X | Y)$: likelihood - the observed data* $P(Y)$: prior - probability of Y given no observations, also known as the base rate* $P(X)$: marginal probability of X - total probability of X for all values of YIn the case where both X and Y are Bernoulli variables (can only take on True / False values), then the marginal probability of X is $$P(X) = P(X | Y) * P(Y) + P(X | \neg Y) * P(\neg Y)$$where $\neg$ means "not". I don't tend to get much value from equations, so let's plug in some numbers into an example of determining if a transaction is fraudulent considering only whether or not the amount of the transactions > $500.Here are the relevant pieces:$$P(\text{Fraudulent} | \text{Amount > \$500}) = \frac{P(\text{Amount > \$500} | \text{Fraudulent}) * P(\text{Fraudulent})}{P(\text{Amount > \$500})}$$* $P(\text{Fraudulent} | \text{Amount > \$500})$: unknown posterior- what we want to find out* $P(\text{Amount > \$500} | \text{Fraudulent}) = 0.6$: likelihood based on past observations* $P(\text{Fraudulent}) = 0.01)$: prior probability of a fraudulent transactionThe marginal probability of a transaction greater than \$500 is must take into account both the situation where the transaction is fraudulent and the situation where the transaction is not fraudulent. * $P(\text{Amount} > \$500) = P(\text{Amount > \$500} | \text{Fraudulent}) * P(\text{Fraudulent}) + P(\text{Amount > \$500} | \neg \text{Fraudulent}) * P(\neg \text{Fraudulent})$This means we need $P(\text{Amount > \$500} | \neg \text{Fraudulent})$ which also must be found from the historical data. Let's say this is 0.2 We now have all the pieces and can carry out the calculation. $$P(\text{Fraudulent} | \text{Amount > \$500}) = \frac{0.6 * 0.01}{(0.6 * 0.01) + (0.2 * 0.99)} = 0.0294 = 2.94\%$$ This is the _posterior_ probability of a fraudulent transaction _conditioned_ on the amount of the transaction greater than \$500. The model we created for fraud is extremely simple - there is only one cause for fraud, the amount of the transaction. What happens when we have multiple causes, each of which in turn may have their own cause? For this, we turn to Bayesian graphical models which allow us to show the cause and effect relationship between variables as well as the associated probabilities.  Bayesian Graphical ModelThe model below is a slightly more complex simulation of fraud. Here, fraud has two causes, if the transaction is greater than \$500, and if the transaction was made online. For now, we are still using only Bernoulli (True / False) variables in our model. A variable that does not have any causes has a prior probability as indicated by the blue tables. Variables that are dependent on other variables have a conditional probability table in orange. The conditional probability table specifies the probability of the dependent variable for all the values of its _parent_ variables.  ![](../images/fraud_model.png)Once we have the model, we can figure out the posterior probability of any variable given any combination of the other variables. For example, to find the marginal probability of an online transaction, we need to sum up the probability for each combination of the parent variables.$$P(Online) = P(Online | Location, Gender) * P(Location) * P(Gender) + P(Online| Location, \neg Gender) * P(Location) * P(\neg Gender) + P(Online | \neg Location, Gender) * P(\neg Location) * P(Gender) + P(Online | \neg Location, \neg Gender) * P(\neg Location) * P(\neg Gender)$$Plugging in the numbers we get: $P(Online) = 0.25 * 0.2 * 0.5 + 0.45 * 0.2 * 0.5 + 0.3 * 0.8 * 0.5 + 0.6 * 0.8 * 0.5 = 0.43$That's a lot of work! We can do a similar operation for a fraudulent transaction and compute the marginal probability:$$P(Fraud) = 0.45 * 0.1 * 0.43 + 0.2 * 0.1 * 0.57 + 0.1 * 0.9 * 0.43 + 0.05 * 0.9 * 0.57 = 0.10518$$Here we see that the probability of a fraudulent transaction knowing known of the other variables is 10.5%. We can also set variables to be fixed and then compute the probabilities. If we set the amount to be less than \$500, then the marginal probability of a fraudulent transaction is $$P(Fraud | Amount) = 0.1 * 1.0 * 0.43 + 0.05 * 1.0 * 0.57 = 0.07635$$The probability we compute has dropped because we have additional evidence. Computing these probabilities by hand is not very fun! We could make a computer program to solve the probabilities for us, which works when we have a limited number of discrete (can only take on a limited number of options) variables. So, for this example currently, we can solve this problem analytically - meaning exactly. However, when we have more variables that can take on any values - continuous variables as almost any quantity in real life is! - then we need an approximate method. The most popular approximate method is Markov Chain Monte Carlo, which draws samples from the posterior in order to approximate the posterior. We set up a probability model, and then have an algorithm sample (draw examples) thousands of times from the model according to the probabilities (prior and conditional) that we have defined. The crazy part is, this method works extremely well! The sampled posterior converges on the true posterior given enough samples. Moreover, we can set any variables as observed and find the probability of any other variable. We can set a transaction to be fraudulent, and then find the probability that the location is San Francisco. This is known as the inverse probability problem. That's enough explanation for now. Let's work through implementing a Bayesian graphical model in PyMC3 and then performing inference on the model.  Bayesian Model in PyMC3PyMC3 is a Python library for performing Bayesian inference on graphical models. It includes a number of MCMC methods for sampling from the posterior including the most efficient, the No-U-Turn Sampler (NUTS). We can translate our above graphical model into code using the following.<jupyter_code>import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
%matplotlib inline

import pymc3 as pm
import theano.tensor as tt

N_SAMPLES = 1000
N_CHAINS = 3
with pm.Model() as model:
        
    # Bernoulli random variables with priors
    gender = pm.Bernoulli('gender', 0.5)
    location = pm.Bernoulli('location', 0.2)

    # Bernoulli deterministic variable with conditional probability 
    online = pm.Bernoulli('online', tt.switch(tt.eq(location, 1), tt.switch(tt.eq(gender, 1), 0.25, 0.45),
                                                                  tt.switch(tt.eq(gender, 1), 0.3, 0.6)))

    # Bernoulli random variable with prior
    amount = pm.Bernoulli('amount', 0.1)

    # Bernoulli deterministic variable with conditional probability
    fraudulent = pm.Bernoulli('fraudulent', tt.switch(tt.eq(amount, 1), tt.switch(tt.eq(online, 1), 0.45, 0.2),
                                                                        tt.switch(tt.eq(online, 1), 0.1, 0.05)))
pm.model_to_graphviz(model)
model.vars
model.test_point
with model:
    trace = pm.sample()
trace['fraudulent'].mean()
trace['online'].mean()
pm.plot_posterior(trace)
with model:
    x = pm.model_graph
pm.summary(trace)
with model:
    pm.model_to_graphviz()
with pm.Model() as model:
        
    # Bernoulli random variables with priors
    gender = pm.Bernoulli('gender', 0.5)
    location = pm.Bernoulli('location', 0.2)

    # Bernoulli deterministic variable with conditional probability 
    online = pm.Bernoulli('online', tt.switch(tt.eq(location, 1), tt.switch(tt.eq(gender, 1), 0.25, 0.45),
                                                                  tt.switch(tt.eq(gender, 1), 0.3, 0.6)))

    # Bernoulli random variable with prior
    amount = pm.Normal('amount', mu = 300, sd = 100)

    # Bernoulli deterministic variable with conditional probability
    fraudulent = pm.Bernoulli('fraudulent', tt.switch(tt.eq(amount, 1), tt.switch(tt.eq(online, 1), 0.45, 0.2),
                                                                        tt.switch(tt.eq(online, 1), 0.1, 0.05)))
with model:
    trace = pm.sample()
    
pm.traceplot(trace, varnames=['amount']);
import theano

values = theano.shared(np.array([100, 150, 200, 250]))
with pm.Model() as model:
    c = pm.Categorical('c', [0.1, 0.2, 0.3, 0.4])
    value = values[c]
    value = pm.Deterministic('value_det', value)
    trace = pm.sample()
trace['value_det'].mean()
trace['c'].mean()
trace['value_det']
value
pm.autocorrplot(trace)
trace.varnames
model.free_RVs
model.deterministics
trace.get_sampler_stats('location')
trace['fraudulent_det'].mean()
model.unobserved_RVs
pm.Deterministic(
model.deterministics
model.vars
model.disc_vars
def make_model(gender_obs=None, location_obs=None, online_obs=None,
               amount_obs=None, fraudulent_obs=None):
    """Make a model for the fraudulent transaction problem."""
    
    with pm.Model() as model:
        
        # Bernoulli random variables with priors
        gender = pm.Bernoulli('gender', 0.5, observed = gender_obs)
        location = pm.Bernoulli('location', 0.2, observed = location_obs)
        
        # Bernoulli deterministic variable with conditional probability 
        online = pm.Bernoulli('online', tt.switch(tt.eq(location, 1), tt.switch(tt.eq(gender, 1), 0.25, 0.45),
                                                                      tt.switch(tt.eq(gender, 1), 0.3, 0.6)),
                              observed = online_obs)
        
        # Bernoulli random variable with prior
        amount = pm.Bernoulli('amount', 0.1, observed = amount_obs)
        
        # Bernoulli deterministic variable with conditional probability
        fraudulent = pm.Bernoulli('fraudulent', tt.switch(tt.eq(amount, 1), tt.switch(tt.eq(online, 1), 0.45, 0.2),
                                                                            tt.switch(tt.eq(online, 1), 0.1, 0.05)),
                                  observed = fraudulent_obs)
        return model
def sample_model(model):
    """Sample from a pymc3 model"""
    
    with model:
        trace = pm.sample(draws=N_SAMPLES, chains=N_CHAINS)
    return trace
fraud_model = make_model()
fraud_trace = sample_model(fraud_model)
fraud_trace['fraudulent'].mean()
pm.traceplot(fraud_trace)
pm.forestplot(fraud_trace)
fraud_trace['gender'].mean()
fraud_trace['online'].mean()
fraud_model = make_model(amount_obs = 0)
fraud_trace = sample_model(fraud_model)
fraud_trace['fraudulent'].mean()
pm.traceplot(fraud_trace)
p_sanfran = fraud_trace['location']
p_sanfran.mean()
p_female = fraud_trace['gender']
p_female.mean()
p_fraud = fraud_trace['fraudulent']
p_fraud.mean()<jupyter_output><empty_output><jupyter_text>Bayesian Linear RegressionWe can also use Probabilistic Programming to perform Bayesian Linear RegressionIn ordinary least squares regression, the data generating process is:$$y = \beta * X + \epsilon$$In Bayesian Linear Regression, the data generating process is a distribution:$$ y \sim N(\mu, (\sigma) ^ 2)$$$$\mu = \alpha + \beta * X$$Where $\beta$ and X are matrices representing the weights and features respectively and $\alpha$ is the intercept.<jupyter_code>from pymc3 import GLM
data = pd.read_csv('../data/bike-sharing/day.csv', parse_dates = ['dteday']).set_index('dteday')
data['date'] = data.index
data.head()
plt.figure(figsize = (8, 6))
plt.plot(data['date'], data['cnt'], marker = '.');
data['day'] = (data.index - data.index[0]).days
data.head()

p = np.polyfit(data['day'], data['cnt'], 8)
y = np.polyval(p, data['day'])
plt.plot(y)
plt.xticks(data.index)
import seaborn as sns


sns.regplot('date', 'cnt', data, fit_reg = False);
formula = ('cnt ~ day + season + mnth + holiday + weekday +  workingday + weathersit + atemp + windspeed')
from pymc3 import GLM
with pm.Model() as bike_model:
    GLM.from_formula(formula, data)
    start = pm.find_MAP()
    trace = pm.sample(tune=1000, njobs=-1, start = start)
pm.traceplot(trace)
def lm(value, sample):
    prediction = sample['Intercept']
    for var in sample.keys():
        prediction += sample[var] * value
    return prediction
trace['season'] * data['season'][0]
len(trace['season'])
row.index
trace.varnames
estimates = []
days = []
for i, (day, row) in enumerate(data.iterrows()):
    estimate = trace['Intercept']
    for var in row.index:
        if var in trace.varnames:
            estimate += trace[var] * row[var]
    estimates.append(estimate)
    days.append(day for _ in range(len(trace[var])))
    print(f'{round(100 * (i / len(data)), 2)}% complete.', end = '\r')
len(days)
len(estimates)
est = pd.DataFrame(estimates)
est.index = data.index
est.head()
import matplotlib.pyplot as plt
%matplotlib inline
fig, ax = plt.subplots(figsize = (10, 10))
data['cnt'].plot(linewidth = 2, linestyle = '--', color = 'k', ax = ax);
est.plot(legend = None, ax = ax);
data['upper'] = np.percentile(est, 95, axis = 1)
data['lower'] = np.percentile(est, 5, axis = 1)

data.head()
plt.figure(figsize = (10, 10))
plt.fill_between(x = data.index, y1=data['lower'], alpha = 10, y2 = data['upper'],step = 'mid')
fig, ax = plt.subplots(figsize = (10, 10))
data['cnt'].plot(ax = ax, label = 'actual')
data['lower'].plot(ax = ax, label = 'lower')
data['upper'].plot(ax = ax, label = 'upper')
plt.legend();
estimate.shape
len(data)
pm.plots.plot_posterior(trace, varnames = ['cnt']);<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup