[Case Study]Bayesian Robust Regression in PyMC3

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Author: Thomas Wiecki

In this blog post I will write about:

• How a few outliers can largely affect the fit of linear regression models.
• How replacing the normal likelihood with Student T distribution produces robust regression.
• How this can easily be done with PyMC3 and its new glm module by passing a family object.
  

本帖隐藏的内容

Bayesian Robust Regression in PyMC3.pdf (439.23 KB)

2016-12-18 07:03:38 上传

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：regression Case study regressio Bayesian regress produces passing family affect easily

1. %matplotlib inline
   
2. 
   
3. import pymc3 as pm
   
4. 
   
5. import matplotlib.pyplot as plt
   
6. import numpy as np
   
7. 
   
8. import theano

复制代码

1. size = 100
   
2. true_intercept = 1
   
3. true_slope = 2
   
4. 
   
5. x = np.linspace(0, 1, size)
   
6. # y = a + b*x
   
7. true_regression_line = true_intercept + true_slope * x
   
8. # add noise
   
9. y = true_regression_line + np.random.normal(scale=.5, size=size)
   
10. 
    
11. # Add outliers
    
12. x_out = np.append(x, [.1, .15, .2])
    
13. y_out = np.append(y, [8, 6, 9])
    
14. 
    
15. data = dict(x=x_out, y=y_out)

复制代码

1. fig = plt.figure(figsize=(7, 7))
   
2. ax = fig.add_subplot(111, xlabel='x', ylabel='y', title='Generated data and underlying model')
   
3. ax.plot(x_out, y_out, 'x', label='sampled data')
   
4. ax.plot(x, true_regression_line, label='true regression line', lw=2.)
   
5. plt.legend(loc=0);

复制代码

1. with pm.Model() as model:
   
2.     pm.glm.glm('y ~ x', data)
   
3.     start = pm.find_MAP()
   
4.     step = pm.NUTS(scaling=start)
   
5.     trace = pm.sample(2000, step, progressbar=False)

复制代码

1. plt.subplot(111, xlabel='x', ylabel='y',
   
2.             title='Posterior predictive regression lines')
   
3. plt.plot(x_out, y_out, 'x', label='data')
   
4. pm.glm.plot_posterior_predictive(trace, samples=100,
   
5.                                  label='posterior predictive regression lines')
   
6. plt.plot(x, true_regression_line,
   
7.          label='true regression line', lw=3., c='y')
   
8. 
   
9. plt.legend(loc=0);

复制代码

1. Lets look at those two distributions to get a feel for them.
   
2. 
   
3. In [7]:
   
4. normal_dist = pm.Normal.dist(mu=0, sd=1)
   
5. t_dist = pm.StudentT.dist(mu=0, lam=1, nu=1)
   
6. x_eval = np.linspace(-8, 8, 300)
   
7. plt.plot(x_eval, theano.tensor.exp(normal_dist.logp(x_eval)).eval(), label='Normal', lw=2.)
   
8. plt.plot(x_eval, theano.tensor.exp(t_dist.logp(x_eval)).eval(), label='Student T', lw=2.)
   
9. plt.xlabel('x')
   
10. plt.ylabel('Probability density')
    
11. plt.legend();

复制代码

1. To define the usage of a T distribution in PyMC3 we can pass a family object -- T -- that specifies that our data is Student T-distributed (see glm.families for more choices). Note that this is the same syntax as R and statsmodels use.
   
2. 
   
3. In [8]:
   
4. with pm.Model() as model_robust:
   
5.     family = pm.glm.families.T()
   
6.     pm.glm.glm('y ~ x', data, family=family)
   
7.     start = pm.find_MAP()
   
8.     step = pm.NUTS(scaling=start)
   
9.     trace_robust = pm.sample(2000, step, progressbar=False)
   
10. 
    
11. plt.figure(figsize=(5, 5))
    
12. plt.plot(x_out, y_out, 'x')
    
13. pm.glm.plot_posterior_predictive(trace_robust,
    
14.                                  label='posterior predictive regression lines')
    
15. plt.plot(x, true_regression_line,
    
16.          label='true regression line', lw=3., c='y')
    
17. plt.legend()

复制代码