M-Estimators for Robust Linear ModelingΒΆ
In [1]:
from __future__ import print_function
from statsmodels.compat import lmap
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import statsmodels.api as sm
- An M-estimator minimizes the function
$$Q(e_i, \rho) = \sum_i~\rho \left (\frac{e_i}{s}\right )$$
where $\rho$ is a symmetric function of the residuals
- The effect of $\rho$ is to reduce the influence of outliers
- $s$ is an estimate of scale.
- The robust estimates $\hat{\beta}$ are computed by the iteratively re-weighted least squares algorithm
- We have several choices available for the weighting functions to be used
In [2]:
norms = sm.robust.norms
In [3]:
def plot_weights(support, weights_func, xlabels, xticks):
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.plot(support, weights_func(support))
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels, fontsize=16)
ax.set_ylim(-.1, 1.1)
return ax
Andrew's Wave¶
In [4]:
help(norms.AndrewWave.weights)
Help on method weights in module statsmodels.robust.norms:
weights(self, z) unbound statsmodels.robust.norms.AndrewWave method
Andrew's wave weighting function for the IRLS algorithm
The psi function scaled by z
Parameters
----------
z : array-like
1d array
Returns
-------
weights : array
weights(z) = sin(z/a)/(z/a) for \|z\| <= a*pi
weights(z) = 0 for \|z\| > a*pi
In [5]:
a = 1.339
support = np.linspace(-np.pi*a, np.pi*a, 100)
andrew = norms.AndrewWave(a=a)
plot_weights(support, andrew.weights, ['$-\pi*a$', '0', '$\pi*a$'], [-np.pi*a, 0, np.pi*a]);
Error in callback <function post_execute at 0xb082b2b0> (for post_execute):
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) /usr/lib/python2.7/dist-packages/matplotlib/pyplot.pyc in post_execute() 145 def post_execute(): 146 if matplotlib.is_interactive(): --> 147 draw_all() 148 149 # IPython >= 2 /usr/lib/python2.7/dist-packages/matplotlib/_pylab_helpers.pyc in draw_all(cls, force) 148 for f_mgr in cls.get_all_fig_managers(): 149 if force or f_mgr.canvas.figure.stale: --> 150 f_mgr.canvas.draw_idle() 151 152 atexit.register(Gcf.destroy_all) /usr/lib/python2.7/dist-packages/matplotlib/backend_bases.pyc in draw_idle(self, *args, **kwargs) 2024 if not self._is_idle_drawing: 2025 with self._idle_draw_cntx(): -> 2026 self.draw(*args, **kwargs) 2027 2028 def draw_cursor(self, event): /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in draw(self) 472 473 try: --> 474 self.figure.draw(self.renderer) 475 finally: 476 RendererAgg.lock.release() /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/figure.pyc in draw(self, renderer) 1157 dsu.sort(key=itemgetter(0)) 1158 for zorder, a, func, args in dsu: -> 1159 func(*args) 1160 1161 renderer.close_group('figure') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in draw(self, renderer, inframe) 2322 2323 for zorder, a in dsu: -> 2324 a.draw(renderer) 2325 2326 renderer.close_group('axes') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in draw(self, renderer, *args, **kwargs) 1106 ticks_to_draw = self._update_ticks(renderer) 1107 ticklabelBoxes, ticklabelBoxes2 = self._get_tick_bboxes(ticks_to_draw, -> 1108 renderer) 1109 1110 for tick in ticks_to_draw: /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in _get_tick_bboxes(self, ticks, renderer) 1056 for tick in ticks: 1057 if tick.label1On and tick.label1.get_visible(): -> 1058 extent = tick.label1.get_window_extent(renderer) 1059 ticklabelBoxes.append(extent) 1060 if tick.label2On and tick.label2.get_visible(): /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in get_window_extent(self, renderer, dpi) 959 raise RuntimeError('Cannot get window extent w/o renderer') 960 --> 961 bbox, info, descent = self._get_layout(self._renderer) 962 x, y = self.get_unitless_position() 963 x, y = self.get_transform().transform_point((x, y)) /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in _get_layout(self, renderer) 350 tmp, lp_h, lp_bl = renderer.get_text_width_height_descent('lp', 351 self._fontproperties, --> 352 ismath=False) 353 offsety = (lp_h - lp_bl) * self._linespacing 354 /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in get_text_width_height_descent(self, s, prop, ismath) 227 fontsize = prop.get_size_in_points() 228 w, h, d = texmanager.get_text_width_height_descent(s, fontsize, --> 229 renderer=self) 230 return w, h, d 231 /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in get_text_width_height_descent(self, tex, fontsize, renderer) 673 else: 674 # use dviread. It sometimes returns a wrong descent. --> 675 dvifile = self.make_dvi(tex, fontsize) 676 dvi = dviread.Dvi(dvifile, 72 * dpi_fraction) 677 try: /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in make_dvi(self, tex, fontsize) 420 'string:\n%s\nHere is the full report generated by ' 421 'LaTeX: \n\n' % repr(tex.encode('unicode_escape')) + --> 422 report)) 423 else: 424 mpl.verbose.report(report, 'debug') RuntimeError: LaTeX was not able to process the following string: 'lp' Here is the full report generated by LaTeX:
<matplotlib.figure.Figure at 0xa7237050>
Hampel's 17A¶
In [6]:
help(norms.Hampel.weights)
Help on method weights in module statsmodels.robust.norms:
weights(self, z) unbound statsmodels.robust.norms.Hampel method
Hampel weighting function for the IRLS algorithm
The psi function scaled by z
Parameters
----------
z : array-like
1d array
Returns
-------
weights : array
weights(z) = 1 for \|z\| <= a
weights(z) = a/\|z\| for a < \|z\| <= b
weights(z) = a*(c - \|z\|)/(\|z\|*(c-b)) for b < \|z\| <= c
weights(z) = 0 for \|z\| > c
In [7]:
c = 8
support = np.linspace(-3*c, 3*c, 1000)
hampel = norms.Hampel(a=2., b=4., c=c)
plot_weights(support, hampel.weights, ['3*c', '0', '3*c'], [-3*c, 0, 3*c]);
Error in callback <function post_execute at 0xb082b2b0> (for post_execute):
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) /usr/lib/python2.7/dist-packages/matplotlib/pyplot.pyc in post_execute() 145 def post_execute(): 146 if matplotlib.is_interactive(): --> 147 draw_all() 148 149 # IPython >= 2 /usr/lib/python2.7/dist-packages/matplotlib/_pylab_helpers.pyc in draw_all(cls, force) 148 for f_mgr in cls.get_all_fig_managers(): 149 if force or f_mgr.canvas.figure.stale: --> 150 f_mgr.canvas.draw_idle() 151 152 atexit.register(Gcf.destroy_all) /usr/lib/python2.7/dist-packages/matplotlib/backend_bases.pyc in draw_idle(self, *args, **kwargs) 2024 if not self._is_idle_drawing: 2025 with self._idle_draw_cntx(): -> 2026 self.draw(*args, **kwargs) 2027 2028 def draw_cursor(self, event): /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in draw(self) 472 473 try: --> 474 self.figure.draw(self.renderer) 475 finally: 476 RendererAgg.lock.release() /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/figure.pyc in draw(self, renderer) 1157 dsu.sort(key=itemgetter(0)) 1158 for zorder, a, func, args in dsu: -> 1159 func(*args) 1160 1161 renderer.close_group('figure') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in draw(self, renderer, inframe) 2322 2323 for zorder, a in dsu: -> 2324 a.draw(renderer) 2325 2326 renderer.close_group('axes') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in draw(self, renderer, *args, **kwargs) 1106 ticks_to_draw = self._update_ticks(renderer) 1107 ticklabelBoxes, ticklabelBoxes2 = self._get_tick_bboxes(ticks_to_draw, -> 1108 renderer) 1109 1110 for tick in ticks_to_draw: /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in _get_tick_bboxes(self, ticks, renderer) 1056 for tick in ticks: 1057 if tick.label1On and tick.label1.get_visible(): -> 1058 extent = tick.label1.get_window_extent(renderer) 1059 ticklabelBoxes.append(extent) 1060 if tick.label2On and tick.label2.get_visible(): /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in get_window_extent(self, renderer, dpi) 959 raise RuntimeError('Cannot get window extent w/o renderer') 960 --> 961 bbox, info, descent = self._get_layout(self._renderer) 962 x, y = self.get_unitless_position() 963 x, y = self.get_transform().transform_point((x, y)) /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in _get_layout(self, renderer) 350 tmp, lp_h, lp_bl = renderer.get_text_width_height_descent('lp', 351 self._fontproperties, --> 352 ismath=False) 353 offsety = (lp_h - lp_bl) * self._linespacing 354 /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in get_text_width_height_descent(self, s, prop, ismath) 227 fontsize = prop.get_size_in_points() 228 w, h, d = texmanager.get_text_width_height_descent(s, fontsize, --> 229 renderer=self) 230 return w, h, d 231 /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in get_text_width_height_descent(self, tex, fontsize, renderer) 673 else: 674 # use dviread. It sometimes returns a wrong descent. --> 675 dvifile = self.make_dvi(tex, fontsize) 676 dvi = dviread.Dvi(dvifile, 72 * dpi_fraction) 677 try: /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in make_dvi(self, tex, fontsize) 420 'string:\n%s\nHere is the full report generated by ' 421 'LaTeX: \n\n' % repr(tex.encode('unicode_escape')) + --> 422 report)) 423 else: 424 mpl.verbose.report(report, 'debug') RuntimeError: LaTeX was not able to process the following string: 'lp' Here is the full report generated by LaTeX:
<matplotlib.figure.Figure at 0xa746da30>
Huber's t¶
In [8]:
help(norms.HuberT.weights)
Help on method weights in module statsmodels.robust.norms:
weights(self, z) unbound statsmodels.robust.norms.HuberT method
Huber's t weighting function for the IRLS algorithm
The psi function scaled by z
Parameters
----------
z : array-like
1d array
Returns
-------
weights : array
weights(z) = 1 for \|z\| <= t
weights(z) = t/\|z\| for \|z\| > t
In [9]:
t = 1.345
support = np.linspace(-3*t, 3*t, 1000)
huber = norms.HuberT(t=t)
plot_weights(support, huber.weights, ['-3*t', '0', '3*t'], [-3*t, 0, 3*t]);
Error in callback <function post_execute at 0xb082b2b0> (for post_execute):
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) /usr/lib/python2.7/dist-packages/matplotlib/pyplot.pyc in post_execute() 145 def post_execute(): 146 if matplotlib.is_interactive(): --> 147 draw_all() 148 149 # IPython >= 2 /usr/lib/python2.7/dist-packages/matplotlib/_pylab_helpers.pyc in draw_all(cls, force) 148 for f_mgr in cls.get_all_fig_managers(): 149 if force or f_mgr.canvas.figure.stale: --> 150 f_mgr.canvas.draw_idle() 151 152 atexit.register(Gcf.destroy_all) /usr/lib/python2.7/dist-packages/matplotlib/backend_bases.pyc in draw_idle(self, *args, **kwargs) 2024 if not self._is_idle_drawing: 2025 with self._idle_draw_cntx(): -> 2026 self.draw(*args, **kwargs) 2027 2028 def draw_cursor(self, event): /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in draw(self) 472 473 try: --> 474 self.figure.draw(self.renderer) 475 finally: 476 RendererAgg.lock.release() /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/figure.pyc in draw(self, renderer) 1157 dsu.sort(key=itemgetter(0)) 1158 for zorder, a, func, args in dsu: -> 1159 func(*args) 1160 1161 renderer.close_group('figure') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in draw(self, renderer, inframe) 2322 2323 for zorder, a in dsu: -> 2324 a.draw(renderer) 2325 2326 renderer.close_group('axes') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in draw(self, renderer, *args, **kwargs) 1106 ticks_to_draw = self._update_ticks(renderer) 1107 ticklabelBoxes, ticklabelBoxes2 = self._get_tick_bboxes(ticks_to_draw, -> 1108 renderer) 1109 1110 for tick in ticks_to_draw: /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in _get_tick_bboxes(self, ticks, renderer) 1056 for tick in ticks: 1057 if tick.label1On and tick.label1.get_visible(): -> 1058 extent = tick.label1.get_window_extent(renderer) 1059 ticklabelBoxes.append(extent) 1060 if tick.label2On and tick.label2.get_visible(): /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in get_window_extent(self, renderer, dpi) 959 raise RuntimeError('Cannot get window extent w/o renderer') 960 --> 961 bbox, info, descent = self._get_layout(self._renderer) 962 x, y = self.get_unitless_position() 963 x, y = self.get_transform().transform_point((x, y)) /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in _get_layout(self, renderer) 350 tmp, lp_h, lp_bl = renderer.get_text_width_height_descent('lp', 351 self._fontproperties, --> 352 ismath=False) 353 offsety = (lp_h - lp_bl) * self._linespacing 354 /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in get_text_width_height_descent(self, s, prop, ismath) 227 fontsize = prop.get_size_in_points() 228 w, h, d = texmanager.get_text_width_height_descent(s, fontsize, --> 229 renderer=self) 230 return w, h, d 231 /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in get_text_width_height_descent(self, tex, fontsize, renderer) 673 else: 674 # use dviread. It sometimes returns a wrong descent. --> 675 dvifile = self.make_dvi(tex, fontsize) 676 dvi = dviread.Dvi(dvifile, 72 * dpi_fraction) 677 try: /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in make_dvi(self, tex, fontsize) 420 'string:\n%s\nHere is the full report generated by ' 421 'LaTeX: \n\n' % repr(tex.encode('unicode_escape')) + --> 422 report)) 423 else: 424 mpl.verbose.report(report, 'debug') RuntimeError: LaTeX was not able to process the following string: 'lp' Here is the full report generated by LaTeX:
<matplotlib.figure.Figure at 0xa721b970>
Least Squares¶
In [10]:
help(norms.LeastSquares.weights)
Help on method weights in module statsmodels.robust.norms:
weights(self, z) unbound statsmodels.robust.norms.LeastSquares method
The least squares estimator weighting function for the IRLS algorithm.
The psi function scaled by the input z
Parameters
----------
z : array-like
1d array
Returns
-------
weights : array
weights(z) = np.ones(z.shape)
In [11]:
support = np.linspace(-3, 3, 1000)
lst_sq = norms.LeastSquares()
plot_weights(support, lst_sq.weights, ['-3', '0', '3'], [-3, 0, 3]);
Error in callback <function post_execute at 0xb082b2b0> (for post_execute):
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) /usr/lib/python2.7/dist-packages/matplotlib/pyplot.pyc in post_execute() 145 def post_execute(): 146 if matplotlib.is_interactive(): --> 147 draw_all() 148 149 # IPython >= 2 /usr/lib/python2.7/dist-packages/matplotlib/_pylab_helpers.pyc in draw_all(cls, force) 148 for f_mgr in cls.get_all_fig_managers(): 149 if force or f_mgr.canvas.figure.stale: --> 150 f_mgr.canvas.draw_idle() 151 152 atexit.register(Gcf.destroy_all) /usr/lib/python2.7/dist-packages/matplotlib/backend_bases.pyc in draw_idle(self, *args, **kwargs) 2024 if not self._is_idle_drawing: 2025 with self._idle_draw_cntx(): -> 2026 self.draw(*args, **kwargs) 2027 2028 def draw_cursor(self, event): /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in draw(self) 472 473 try: --> 474 self.figure.draw(self.renderer) 475 finally: 476 RendererAgg.lock.release() /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/figure.pyc in draw(self, renderer) 1157 dsu.sort(key=itemgetter(0)) 1158 for zorder, a, func, args in dsu: -> 1159 func(*args) 1160 1161 renderer.close_group('figure') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in draw(self, renderer, inframe) 2322 2323 for zorder, a in dsu: -> 2324 a.draw(renderer) 2325 2326 renderer.close_group('axes') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in draw(self, renderer, *args, **kwargs) 1106 ticks_to_draw = self._update_ticks(renderer) 1107 ticklabelBoxes, ticklabelBoxes2 = self._get_tick_bboxes(ticks_to_draw, -> 1108 renderer) 1109 1110 for tick in ticks_to_draw: /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in _get_tick_bboxes(self, ticks, renderer) 1056 for tick in ticks: 1057 if tick.label1On and tick.label1.get_visible(): -> 1058 extent = tick.label1.get_window_extent(renderer) 1059 ticklabelBoxes.append(extent) 1060 if tick.label2On and tick.label2.get_visible(): /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in get_window_extent(self, renderer, dpi) 959 raise RuntimeError('Cannot get window extent w/o renderer') 960 --> 961 bbox, info, descent = self._get_layout(self._renderer) 962 x, y = self.get_unitless_position() 963 x, y = self.get_transform().transform_point((x, y)) /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in _get_layout(self, renderer) 350 tmp, lp_h, lp_bl = renderer.get_text_width_height_descent('lp', 351 self._fontproperties, --> 352 ismath=False) 353 offsety = (lp_h - lp_bl) * self._linespacing 354 /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in get_text_width_height_descent(self, s, prop, ismath) 227 fontsize = prop.get_size_in_points() 228 w, h, d = texmanager.get_text_width_height_descent(s, fontsize, --> 229 renderer=self) 230 return w, h, d 231 /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in get_text_width_height_descent(self, tex, fontsize, renderer) 673 else: 674 # use dviread. It sometimes returns a wrong descent. --> 675 dvifile = self.make_dvi(tex, fontsize) 676 dvi = dviread.Dvi(dvifile, 72 * dpi_fraction) 677 try: /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in make_dvi(self, tex, fontsize) 420 'string:\n%s\nHere is the full report generated by ' 421 'LaTeX: \n\n' % repr(tex.encode('unicode_escape')) + --> 422 report)) 423 else: 424 mpl.verbose.report(report, 'debug') RuntimeError: LaTeX was not able to process the following string: 'lp' Here is the full report generated by LaTeX:
<matplotlib.figure.Figure at 0xa7378d70>
Ramsay's Ea¶
In [12]:
help(norms.RamsayE.weights)
Help on method weights in module statsmodels.robust.norms:
weights(self, z) unbound statsmodels.robust.norms.RamsayE method
Ramsay's Ea weighting function for the IRLS algorithm
The psi function scaled by z
Parameters
----------
z : array-like
1d array
Returns
-------
weights : array
weights(z) = exp(-a*\|z\|)
In [13]:
a = .3
support = np.linspace(-3*a, 3*a, 1000)
ramsay = norms.RamsayE(a=a)
plot_weights(support, ramsay.weights, ['-3*a', '0', '3*a'], [-3*a, 0, 3*a]);
Error in callback <function post_execute at 0xb082b2b0> (for post_execute):
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) /usr/lib/python2.7/dist-packages/matplotlib/pyplot.pyc in post_execute() 145 def post_execute(): 146 if matplotlib.is_interactive(): --> 147 draw_all() 148 149 # IPython >= 2 /usr/lib/python2.7/dist-packages/matplotlib/_pylab_helpers.pyc in draw_all(cls, force) 148 for f_mgr in cls.get_all_fig_managers(): 149 if force or f_mgr.canvas.figure.stale: --> 150 f_mgr.canvas.draw_idle() 151 152 atexit.register(Gcf.destroy_all) /usr/lib/python2.7/dist-packages/matplotlib/backend_bases.pyc in draw_idle(self, *args, **kwargs) 2024 if not self._is_idle_drawing: 2025 with self._idle_draw_cntx(): -> 2026 self.draw(*args, **kwargs) 2027 2028 def draw_cursor(self, event): /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in draw(self) 472 473 try: --> 474 self.figure.draw(self.renderer) 475 finally: 476 RendererAgg.lock.release() /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/figure.pyc in draw(self, renderer) 1157 dsu.sort(key=itemgetter(0)) 1158 for zorder, a, func, args in dsu: -> 1159 func(*args) 1160 1161 renderer.close_group('figure') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in draw(self, renderer, inframe) 2322 2323 for zorder, a in dsu: -> 2324 a.draw(renderer) 2325 2326 renderer.close_group('axes') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in draw(self, renderer, *args, **kwargs) 1106 ticks_to_draw = self._update_ticks(renderer) 1107 ticklabelBoxes, ticklabelBoxes2 = self._get_tick_bboxes(ticks_to_draw, -> 1108 renderer) 1109 1110 for tick in ticks_to_draw: /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in _get_tick_bboxes(self, ticks, renderer) 1056 for tick in ticks: 1057 if tick.label1On and tick.label1.get_visible(): -> 1058 extent = tick.label1.get_window_extent(renderer) 1059 ticklabelBoxes.append(extent) 1060 if tick.label2On and tick.label2.get_visible(): /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in get_window_extent(self, renderer, dpi) 959 raise RuntimeError('Cannot get window extent w/o renderer') 960 --> 961 bbox, info, descent = self._get_layout(self._renderer) 962 x, y = self.get_unitless_position() 963 x, y = self.get_transform().transform_point((x, y)) /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in _get_layout(self, renderer) 350 tmp, lp_h, lp_bl = renderer.get_text_width_height_descent('lp', 351 self._fontproperties, --> 352 ismath=False) 353 offsety = (lp_h - lp_bl) * self._linespacing 354 /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in get_text_width_height_descent(self, s, prop, ismath) 227 fontsize = prop.get_size_in_points() 228 w, h, d = texmanager.get_text_width_height_descent(s, fontsize, --> 229 renderer=self) 230 return w, h, d 231 /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in get_text_width_height_descent(self, tex, fontsize, renderer) 673 else: 674 # use dviread. It sometimes returns a wrong descent. --> 675 dvifile = self.make_dvi(tex, fontsize) 676 dvi = dviread.Dvi(dvifile, 72 * dpi_fraction) 677 try: /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in make_dvi(self, tex, fontsize) 420 'string:\n%s\nHere is the full report generated by ' 421 'LaTeX: \n\n' % repr(tex.encode('unicode_escape')) + --> 422 report)) 423 else: 424 mpl.verbose.report(report, 'debug') RuntimeError: LaTeX was not able to process the following string: 'lp' Here is the full report generated by LaTeX:
<matplotlib.figure.Figure at 0xa72b76d0>
Trimmed Mean¶
In [14]:
help(norms.TrimmedMean.weights)
Help on method weights in module statsmodels.robust.norms:
weights(self, z) unbound statsmodels.robust.norms.TrimmedMean method
Least trimmed mean weighting function for the IRLS algorithm
The psi function scaled by z
Parameters
----------
z : array-like
1d array
Returns
-------
weights : array
weights(z) = 1 for \|z\| <= c
weights(z) = 0 for \|z\| > c
In [15]:
c = 2
support = np.linspace(-3*c, 3*c, 1000)
trimmed = norms.TrimmedMean(c=c)
plot_weights(support, trimmed.weights, ['-3*c', '0', '3*c'], [-3*c, 0, 3*c]);
Error in callback <function post_execute at 0xb082b2b0> (for post_execute):
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) /usr/lib/python2.7/dist-packages/matplotlib/pyplot.pyc in post_execute() 145 def post_execute(): 146 if matplotlib.is_interactive(): --> 147 draw_all() 148 149 # IPython >= 2 /usr/lib/python2.7/dist-packages/matplotlib/_pylab_helpers.pyc in draw_all(cls, force) 148 for f_mgr in cls.get_all_fig_managers(): 149 if force or f_mgr.canvas.figure.stale: --> 150 f_mgr.canvas.draw_idle() 151 152 atexit.register(Gcf.destroy_all) /usr/lib/python2.7/dist-packages/matplotlib/backend_bases.pyc in draw_idle(self, *args, **kwargs) 2024 if not self._is_idle_drawing: 2025 with self._idle_draw_cntx(): -> 2026 self.draw(*args, **kwargs) 2027 2028 def draw_cursor(self, event): /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in draw(self) 472 473 try: --> 474 self.figure.draw(self.renderer) 475 finally: 476 RendererAgg.lock.release() /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/figure.pyc in draw(self, renderer) 1157 dsu.sort(key=itemgetter(0)) 1158 for zorder, a, func, args in dsu: -> 1159 func(*args) 1160 1161 renderer.close_group('figure') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in draw(self, renderer, inframe) 2322 2323 for zorder, a in dsu: -> 2324 a.draw(renderer) 2325 2326 renderer.close_group('axes') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in draw(self, renderer, *args, **kwargs) 1106 ticks_to_draw = self._update_ticks(renderer) 1107 ticklabelBoxes, ticklabelBoxes2 = self._get_tick_bboxes(ticks_to_draw, -> 1108 renderer) 1109 1110 for tick in ticks_to_draw: /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in _get_tick_bboxes(self, ticks, renderer) 1056 for tick in ticks: 1057 if tick.label1On and tick.label1.get_visible(): -> 1058 extent = tick.label1.get_window_extent(renderer) 1059 ticklabelBoxes.append(extent) 1060 if tick.label2On and tick.label2.get_visible(): /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in get_window_extent(self, renderer, dpi) 959 raise RuntimeError('Cannot get window extent w/o renderer') 960 --> 961 bbox, info, descent = self._get_layout(self._renderer) 962 x, y = self.get_unitless_position() 963 x, y = self.get_transform().transform_point((x, y)) /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in _get_layout(self, renderer) 350 tmp, lp_h, lp_bl = renderer.get_text_width_height_descent('lp', 351 self._fontproperties, --> 352 ismath=False) 353 offsety = (lp_h - lp_bl) * self._linespacing 354 /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in get_text_width_height_descent(self, s, prop, ismath) 227 fontsize = prop.get_size_in_points() 228 w, h, d = texmanager.get_text_width_height_descent(s, fontsize, --> 229 renderer=self) 230 return w, h, d 231 /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in get_text_width_height_descent(self, tex, fontsize, renderer) 673 else: 674 # use dviread. It sometimes returns a wrong descent. --> 675 dvifile = self.make_dvi(tex, fontsize) 676 dvi = dviread.Dvi(dvifile, 72 * dpi_fraction) 677 try: /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in make_dvi(self, tex, fontsize) 420 'string:\n%s\nHere is the full report generated by ' 421 'LaTeX: \n\n' % repr(tex.encode('unicode_escape')) + --> 422 report)) 423 else: 424 mpl.verbose.report(report, 'debug') RuntimeError: LaTeX was not able to process the following string: 'lp' Here is the full report generated by LaTeX:
<matplotlib.figure.Figure at 0xa73f6650>
Tukey's Biweight¶
In [16]:
help(norms.TukeyBiweight.weights)
Help on method weights in module statsmodels.robust.norms:
weights(self, z) unbound statsmodels.robust.norms.TukeyBiweight method
Tukey's biweight weighting function for the IRLS algorithm
The psi function scaled by z
Parameters
----------
z : array-like
1d array
Returns
-------
weights : array
psi(z) = (1 - (z/c)**2)**2 for \|z\| <= R
psi(z) = 0 for \|z\| > R
In [17]:
c = 4.685
support = np.linspace(-3*c, 3*c, 1000)
tukey = norms.TukeyBiweight(c=c)
plot_weights(support, tukey.weights, ['-3*c', '0', '3*c'], [-3*c, 0, 3*c]);
Error in callback <function post_execute at 0xb082b2b0> (for post_execute):
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) /usr/lib/python2.7/dist-packages/matplotlib/pyplot.pyc in post_execute() 145 def post_execute(): 146 if matplotlib.is_interactive(): --> 147 draw_all() 148 149 # IPython >= 2 /usr/lib/python2.7/dist-packages/matplotlib/_pylab_helpers.pyc in draw_all(cls, force) 148 for f_mgr in cls.get_all_fig_managers(): 149 if force or f_mgr.canvas.figure.stale: --> 150 f_mgr.canvas.draw_idle() 151 152 atexit.register(Gcf.destroy_all) /usr/lib/python2.7/dist-packages/matplotlib/backend_bases.pyc in draw_idle(self, *args, **kwargs) 2024 if not self._is_idle_drawing: 2025 with self._idle_draw_cntx(): -> 2026 self.draw(*args, **kwargs) 2027 2028 def draw_cursor(self, event): /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in draw(self) 472 473 try: --> 474 self.figure.draw(self.renderer) 475 finally: 476 RendererAgg.lock.release() /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/figure.pyc in draw(self, renderer) 1157 dsu.sort(key=itemgetter(0)) 1158 for zorder, a, func, args in dsu: -> 1159 func(*args) 1160 1161 renderer.close_group('figure') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in draw(self, renderer, inframe) 2322 2323 for zorder, a in dsu: -> 2324 a.draw(renderer) 2325 2326 renderer.close_group('axes') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in draw(self, renderer, *args, **kwargs) 1106 ticks_to_draw = self._update_ticks(renderer) 1107 ticklabelBoxes, ticklabelBoxes2 = self._get_tick_bboxes(ticks_to_draw, -> 1108 renderer) 1109 1110 for tick in ticks_to_draw: /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in _get_tick_bboxes(self, ticks, renderer) 1056 for tick in ticks: 1057 if tick.label1On and tick.label1.get_visible(): -> 1058 extent = tick.label1.get_window_extent(renderer) 1059 ticklabelBoxes.append(extent) 1060 if tick.label2On and tick.label2.get_visible(): /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in get_window_extent(self, renderer, dpi) 959 raise RuntimeError('Cannot get window extent w/o renderer') 960 --> 961 bbox, info, descent = self._get_layout(self._renderer) 962 x, y = self.get_unitless_position() 963 x, y = self.get_transform().transform_point((x, y)) /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in _get_layout(self, renderer) 350 tmp, lp_h, lp_bl = renderer.get_text_width_height_descent('lp', 351 self._fontproperties, --> 352 ismath=False) 353 offsety = (lp_h - lp_bl) * self._linespacing 354 /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in get_text_width_height_descent(self, s, prop, ismath) 227 fontsize = prop.get_size_in_points() 228 w, h, d = texmanager.get_text_width_height_descent(s, fontsize, --> 229 renderer=self) 230 return w, h, d 231 /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in get_text_width_height_descent(self, tex, fontsize, renderer) 673 else: 674 # use dviread. It sometimes returns a wrong descent. --> 675 dvifile = self.make_dvi(tex, fontsize) 676 dvi = dviread.Dvi(dvifile, 72 * dpi_fraction) 677 try: /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in make_dvi(self, tex, fontsize) 420 'string:\n%s\nHere is the full report generated by ' 421 'LaTeX: \n\n' % repr(tex.encode('unicode_escape')) + --> 422 report)) 423 else: 424 mpl.verbose.report(report, 'debug') RuntimeError: LaTeX was not able to process the following string: 'lp' Here is the full report generated by LaTeX:
<matplotlib.figure.Figure at 0xa7388e10>
Scale Estimators¶
- Robust estimates of the location
In [18]:
x = np.array([1, 2, 3, 4, 500])
- The mean is not a robust estimator of location
In [19]:
x.mean()
Out[19]:
102.0
- The median, on the other hand, is a robust estimator with a breakdown point of 50%
In [20]:
np.median(x)
Out[20]:
3.0
- Analagously for the scale
- The standard deviation is not robust
In [21]:
x.std()
Out[21]:
199.00251254695254
Median Absolute Deviation
$$ median_i |X_i - median_j(X_j)|) $$
Standardized Median Absolute Deviation is a consistent estimator for $\hat{\sigma}$
$$\hat{\sigma}=K \cdot MAD$$
where $K$ depends on the distribution. For the normal distribution for example,
$$K = \Phi^{-1}(.75)$$
In [22]:
stats.norm.ppf(.75)
Out[22]:
0.67448975019608171
In [23]:
print(x)
[ 1 2 3 4 500]
In [24]:
sm.robust.scale.stand_mad(x)
/statsmodels-0.6.1/debian/python-statsmodels/usr/lib/python2.7/dist-packages/statsmodels/robust/scale.py:49: FutureWarning: stand_mad is deprecated and will be removed in 0.7.0. Use mad instead. "instead.", FutureWarning)
Out[24]:
1.482602218505602
In [25]:
np.array([1,2,3,4,5.]).std()
Out[25]:
1.4142135623730951
- The default for Robust Linear Models is MAD
- another popular choice is Huber's proposal 2
In [26]:
np.random.seed(12345)
fat_tails = stats.t(6).rvs(40)
In [27]:
kde = sm.nonparametric.KDE(fat_tails)
kde.fit()
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.plot(kde.support, kde.density);
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) <ipython-input-240-605f8dff1422> in <module>() ----> 1 kde = sm.nonparametric.KDE(fat_tails) 2 kde.fit() 3 fig = plt.figure(figsize=(12,8)) 4 ax = fig.add_subplot(111) 5 ax.plot(kde.support, kde.density); AttributeError: 'module' object has no attribute 'KDE'
In [28]:
print(fat_tails.mean(), fat_tails.std())
0.0688231044811 1.34716332297
In [29]:
print(stats.norm.fit(fat_tails))
(0.068823104481087499, 1.3471633229698652)
In [30]:
print(stats.t.fit(fat_tails, f0=6))
(6, 0.039009187170278181, 1.0564230978488927)
In [31]:
huber = sm.robust.scale.Huber()
loc, scale = huber(fat_tails)
print(loc, scale)
0.0404898433327 1.15571400476
In [32]:
sm.robust.stand_mad(fat_tails)
Out[32]:
1.1153350011654151
In [33]:
sm.robust.stand_mad(fat_tails, c=stats.t(6).ppf(.75))
Out[33]:
1.0483916565928972
In [34]:
sm.robust.scale.mad(fat_tails)
Out[34]:
1.1153350011654151
Duncan's Occupational Prestige data - M-estimation for outliers¶
In [35]:
from statsmodels.graphics.api import abline_plot
from statsmodels.formula.api import ols, rlm
In [36]:
prestige = sm.datasets.get_rdataset("Duncan", "car", cache=True).data
In [37]:
print(prestige.head(10))
type income education prestige
accountant prof 62 86 82
pilot prof 72 76 83
architect prof 75 92 90
author prof 55 90 76
chemist prof 64 86 90
minister prof 21 84 87
professor prof 64 93 93
dentist prof 80 100 90
reporter wc 67 87 52
engineer prof 72 86 88
In [38]:
fig = plt.figure(figsize=(12,12))
ax1 = fig.add_subplot(211, xlabel='Income', ylabel='Prestige')
ax1.scatter(prestige.income, prestige.prestige)
xy_outlier = prestige.ix['minister'][['income','prestige']]
ax1.annotate('Minister', xy_outlier, xy_outlier+1, fontsize=16)
ax2 = fig.add_subplot(212, xlabel='Education',
ylabel='Prestige')
ax2.scatter(prestige.education, prestige.prestige);
Error in callback <function post_execute at 0xb082b2b0> (for post_execute):
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) /usr/lib/python2.7/dist-packages/matplotlib/pyplot.pyc in post_execute() 145 def post_execute(): 146 if matplotlib.is_interactive(): --> 147 draw_all() 148 149 # IPython >= 2 /usr/lib/python2.7/dist-packages/matplotlib/_pylab_helpers.pyc in draw_all(cls, force) 148 for f_mgr in cls.get_all_fig_managers(): 149 if force or f_mgr.canvas.figure.stale: --> 150 f_mgr.canvas.draw_idle() 151 152 atexit.register(Gcf.destroy_all) /usr/lib/python2.7/dist-packages/matplotlib/backend_bases.pyc in draw_idle(self, *args, **kwargs) 2024 if not self._is_idle_drawing: 2025 with self._idle_draw_cntx(): -> 2026 self.draw(*args, **kwargs) 2027 2028 def draw_cursor(self, event): /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in draw(self) 472 473 try: --> 474 self.figure.draw(self.renderer) 475 finally: 476 RendererAgg.lock.release() /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/figure.pyc in draw(self, renderer) 1157 dsu.sort(key=itemgetter(0)) 1158 for zorder, a, func, args in dsu: -> 1159 func(*args) 1160 1161 renderer.close_group('figure') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in draw(self, renderer, inframe) 2322 2323 for zorder, a in dsu: -> 2324 a.draw(renderer) 2325 2326 renderer.close_group('axes') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in draw(self, renderer, *args, **kwargs) 1106 ticks_to_draw = self._update_ticks(renderer) 1107 ticklabelBoxes, ticklabelBoxes2 = self._get_tick_bboxes(ticks_to_draw, -> 1108 renderer) 1109 1110 for tick in ticks_to_draw: /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in _get_tick_bboxes(self, ticks, renderer) 1056 for tick in ticks: 1057 if tick.label1On and tick.label1.get_visible(): -> 1058 extent = tick.label1.get_window_extent(renderer) 1059 ticklabelBoxes.append(extent) 1060 if tick.label2On and tick.label2.get_visible(): /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in get_window_extent(self, renderer, dpi) 959 raise RuntimeError('Cannot get window extent w/o renderer') 960 --> 961 bbox, info, descent = self._get_layout(self._renderer) 962 x, y = self.get_unitless_position() 963 x, y = self.get_transform().transform_point((x, y)) /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in _get_layout(self, renderer) 350 tmp, lp_h, lp_bl = renderer.get_text_width_height_descent('lp', 351 self._fontproperties, --> 352 ismath=False) 353 offsety = (lp_h - lp_bl) * self._linespacing 354 /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in get_text_width_height_descent(self, s, prop, ismath) 227 fontsize = prop.get_size_in_points() 228 w, h, d = texmanager.get_text_width_height_descent(s, fontsize, --> 229 renderer=self) 230 return w, h, d 231 /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in get_text_width_height_descent(self, tex, fontsize, renderer) 673 else: 674 # use dviread. It sometimes returns a wrong descent. --> 675 dvifile = self.make_dvi(tex, fontsize) 676 dvi = dviread.Dvi(dvifile, 72 * dpi_fraction) 677 try: /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in make_dvi(self, tex, fontsize) 420 'string:\n%s\nHere is the full report generated by ' 421 'LaTeX: \n\n' % repr(tex.encode('unicode_escape')) + --> 422 report)) 423 else: 424 mpl.verbose.report(report, 'debug') RuntimeError: LaTeX was not able to process the following string: 'lp' Here is the full report generated by LaTeX:
<matplotlib.figure.Figure at 0xa719f3d0>
In [39]:
ols_model = ols('prestige ~ income + education', prestige).fit()
print(ols_model.summary())
OLS Regression Results
==============================================================================
Dep. Variable: prestige R-squared: 0.828
Model: OLS Adj. R-squared: 0.820
Method: Least Squares F-statistic: 101.2
Date: Mon, 04 Apr 2016 Prob (F-statistic): 8.65e-17
Time: 23:00:11 Log-Likelihood: -178.98
No. Observations: 45 AIC: 364.0
Df Residuals: 42 BIC: 369.4
Df Model: 2
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept -6.0647 4.272 -1.420 0.163 -14.686 2.556
income 0.5987 0.120 5.003 0.000 0.357 0.840
education 0.5458 0.098 5.555 0.000 0.348 0.744
==============================================================================
Omnibus: 1.279 Durbin-Watson: 1.458
Prob(Omnibus): 0.528 Jarque-Bera (JB): 0.520
Skew: 0.155 Prob(JB): 0.771
Kurtosis: 3.426 Cond. No. 163.
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
In [40]:
infl = ols_model.get_influence()
student = infl.summary_frame()['student_resid']
print(student)
accountant 0.303900 pilot 0.340920 architect 0.072256 author 0.000711 chemist 0.826578 minister 3.134519 professor 0.768277 dentist -0.498082 reporter -2.397022 engineer 0.306225 undertaker -0.187339 lawyer -0.303082 physician 0.355687 welfare.worker -0.411406 teacher 0.050510 conductor -1.704032 contractor 2.043805 factory.owner 1.602429 store.manager 0.142425 banker 0.508388 bookkeeper -0.902388 mail.carrier -1.433249 insurance.agent -1.930919 store.clerk -1.760491 carpenter 1.068858 electrician 0.731949 RR.engineer 0.808922 machinist 1.887047 auto.repairman 0.522735 plumber -0.377954 gas.stn.attendant -0.666596 coal.miner 1.018527 streetcar.motorman -1.104485 taxi.driver 0.023322 truck.driver -0.129227 machine.operator 0.499922 barber 0.173805 bartender -0.902422 shoe.shiner -0.429357 cook 0.127207 soda.clerk -0.883095 watchman -0.513502 janitor -0.079890 policeman 0.078847 waiter -0.475972 Name: student_resid, dtype: float64
In [41]:
print(student.ix[np.abs(student) > 2])
minister 3.134519 reporter -2.397022 contractor 2.043805 Name: student_resid, dtype: float64
In [42]:
print(infl.summary_frame().ix['minister'])
dfb_Intercept 0.144937 dfb_income -1.220939 dfb_education 1.263019 cooks_d 0.566380 dffits 1.433935 dffits_internal 1.303510 hat_diag 0.173058 standard_resid 2.849416 student_resid 3.134519 Name: minister, dtype: float64
In [43]:
sidak = ols_model.outlier_test('sidak')
sidak.sort('unadj_p', inplace=True)
print(sidak)
student_resid unadj_p sidak(p)
minister 3.134519 0.003177 0.133421
reporter -2.397022 0.021170 0.618213
contractor 2.043805 0.047433 0.887721
insurance.agent -1.930919 0.060428 0.939485
machinist 1.887047 0.066248 0.954247
store.clerk -1.760491 0.085783 0.982331
conductor -1.704032 0.095944 0.989315
factory.owner 1.602429 0.116738 0.996250
mail.carrier -1.433249 0.159369 0.999595
streetcar.motorman -1.104485 0.275823 1.000000
carpenter 1.068858 0.291386 1.000000
coal.miner 1.018527 0.314400 1.000000
bartender -0.902422 0.372104 1.000000
bookkeeper -0.902388 0.372122 1.000000
soda.clerk -0.883095 0.382334 1.000000
chemist 0.826578 0.413261 1.000000
RR.engineer 0.808922 0.423229 1.000000
professor 0.768277 0.446725 1.000000
electrician 0.731949 0.468363 1.000000
gas.stn.attendant -0.666596 0.508764 1.000000
auto.repairman 0.522735 0.603972 1.000000
watchman -0.513502 0.610357 1.000000
banker 0.508388 0.613906 1.000000
machine.operator 0.499922 0.619802 1.000000
dentist -0.498082 0.621088 1.000000
waiter -0.475972 0.636621 1.000000
shoe.shiner -0.429357 0.669912 1.000000
welfare.worker -0.411406 0.682918 1.000000
plumber -0.377954 0.707414 1.000000
physician 0.355687 0.723898 1.000000
pilot 0.340920 0.734905 1.000000
engineer 0.306225 0.760983 1.000000
accountant 0.303900 0.762741 1.000000
lawyer -0.303082 0.763360 1.000000
undertaker -0.187339 0.852319 1.000000
barber 0.173805 0.862874 1.000000
store.manager 0.142425 0.887442 1.000000
truck.driver -0.129227 0.897810 1.000000
cook 0.127207 0.899399 1.000000
janitor -0.079890 0.936713 1.000000
policeman 0.078847 0.937538 1.000000
architect 0.072256 0.942750 1.000000
teacher 0.050510 0.959961 1.000000
taxi.driver 0.023322 0.981507 1.000000
author 0.000711 0.999436 1.000000
In [44]:
fdr = ols_model.outlier_test('fdr_bh')
fdr.sort('unadj_p', inplace=True)
print(fdr)
student_resid unadj_p fdr_bh(p)
minister 3.134519 0.003177 0.142974
reporter -2.397022 0.021170 0.476332
contractor 2.043805 0.047433 0.596233
insurance.agent -1.930919 0.060428 0.596233
machinist 1.887047 0.066248 0.596233
store.clerk -1.760491 0.085783 0.616782
conductor -1.704032 0.095944 0.616782
factory.owner 1.602429 0.116738 0.656653
mail.carrier -1.433249 0.159369 0.796844
streetcar.motorman -1.104485 0.275823 0.999436
carpenter 1.068858 0.291386 0.999436
coal.miner 1.018527 0.314400 0.999436
bartender -0.902422 0.372104 0.999436
bookkeeper -0.902388 0.372122 0.999436
soda.clerk -0.883095 0.382334 0.999436
chemist 0.826578 0.413261 0.999436
RR.engineer 0.808922 0.423229 0.999436
professor 0.768277 0.446725 0.999436
electrician 0.731949 0.468363 0.999436
gas.stn.attendant -0.666596 0.508764 0.999436
auto.repairman 0.522735 0.603972 0.999436
watchman -0.513502 0.610357 0.999436
banker 0.508388 0.613906 0.999436
machine.operator 0.499922 0.619802 0.999436
dentist -0.498082 0.621088 0.999436
waiter -0.475972 0.636621 0.999436
shoe.shiner -0.429357 0.669912 0.999436
welfare.worker -0.411406 0.682918 0.999436
plumber -0.377954 0.707414 0.999436
physician 0.355687 0.723898 0.999436
pilot 0.340920 0.734905 0.999436
engineer 0.306225 0.760983 0.999436
accountant 0.303900 0.762741 0.999436
lawyer -0.303082 0.763360 0.999436
undertaker -0.187339 0.852319 0.999436
barber 0.173805 0.862874 0.999436
store.manager 0.142425 0.887442 0.999436
truck.driver -0.129227 0.897810 0.999436
cook 0.127207 0.899399 0.999436
janitor -0.079890 0.936713 0.999436
policeman 0.078847 0.937538 0.999436
architect 0.072256 0.942750 0.999436
teacher 0.050510 0.959961 0.999436
taxi.driver 0.023322 0.981507 0.999436
author 0.000711 0.999436 0.999436
In [45]:
rlm_model = rlm('prestige ~ income + education', prestige).fit()
print(rlm_model.summary())
Robust linear Model Regression Results
==============================================================================
Dep. Variable: prestige No. Observations: 45
Model: RLM Df Residuals: 42
Method: IRLS Df Model: 2
Norm: HuberT
Scale Est.: mad
Cov Type: H1
Date: Mon, 04 Apr 2016
Time: 23:00:13
No. Iterations: 18
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept -7.1107 3.879 -1.833 0.067 -14.713 0.492
income 0.7015 0.109 6.456 0.000 0.489 0.914
education 0.4854 0.089 5.441 0.000 0.311 0.660
==============================================================================
If the model instance has been used for another fit with different fit
parameters, then the fit options might not be the correct ones anymore .
In [46]:
print(rlm_model.weights)
accountant 1.000000 pilot 1.000000 architect 1.000000 author 1.000000 chemist 1.000000 minister 0.344596 professor 1.000000 dentist 1.000000 reporter 0.441669 engineer 1.000000 undertaker 1.000000 lawyer 1.000000 physician 1.000000 welfare.worker 1.000000 teacher 1.000000 conductor 0.538445 contractor 0.552262 factory.owner 0.706169 store.manager 1.000000 banker 1.000000 bookkeeper 1.000000 mail.carrier 0.690764 insurance.agent 0.533499 store.clerk 0.618656 carpenter 0.935848 electrician 1.000000 RR.engineer 1.000000 machinist 0.570360 auto.repairman 1.000000 plumber 1.000000 gas.stn.attendant 1.000000 coal.miner 0.963821 streetcar.motorman 0.832870 taxi.driver 1.000000 truck.driver 1.000000 machine.operator 1.000000 barber 1.000000 bartender 1.000000 shoe.shiner 1.000000 cook 1.000000 soda.clerk 1.000000 watchman 1.000000 janitor 1.000000 policeman 1.000000 waiter 1.000000 dtype: float64
Hertzprung Russell data for Star Cluster CYG 0B1 - Leverage Points¶
- Data is on the luminosity and temperature of 47 stars in the direction of Cygnus.
In [47]:
dta = sm.datasets.get_rdataset("starsCYG", "robustbase", cache=True).data
In [48]:
from matplotlib.patches import Ellipse
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111, xlabel='log(Temp)', ylabel='log(Light)', title='Hertzsprung-Russell Diagram of Star Cluster CYG OB1')
ax.scatter(*dta.values.T)
# highlight outliers
e = Ellipse((3.5, 6), .2, 1, alpha=.25, color='r')
ax.add_patch(e);
ax.annotate('Red giants', xy=(3.6, 6), xytext=(3.8, 6),
arrowprops=dict(facecolor='black', shrink=0.05, width=2),
horizontalalignment='left', verticalalignment='bottom',
clip_on=True, # clip to the axes bounding box
fontsize=16,
)
# annotate these with their index
for i,row in dta.ix[dta['log.Te'] < 3.8].iterrows():
ax.annotate(i, row, row + .01, fontsize=14)
xlim, ylim = ax.get_xlim(), ax.get_ylim()
Error in callback <function post_execute at 0xb082b2b0> (for post_execute):
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) /usr/lib/python2.7/dist-packages/matplotlib/pyplot.pyc in post_execute() 145 def post_execute(): 146 if matplotlib.is_interactive(): --> 147 draw_all() 148 149 # IPython >= 2 /usr/lib/python2.7/dist-packages/matplotlib/_pylab_helpers.pyc in draw_all(cls, force) 148 for f_mgr in cls.get_all_fig_managers(): 149 if force or f_mgr.canvas.figure.stale: --> 150 f_mgr.canvas.draw_idle() 151 152 atexit.register(Gcf.destroy_all) /usr/lib/python2.7/dist-packages/matplotlib/backend_bases.pyc in draw_idle(self, *args, **kwargs) 2024 if not self._is_idle_drawing: 2025 with self._idle_draw_cntx(): -> 2026 self.draw(*args, **kwargs) 2027 2028 def draw_cursor(self, event): /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in draw(self) 472 473 try: --> 474 self.figure.draw(self.renderer) 475 finally: 476 RendererAgg.lock.release() /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/figure.pyc in draw(self, renderer) 1157 dsu.sort(key=itemgetter(0)) 1158 for zorder, a, func, args in dsu: -> 1159 func(*args) 1160 1161 renderer.close_group('figure') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in draw(self, renderer, inframe) 2322 2323 for zorder, a in dsu: -> 2324 a.draw(renderer) 2325 2326 renderer.close_group('axes') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in draw(self, renderer, *args, **kwargs) 1106 ticks_to_draw = self._update_ticks(renderer) 1107 ticklabelBoxes, ticklabelBoxes2 = self._get_tick_bboxes(ticks_to_draw, -> 1108 renderer) 1109 1110 for tick in ticks_to_draw: /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in _get_tick_bboxes(self, ticks, renderer) 1056 for tick in ticks: 1057 if tick.label1On and tick.label1.get_visible(): -> 1058 extent = tick.label1.get_window_extent(renderer) 1059 ticklabelBoxes.append(extent) 1060 if tick.label2On and tick.label2.get_visible(): /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in get_window_extent(self, renderer, dpi) 959 raise RuntimeError('Cannot get window extent w/o renderer') 960 --> 961 bbox, info, descent = self._get_layout(self._renderer) 962 x, y = self.get_unitless_position() 963 x, y = self.get_transform().transform_point((x, y)) /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in _get_layout(self, renderer) 350 tmp, lp_h, lp_bl = renderer.get_text_width_height_descent('lp', 351 self._fontproperties, --> 352 ismath=False) 353 offsety = (lp_h - lp_bl) * self._linespacing 354 /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in get_text_width_height_descent(self, s, prop, ismath) 227 fontsize = prop.get_size_in_points() 228 w, h, d = texmanager.get_text_width_height_descent(s, fontsize, --> 229 renderer=self) 230 return w, h, d 231 /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in get_text_width_height_descent(self, tex, fontsize, renderer) 673 else: 674 # use dviread. It sometimes returns a wrong descent. --> 675 dvifile = self.make_dvi(tex, fontsize) 676 dvi = dviread.Dvi(dvifile, 72 * dpi_fraction) 677 try: /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in make_dvi(self, tex, fontsize) 420 'string:\n%s\nHere is the full report generated by ' 421 'LaTeX: \n\n' % repr(tex.encode('unicode_escape')) + --> 422 report)) 423 else: 424 mpl.verbose.report(report, 'debug') RuntimeError: LaTeX was not able to process the following string: 'lp' Here is the full report generated by LaTeX:
<matplotlib.figure.Figure at 0xa7410530>
In [49]:
from IPython.display import Image
Image(filename='star_diagram.png')
Out[49]:
In [50]:
y = dta['log.light']
X = sm.add_constant(dta['log.Te'], prepend=True)
ols_model = sm.OLS(y, X).fit()
abline_plot(model_results=ols_model, ax=ax)
Out[50]:
<matplotlib.figure.Figure at 0xa7410530>
In [51]:
rlm_mod = sm.RLM(y, X, sm.robust.norms.TrimmedMean(.5)).fit()
abline_plot(model_results=rlm_mod, ax=ax, color='red')
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) <ipython-input-264-9210cc7d7045> in <module>() 1 rlm_mod = sm.RLM(y, X, sm.robust.norms.TrimmedMean(.5)).fit() ----> 2 abline_plot(model_results=rlm_mod, ax=ax, color='red') /statsmodels-0.6.1/debian/python-statsmodels/usr/lib/python2.7/dist-packages/statsmodels/graphics/regressionplots.pyc in abline_plot(intercept, slope, horiz, vert, model_results, ax, **kwargs) 665 666 data_y = [x[0]*slope+intercept, x[1]*slope+intercept] --> 667 ax.set_xlim(x) 668 #ax.set_ylim(y) 669 /usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in set_xlim(self, left, right, emit, auto, **kw) 2794 2795 if emit: -> 2796 self.callbacks.process('xlim_changed', self) 2797 # Call all of the other x-axes that are shared with this one 2798 for other in self._shared_x_axes.get_siblings(self): /usr/lib/python2.7/dist-packages/matplotlib/cbook.pyc in process(self, s, *args, **kwargs) 561 for cid, proxy in list(six.iteritems(self.callbacks[s])): 562 try: --> 563 proxy(*args, **kwargs) 564 except ReferenceError: 565 self._remove_proxy(proxy) /usr/lib/python2.7/dist-packages/matplotlib/cbook.pyc in __call__(self, *args, **kwargs) 428 mtd = self.func 429 # invoke the callable and return the result --> 430 return mtd(*args, **kwargs) 431 432 def __eq__(self, other): /statsmodels-0.6.1/debian/python-statsmodels/usr/lib/python2.7/dist-packages/statsmodels/graphics/regressionplots.pyc in update_datalim(self, ax) 681 y = [x[0]*slope+intercept, x[1]*slope+intercept] 682 abline.set_data(x, y) --> 683 ax.figure.canvas.draw() 684 #TODO: how to intercept something like a margins call and adjust? 685 /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in draw(self) 472 473 try: --> 474 self.figure.draw(self.renderer) 475 finally: 476 RendererAgg.lock.release() /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/figure.pyc in draw(self, renderer) 1157 dsu.sort(key=itemgetter(0)) 1158 for zorder, a, func, args in dsu: -> 1159 func(*args) 1160 1161 renderer.close_group('figure') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in draw(self, renderer, inframe) 2322 2323 for zorder, a in dsu: -> 2324 a.draw(renderer) 2325 2326 renderer.close_group('axes') /usr/lib/python2.7/dist-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs) 59 def draw_wrapper(artist, renderer, *args, **kwargs): 60 before(artist, renderer) ---> 61 draw(artist, renderer, *args, **kwargs) 62 after(artist, renderer) 63 /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in draw(self, renderer, *args, **kwargs) 1106 ticks_to_draw = self._update_ticks(renderer) 1107 ticklabelBoxes, ticklabelBoxes2 = self._get_tick_bboxes(ticks_to_draw, -> 1108 renderer) 1109 1110 for tick in ticks_to_draw: /usr/lib/python2.7/dist-packages/matplotlib/axis.pyc in _get_tick_bboxes(self, ticks, renderer) 1056 for tick in ticks: 1057 if tick.label1On and tick.label1.get_visible(): -> 1058 extent = tick.label1.get_window_extent(renderer) 1059 ticklabelBoxes.append(extent) 1060 if tick.label2On and tick.label2.get_visible(): /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in get_window_extent(self, renderer, dpi) 959 raise RuntimeError('Cannot get window extent w/o renderer') 960 --> 961 bbox, info, descent = self._get_layout(self._renderer) 962 x, y = self.get_unitless_position() 963 x, y = self.get_transform().transform_point((x, y)) /usr/lib/python2.7/dist-packages/matplotlib/text.pyc in _get_layout(self, renderer) 350 tmp, lp_h, lp_bl = renderer.get_text_width_height_descent('lp', 351 self._fontproperties, --> 352 ismath=False) 353 offsety = (lp_h - lp_bl) * self._linespacing 354 /usr/lib/python2.7/dist-packages/matplotlib/backends/backend_agg.pyc in get_text_width_height_descent(self, s, prop, ismath) 227 fontsize = prop.get_size_in_points() 228 w, h, d = texmanager.get_text_width_height_descent(s, fontsize, --> 229 renderer=self) 230 return w, h, d 231 /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in get_text_width_height_descent(self, tex, fontsize, renderer) 673 else: 674 # use dviread. It sometimes returns a wrong descent. --> 675 dvifile = self.make_dvi(tex, fontsize) 676 dvi = dviread.Dvi(dvifile, 72 * dpi_fraction) 677 try: /usr/lib/python2.7/dist-packages/matplotlib/texmanager.pyc in make_dvi(self, tex, fontsize) 420 'string:\n%s\nHere is the full report generated by ' 421 'LaTeX: \n\n' % repr(tex.encode('unicode_escape')) + --> 422 report)) 423 else: 424 mpl.verbose.report(report, 'debug') RuntimeError: LaTeX was not able to process the following string: 'lp' Here is the full report generated by LaTeX:
- Why? Because M-estimators are not robust to leverage points.
In [52]:
infl = ols_model.get_influence()
In [53]:
h_bar = 2*(ols_model.df_model + 1 )/ols_model.nobs
hat_diag = infl.summary_frame()['hat_diag']
hat_diag.ix[hat_diag > h_bar]
Out[53]:
10 0.194103 19 0.194103 29 0.198344 33 0.194103 Name: hat_diag, dtype: float64
In [54]:
sidak2 = ols_model.outlier_test('sidak')
sidak2.sort('unadj_p', inplace=True)
print(sidak2)
student_resid unadj_p sidak(p)
16 -2.049393 0.046415 0.892872
13 -2.035329 0.047868 0.900286
33 1.905847 0.063216 0.953543
18 -1.575505 0.122304 0.997826
1 1.522185 0.135118 0.998911
3 1.522185 0.135118 0.998911
21 -1.450418 0.154034 0.999615
17 -1.426675 0.160731 0.999735
29 1.388520 0.171969 0.999859
14 -1.374733 0.176175 0.999889
35 1.346543 0.185023 0.999933
34 -1.272159 0.209999 0.999985
28 -1.186946 0.241618 0.999998
20 -1.150621 0.256103 0.999999
44 1.134779 0.262612 0.999999
39 1.091886 0.280826 1.000000
19 1.064878 0.292740 1.000000
6 -1.026873 0.310093 1.000000
30 -1.009096 0.318446 1.000000
22 -0.979768 0.332557 1.000000
8 0.961218 0.341695 1.000000
5 0.913802 0.365801 1.000000
11 0.871997 0.387943 1.000000
12 0.856685 0.396261 1.000000
46 -0.833923 0.408829 1.000000
10 0.743920 0.460879 1.000000
42 0.727179 0.470968 1.000000
15 -0.689258 0.494280 1.000000
43 0.688272 0.494895 1.000000
7 0.655712 0.515424 1.000000
40 -0.646396 0.521381 1.000000
26 -0.640978 0.524862 1.000000
25 -0.545561 0.588123 1.000000
37 0.472819 0.638680 1.000000
32 0.472819 0.638680 1.000000
38 0.462187 0.646225 1.000000
0 0.430686 0.668799 1.000000
31 0.341726 0.734184 1.000000
36 0.318911 0.751303 1.000000
4 0.307890 0.759619 1.000000
9 0.235114 0.815211 1.000000
41 0.187732 0.851950 1.000000
2 -0.182093 0.856346 1.000000
23 -0.156014 0.876736 1.000000
27 -0.147406 0.883485 1.000000
24 0.065195 0.948314 1.000000
45 0.045675 0.963776 1.000000
In [55]:
fdr2 = ols_model.outlier_test('fdr_bh')
fdr2.sort('unadj_p', inplace=True)
print(fdr2)
student_resid unadj_p fdr_bh(p)
16 -2.049393 0.046415 0.764747
13 -2.035329 0.047868 0.764747
33 1.905847 0.063216 0.764747
18 -1.575505 0.122304 0.764747
1 1.522185 0.135118 0.764747
3 1.522185 0.135118 0.764747
21 -1.450418 0.154034 0.764747
17 -1.426675 0.160731 0.764747
29 1.388520 0.171969 0.764747
14 -1.374733 0.176175 0.764747
35 1.346543 0.185023 0.764747
34 -1.272159 0.209999 0.764747
28 -1.186946 0.241618 0.764747
20 -1.150621 0.256103 0.764747
44 1.134779 0.262612 0.764747
39 1.091886 0.280826 0.764747
19 1.064878 0.292740 0.764747
6 -1.026873 0.310093 0.764747
30 -1.009096 0.318446 0.764747
22 -0.979768 0.332557 0.764747
8 0.961218 0.341695 0.764747
5 0.913802 0.365801 0.768599
11 0.871997 0.387943 0.768599
12 0.856685 0.396261 0.768599
46 -0.833923 0.408829 0.768599
10 0.743920 0.460879 0.770890
42 0.727179 0.470968 0.770890
15 -0.689258 0.494280 0.770890
43 0.688272 0.494895 0.770890
7 0.655712 0.515424 0.770890
40 -0.646396 0.521381 0.770890
26 -0.640978 0.524862 0.770890
25 -0.545561 0.588123 0.837630
37 0.472819 0.638680 0.843682
32 0.472819 0.638680 0.843682
38 0.462187 0.646225 0.843682
0 0.430686 0.668799 0.849556
31 0.341726 0.734184 0.892552
36 0.318911 0.751303 0.892552
4 0.307890 0.759619 0.892552
9 0.235114 0.815211 0.922751
41 0.187732 0.851950 0.922751
2 -0.182093 0.856346 0.922751
23 -0.156014 0.876736 0.922751
27 -0.147406 0.883485 0.922751
24 0.065195 0.948314 0.963776
45 0.045675 0.963776 0.963776
- Let's delete that line
In [56]:
del ax.lines[-1]
In [57]:
weights = np.ones(len(X))
weights[X[X['log.Te'] < 3.8].index.values - 1] = 0
wls_model = sm.WLS(y, X, weights=weights).fit()
abline_plot(model_results=wls_model, ax=ax, color='green')
--------------------------------------------------------------------------- IndexError Traceback (most recent call last) <ipython-input-270-e3ec53a40864> in <module>() 2 weights[X[X['log.Te'] < 3.8].index.values - 1] = 0 3 wls_model = sm.WLS(y, X, weights=weights).fit() ----> 4 abline_plot(model_results=wls_model, ax=ax, color='green') /statsmodels-0.6.1/debian/python-statsmodels/usr/lib/python2.7/dist-packages/statsmodels/graphics/regressionplots.pyc in abline_plot(intercept, slope, horiz, vert, model_results, ax, **kwargs) 665 666 data_y = [x[0]*slope+intercept, x[1]*slope+intercept] --> 667 ax.set_xlim(x) 668 #ax.set_ylim(y) 669 /usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in set_xlim(self, left, right, emit, auto, **kw) 2794 2795 if emit: -> 2796 self.callbacks.process('xlim_changed', self) 2797 # Call all of the other x-axes that are shared with this one 2798 for other in self._shared_x_axes.get_siblings(self): /usr/lib/python2.7/dist-packages/matplotlib/cbook.pyc in process(self, s, *args, **kwargs) 561 for cid, proxy in list(six.iteritems(self.callbacks[s])): 562 try: --> 563 proxy(*args, **kwargs) 564 except ReferenceError: 565 self._remove_proxy(proxy) /usr/lib/python2.7/dist-packages/matplotlib/cbook.pyc in __call__(self, *args, **kwargs) 428 mtd = self.func 429 # invoke the callable and return the result --> 430 return mtd(*args, **kwargs) 431 432 def __eq__(self, other): /statsmodels-0.6.1/debian/python-statsmodels/usr/lib/python2.7/dist-packages/statsmodels/graphics/regressionplots.pyc in update_datalim(self, ax) 677 children = ax.get_children() 678 abline = [children[i] for i in range(len(children)) --> 679 if isinstance(children[i], ABLine2D)][0] 680 x = ax.get_xlim() 681 y = [x[0]*slope+intercept, x[1]*slope+intercept] IndexError: list index out of range
- MM estimators are good for this type of problem, unfortunately, we don't yet have these yet.
- It's being worked on, but it gives a good excuse to look at the R cell magics in the notebook.
In [58]:
yy = y.values[:,None]
xx = X['log.Te'].values[:,None]
In [59]:
%load_ext rmagic
%R library(robustbase)
%Rpush yy xx
%R mod <- lmrob(yy ~ xx);
%R params <- mod$coefficients;
%Rpull params
--------------------------------------------------------------------------- ImportError Traceback (most recent call last) <ipython-input-272-53d6967d0401> in <module>() ----> 1 get_ipython().magic(u'load_ext rmagic') 2 3 get_ipython().magic(u'R library(robustbase)') 4 get_ipython().magic(u'Rpush yy xx') 5 get_ipython().magic(u'R mod <- lmrob(yy ~ xx);') /usr/lib/python2.7/dist-packages/IPython/core/interactiveshell.pyc in magic(self, arg_s) 2203 magic_name, _, magic_arg_s = arg_s.partition(' ') 2204 magic_name = magic_name.lstrip(prefilter.ESC_MAGIC) -> 2205 return self.run_line_magic(magic_name, magic_arg_s) 2206 2207 #------------------------------------------------------------------------- /usr/lib/python2.7/dist-packages/IPython/core/interactiveshell.pyc in run_line_magic(self, magic_name, line) 2124 kwargs['local_ns'] = sys._getframe(stack_depth).f_locals 2125 with self.builtin_trap: -> 2126 result = fn(*args,**kwargs) 2127 return result 2128 <decorator-gen-64> in load_ext(self, module_str) /usr/lib/python2.7/dist-packages/IPython/core/magic.pyc in <lambda>(f, *a, **k) 191 # but it's overkill for just that one bit of state. 192 def magic_deco(arg): --> 193 call = lambda f, *a, **k: f(*a, **k) 194 195 if callable(arg): /usr/lib/python2.7/dist-packages/IPython/core/magics/extension.pyc in load_ext(self, module_str) 61 if not module_str: 62 raise UsageError('Missing module name.') ---> 63 res = self.shell.extension_manager.load_extension(module_str) 64 65 if res == 'already loaded': /usr/lib/python2.7/dist-packages/IPython/core/extensions.pyc in load_extension(self, module_str) 96 if module_str not in sys.modules: 97 with prepended_to_syspath(self.ipython_extension_dir): ---> 98 __import__(module_str) 99 mod = sys.modules[module_str] 100 if self._call_load_ipython_extension(mod): /usr/lib/python2.7/dist-packages/IPython/extensions/rmagic.py in <module>() 54 import numpy as np 55 ---> 56 import rpy2.rinterface as ri 57 import rpy2.robjects as ro 58 try: ImportError: No module named rpy2.rinterface
In [60]:
%R print(mod)
ERROR: Line magic function `%R` not found.
In [61]:
print(params)
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-274-73ac4b936803> in <module>() ----> 1 print(params) NameError: name 'params' is not defined
In [62]:
abline_plot(intercept=params[0], slope=params[1], ax=ax, color='green')
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-275-e1a9f35b3320> in <module>() ----> 1 abline_plot(intercept=params[0], slope=params[1], ax=ax, color='green') NameError: name 'params' is not defined
Exercise: Breakdown points of M-estimator¶
In [63]:
np.random.seed(12345)
nobs = 200
beta_true = np.array([3, 1, 2.5, 3, -4])
X = np.random.uniform(-20,20, size=(nobs, len(beta_true)-1))
# stack a constant in front
X = sm.add_constant(X, prepend=True) # np.c_[np.ones(nobs), X]
mc_iter = 500
contaminate = .25 # percentage of response variables to contaminate
In [64]:
all_betas = []
for i in range(mc_iter):
y = np.dot(X, beta_true) + np.random.normal(size=200)
random_idx = np.random.randint(0, nobs, size=int(contaminate * nobs))
y[random_idx] = np.random.uniform(-750, 750)
beta_hat = sm.RLM(y, X).fit().params
all_betas.append(beta_hat)
In [65]:
all_betas = np.asarray(all_betas)
se_loss = lambda x : np.linalg.norm(x, ord=2)**2
se_beta = lmap(se_loss, all_betas - beta_true)
Squared error loss¶
In [66]:
np.array(se_beta).mean()
Out[66]:
0.44502948730686182
In [67]:
all_betas.mean(0)
Out[67]:
array([ 2.9971, 0.999 , 2.4991, 2.9971, -3.9963])
In [68]:
beta_true
Out[68]:
array([ 3. , 1. , 2.5, 3. , -4. ])
In [69]:
se_loss(all_betas.mean(0) - beta_true)
Out[69]:
3.2360913286754188e-05