Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1# -*- coding: utf-8 -*- 

2""" 

3Created on Sun May 10 08:23:48 2015 

4 

5Author: Josef Perktold 

6License: BSD-3 

7""" 

8 

9import numpy as np 

10from ._penalties import NonePenalty 

11from statsmodels.tools.numdiff import approx_fprime_cs, approx_fprime 

12 

13 

14class PenalizedMixin(object): 

15 """Mixin class for Maximum Penalized Likelihood 

16 

17 Parameters 

18 ---------- 

19 args and kwds for the model super class 

20 penal : None or instance of Penalized function class 

21 If penal is None, then NonePenalty is used. 

22 pen_weight : float or None 

23 factor for weighting the penalization term. 

24 If None, then pen_weight is set to nobs. 

25 

26 

27 TODO: missing **kwds or explicit keywords 

28 

29 TODO: do we adjust the inherited docstrings? 

30 We would need templating to add the penalization parameters 

31 """ 

32 

33 def __init__(self, *args, **kwds): 

34 

35 # pop extra kwds before calling super 

36 self.penal = kwds.pop('penal', None) 

37 self.pen_weight = kwds.pop('pen_weight', None) 

38 

39 super(PenalizedMixin, self).__init__(*args, **kwds) 

40 

41 # TODO: define pen_weight as average pen_weight? i.e. per observation 

42 # I would have prefered len(self.endog) * kwds.get('pen_weight', 1) 

43 # or use pen_weight_factor in signature 

44 if self.pen_weight is None: 

45 self.pen_weight = len(self.endog) 

46 

47 if self.penal is None: 

48 # unpenalized by default 

49 self.penal = NonePenalty() 

50 self.pen_weight = 0 

51 

52 self._init_keys.extend(['penal', 'pen_weight']) 

53 self._null_drop_keys = getattr(self, '_null_drop_keys', []) 

54 self._null_drop_keys.extend(['penal']) 

55 

56 def _handle_scale(self, params, scale=None, **kwds): 

57 

58 if scale is None: 

59 # special handling for GLM 

60 if hasattr(self, 'scaletype'): 

61 mu = self.predict(params) 

62 scale = self.estimate_scale(mu) 

63 else: 

64 scale = 1 

65 

66 return scale 

67 

68 def loglike(self, params, pen_weight=None, **kwds): 

69 """ 

70 Log-likelihood of model at params 

71 """ 

72 if pen_weight is None: 

73 pen_weight = self.pen_weight 

74 

75 llf = super(PenalizedMixin, self).loglike(params, **kwds) 

76 if pen_weight != 0: 

77 scale = self._handle_scale(params, **kwds) 

78 llf -= 1/scale * pen_weight * self.penal.func(params) 

79 

80 return llf 

81 

82 def loglikeobs(self, params, pen_weight=None, **kwds): 

83 """ 

84 Log-likelihood of model observations at params 

85 """ 

86 if pen_weight is None: 

87 pen_weight = self.pen_weight 

88 

89 llf = super(PenalizedMixin, self).loglikeobs(params, **kwds) 

90 nobs_llf = float(llf.shape[0]) 

91 

92 if pen_weight != 0: 

93 scale = self._handle_scale(params, **kwds) 

94 llf -= 1/scale * pen_weight / nobs_llf * self.penal.func(params) 

95 

96 return llf 

97 

98 def score_numdiff(self, params, pen_weight=None, method='fd', **kwds): 

99 """score based on finite difference derivative 

100 """ 

101 if pen_weight is None: 

102 pen_weight = self.pen_weight 

103 

104 loglike = lambda p: self.loglike(p, pen_weight=pen_weight, **kwds) 

105 

106 if method == 'cs': 

107 return approx_fprime_cs(params, loglike) 

108 elif method == 'fd': 

109 return approx_fprime(params, loglike, centered=True) 

110 else: 

111 raise ValueError('method not recognized, should be "fd" or "cs"') 

112 

113 def score(self, params, pen_weight=None, **kwds): 

114 """ 

115 Gradient of model at params 

116 """ 

117 if pen_weight is None: 

118 pen_weight = self.pen_weight 

119 

120 sc = super(PenalizedMixin, self).score(params, **kwds) 

121 if pen_weight != 0: 

122 scale = self._handle_scale(params, **kwds) 

123 sc -= 1/scale * pen_weight * self.penal.deriv(params) 

124 

125 return sc 

126 

127 def score_obs(self, params, pen_weight=None, **kwds): 

128 """ 

129 Gradient of model observations at params 

130 """ 

131 if pen_weight is None: 

132 pen_weight = self.pen_weight 

133 

134 sc = super(PenalizedMixin, self).score_obs(params, **kwds) 

135 nobs_sc = float(sc.shape[0]) 

136 if pen_weight != 0: 

137 scale = self._handle_scale(params, **kwds) 

138 sc -= 1/scale * pen_weight / nobs_sc * self.penal.deriv(params) 

139 

140 return sc 

141 

142 def hessian_numdiff(self, params, pen_weight=None, **kwds): 

143 """hessian based on finite difference derivative 

144 """ 

145 if pen_weight is None: 

146 pen_weight = self.pen_weight 

147 loglike = lambda p: self.loglike(p, pen_weight=pen_weight, **kwds) 

148 

149 from statsmodels.tools.numdiff import approx_hess 

150 return approx_hess(params, loglike) 

151 

152 def hessian(self, params, pen_weight=None, **kwds): 

153 """ 

154 Hessian of model at params 

155 """ 

156 if pen_weight is None: 

157 pen_weight = self.pen_weight 

158 

159 hess = super(PenalizedMixin, self).hessian(params, **kwds) 

160 if pen_weight != 0: 

161 scale = self._handle_scale(params, **kwds) 

162 h = self.penal.deriv2(params) 

163 if h.ndim == 1: 

164 hess -= 1/scale * np.diag(pen_weight * h) 

165 else: 

166 hess -= 1/scale * pen_weight * h 

167 

168 return hess 

169 

170 def fit(self, method=None, trim=None, **kwds): 

171 """minimize negative penalized log-likelihood 

172 

173 Parameters 

174 ---------- 

175 method : None or str 

176 Method specifies the scipy optimizer as in nonlinear MLE models. 

177 trim : {bool, float} 

178 Default is False or None, which uses no trimming. 

179 If trim is True or a float, then small parameters are set to zero. 

180 If True, then a default threshold is used. If trim is a float, then 

181 it will be used as threshold. 

182 The default threshold is currently 1e-4, but it will change in 

183 future and become penalty function dependent. 

184 kwds : extra keyword arguments 

185 This keyword arguments are treated in the same way as in the 

186 fit method of the underlying model class. 

187 Specifically, additional optimizer keywords and cov_type related 

188 keywords can be added. 

189 """ 

190 # If method is None, then we choose a default method ourselves 

191 

192 # TODO: temporary hack, need extra fit kwds 

193 # we need to rule out fit methods in a model that will not work with 

194 # penalization 

195 if hasattr(self, 'family'): # assume this identifies GLM 

196 kwds.update({'max_start_irls' : 0}) 

197 

198 # currently we use `bfgs` by default 

199 if method is None: 

200 method = 'bfgs' 

201 

202 if trim is None: 

203 trim = False 

204 

205 res = super(PenalizedMixin, self).fit(method=method, **kwds) 

206 

207 if trim is False: 

208 # note boolean check for "is False", not "False_like" 

209 return res 

210 else: 

211 if trim is True: 

212 trim = 1e-4 # trim threshold 

213 # TODO: make it penal function dependent 

214 # temporary standin, only checked for Poisson and GLM, 

215 # and is computationally inefficient 

216 drop_index = np.nonzero(np.abs(res.params) < trim)[0] 

217 keep_index = np.nonzero(np.abs(res.params) > trim)[0] 

218 

219 if drop_index.any(): 

220 # TODO: do we need to add results attributes? 

221 res_aux = self._fit_zeros(keep_index, **kwds) 

222 return res_aux 

223 else: 

224 return res