<jupyter_start><jupyter_text>T-Student regression---In this notebook we analyze how to adjust a *bayesian linear regression* code when the outcome has heavy tail and the hypotesis of *normal distribution* are not vadid. We use a synthetic dataset in order to see which is the impact of the right model on the results. The code used is *Python* with a specific bayesian framework **PyMC3**.<jupyter_code>import warnings
warnings.filterwarnings("ignore")
from scipy import stats
import arviz as az
import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import seaborn as sns
import pandas as pd
import random
from theano import shared
from sklearn import preprocessing
print('Running on PyMC3 v{}'.format(pm.__version__))<jupyter_output>Running on PyMC3 v3.10.0<jupyter_text>Let's simulate a riproducible synthetic dataset. By default we don't have mean and variance for *T-Student distribution* but we can extend to the usual three parameter **location-scale** family as follow:$X = \mu + \sigma T$, see [here](https://en.wikipedia.org/wiki/Student%27s_t-distributionGeneralized_Student's_t-distribution)<jupyter_code>random.seed(1990)
N  = 359   # Sample size
b0 = 3.8   # b0 coefficient (intercept)
b1 = 2.4   # b1 coefficient (slope)

# Simulate indipendent variable
x = np.linspace(-3,4,N)

# Simulate the output
y = b0 + b1*x + np.random.normal(0,3,N) + [*np.repeat([0], 300, axis=0), *20*np.random.binomial(1,0.3,59)]
plt.hist(y)
plt.title('Histogram Outcome')
plt.show()<jupyter_output><empty_output><jupyter_text>Gaussian Model<jupyter_code>with pm.Model() as model_gaussian:
    # Prior distribution
    sigma = pm.HalfCauchy("sigma", beta=10, testval=1.0)
    beta0 = pm.Normal("beta0", 0, sigma=20)
    beta1 = pm.Normal("beta1", 0, sigma=20)
    
    # Likelihood
    lik = pm.Normal('y', mu = beta0 + beta1 * x, sigma=sigma, observed=y)
        
    # draw 1000 posterior samples using NUTS sampling + 1000 of warmup on 4 chain
    trace_gaussian = pm.sample(1000, tune=1000)
az.plot_trace(trace_gaussian);<jupyter_output><empty_output><jupyter_text>T-Student Model<jupyter_code>with pm.Model() as model_tstudent:
    # Prior distribution
    sigma = pm.HalfCauchy("sigma", beta=10, testval=1.0)
    beta0 = pm.Normal("beta0", 0, sigma=20)
    beta1 = pm.Normal("beta1", 0, sigma=20)
    n = pm.Gamma('n', alpha = 2, beta = 0.1)
    
    # Likelihood
    lik = pm.StudentT('y', nu=n, mu = beta0 + beta1 * x, sigma=sigma, observed=y)
        
    # draw 1000 posterior samples using NUTS sampling + 1000 of warmup on 4 chain
    trace_tstudent = pm.sample(1000, tune=1000)    
az.plot_trace(trace_tstudent);<jupyter_output><empty_output><jupyter_text>Results<jupyter_code>az.summary(trace_gaussian)
az.summary(trace_tstudent)
fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(111, xlabel="x", ylabel="y", title="Generated data and underlying model")
ax.plot(x, y, "x", label="sampled data")
ax.plot(x, b0 + b1*x, label="true regression line", lw=2.0)
ax.plot(x, trace_gaussian["beta0"].mean() + trace_gaussian["beta1"].mean()*x, label="regression line", lw=2.0, color = "r")
ax.plot(x, trace_tstudent["beta0"].mean() + trace_tstudent["beta1"].mean()*x, label="t-regression line", lw=2.0, color = "m")
plt.legend(loc=0);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13235742604452577, 'fraction_numerical': 0.03964623360780726, 'mean_word_length': 5.007326007326007, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 20, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '496674', 'license': 'Apache-2.0', 'n_tokens_mistral': 1238, 'n_tokens_neox': 1142, 'n_words': 312, 'path': 't_regression.ipynb', 'repo_name': 'nicolofp/T-Student-Regression', 'size': '3279'}