Changeset 10859:08a0a35c1687 in orange for Orange/regression/lasso.py
 Timestamp:
 05/04/12 18:07:12 (2 years ago)
 Branch:
 default
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Orange/regression/lasso.py
r10580 r10859 1 from numpy import dot, std, array, zeros, maximum, sqrt, sign, log, abs, \ 2 ascontiguousarray, random as rnd 3 from scipy.linalg import norm, eigh 4 1 5 import Orange 2 import numpy3 4 from Orange.regression import base5 6 6 from Orange.utils import deprecated_members, deprecated_keywords 7 7 8 def center(X): 9 """Centers the data, i.e. subtracts the column means. 10 Returns the centered data and the mean.11 12 :param X: the data arry13 :type table: :class:`numpy.array`8 9 def get_bootstrap_sample(data): 10 """Generate a bootstrap sample of a given data set. 11 12 :param data: the original data sample 13 :type data: :class:`Orange.data.Table` 14 14 """ 15 mu = X.mean(axis=0) 16 return X  mu, mu 17 18 def standardize(X): 19 """Standardizes the data, i.e. subtracts the column means and divide by 20 standard deviation. 21 Returns the centered data, the mean, and deviations. 22 23 :param X: the data arry 24 :type table: :class:`numpy.array` 25 """ 26 mu = numpy.mean(X, axis=0) 27 std = numpy.std(X, axis=0) 28 return (X  mu) / std, mu, std 29 30 def get_bootstrap_sample(table): 31 """Generates boostrap sample from an Orange Example Table 32 and stores it in a new :class:`Orange.data.Table` object 33 34 :param table: the original data sample 35 :type table: :class:`Orange.data.Table` 36 """ 37 n = len(table) 38 bootTable = Orange.data.Table(table.domain) 39 for i in range(n): 40 id = numpy.random.randint(0, n) 41 bootTable.append(table[id]) 42 return bootTable 43 44 def permute_responses(table): 45 """ Permutes values of the class (response) variable. 15 n = len(data) 16 bootstrap = Orange.data.Table(data.domain) 17 for id in rnd.randint(0, n, n): 18 bootstrap.append(data[id]) 19 return bootstrap 20 21 def permute_responses(data): 22 """Permute values of the class (response) variable. 46 23 The independence between independent variables and the response 47 24 is obtained but the distribution of the response variable is kept. 48 25 49 :param table: the original data sample50 :type table: :class:`Orange.data.Table`26 :param data: Original data. 27 :type data: :class:`Orange.data.Table` 51 28 """ 52 n = len(table) 53 perm = numpy.random.permutation(n) 54 permTable = Orange.data.Table(table.domain, table) 55 for i, ins in enumerate(table): 56 permTable[i].set_class(table[perm[i]].get_class()) 57 return permTable 58 59 class LassoRegressionLearner(base.BaseRegressionLearner): 60 """Fits the lasso regression model, i.e. learns the regression parameters 61 The class is derived from 62 :class:`Orange.regression.base.BaseRegressionLearner` 63 which is used for preprocessing the data (continuization and imputation) 64 before fitting the regression parameters 65 29 n = len(data) 30 perm = rnd.permutation(n) 31 perm_data = Orange.data.Table(data.domain, data) 32 for i, ins in enumerate(data): 33 perm_data[i].set_class(data[perm[i]].get_class()) 34 return perm_data 35 36 class LassoRegressionLearner(Orange.regression.base.BaseRegressionLearner): 37 """Fits the lasso regression model using FISTA 38 (Fast Iterative ShrinkageThresholding Algorithm). 66 39 """ 67 40 68 def __init__(self, name='lasso regression', t=1, s=None, tol=0.001, \ 69 n_boot=100, n_perm=100, imputer=None, continuizer=None): 70 """ 71 :param name: name of the linear model, default 'lasso regression' 72 :type name: string 41 def __init__(self, lasso_lambda=0.1, max_iter=20000, eps=1e6, 42 n_boot=0, n_perm=0, imputer=None, continuizer=None, 43 name='Lasso'): 44 """ 45 :param lasso_lambda: Regularization parameter. 46 :type lasso_lambda: float 47 48 :param max_iter: Maximum number of iterations for 49 the optimization method. 50 :type max_iter: int 51 52 :param eps: Stop optimization when improvements are lower than eps. 53 :type eps: float 54 55 :param n_boot: Number of bootstrap samples used for nonparametric 56 estimation of standard errors. 57 :type n_boot: int 58 59 :param n_perm: Number of permuations used for nonparametric 60 estimation of pvalues. 61 :type n_perm: int 62 63 :param name: Learner name. 64 :type name: str 73 65 74 :param t: tuning parameter, upper bound for the L1norm of the 75 regression coefficients 76 :type t: float 77 78 :param s: An alternative way to specify the tuning parameter ``t``. 79 Here ``t`` is taken to be t = s * sum(abs(B)) where B are the 80 coefficients of an ordinary least square linear fit. ``t`` parameter is ignored if ``s`` is specified (by default it 81 is None). 82 :type s: float 83 84 :param tol: tolerance parameter, regression coefficients 85 (absoulute value) under tol are set to 0, 86 default=0.001 87 :type tol: float 88 89 :param n_boot: number of bootstrap samples used for nonparametric 90 estimation of standard errors 91 :type n_boot: int 92 93 :param n_perm: number of permuations used for nonparametric 94 estimation of pvalues 95 :type n_perm: int 96 97 """ 98 99 self.name = name 100 self.t = t 101 self.s = s 102 self.tol = tol 66 """ 67 self.lasso_lambda = lasso_lambda 68 self.max_iter = max_iter 69 self.eps = eps 103 70 self.n_boot = n_boot 104 71 self.n_perm = n_perm 105 72 self.set_imputer(imputer=imputer) 106 73 self.set_continuizer(continuizer=continuizer) 107 108 109 def __call__(self, table, weight=None): 110 """ 111 :param table: data instances. 112 :type table: :class:`Orange.data.Table` 113 :param weight: the weights for instances. Default: None, i.e. 114 all data instances are eqaully important in fitting 115 the regression parameters 116 :type weight: None or list of Orange.feature.Continuous 117 which stores weights for instances 74 self.name = name 75 76 def get_lipschitz(self, X): 77 """Return the Lipschitz constant of :math:`\\nabla f`, 78 where :math:`f(w) = \\frac{1}{2}Xwy^2`. 79 """ 80 n, m = X.shape 81 if n > m: 82 X = ascontiguousarray(X.T) 83 k = min(n, m)  1 84 eigvals = eigh(dot(X, X.T), eigvals_only=True, eigvals=(k, k)) 85 return eigvals[1] 86 87 def fista(self, X, y, l, lipschitz, w_init=None): 88 """Fast Iterative ShrinkageThresholding Algorithm (FISTA).""" 89 z = w_old = zeros(X.shape[1]) if w_init is None else w_init 90 t_old, obj_old = 1, 1e400 91 XT = ascontiguousarray(X.T) 92 for i in range(self.max_iter): 93 z = 1. / lipschitz * dot(XT, dot(X, z)  y) 94 w = maximum(0, abs(z)  l / lipschitz) * sign(z) 95 t = (1 + sqrt(1 + 4 * t_old**2)) / 2 96 z = w + (t_old  1) / t * (w  w_old) 97 obj = ((y  dot(X, w))**2).sum() + l * norm(w, 1) 98 if abs(obj_old  obj) / obj < self.eps: 99 stop += 1 100 if obj < obj_old and stop > log(i + 1): 101 break 102 else: 103 stop = 0 104 w_old, t_old = w, t 105 obj_old = obj 106 return w 107 108 def __call__(self, data, weight=None): 109 """ 110 :param data: Training data. 111 :type data: :class:`Orange.data.Table` 112 :param weight: Weights for instances. Not implemented yet. 118 113 119 114 """ 120 115 # dicrete values are continuized 121 table = self.continuize_table(table)116 data = self.continuize_table(data) 122 117 # missing values are imputed 123 table = self.impute_table(table) 124 125 domain = table.domain 126 X, y, w = table.to_numpy() 127 n, m = numpy.shape(X) 128 129 X, mu_x, sigma_x = standardize(X) 130 y, coef0 = center(y) 131 132 t = self.t 133 134 if self.s is not None: 135 beta_full, rss, _, _ = numpy.linalg.lstsq(X, y) 136 t = self.s * numpy.sum(numpy.abs(beta_full)) 137 print "t =", t 138 139 import scipy.optimize 140 141 # objective function to be minimized 142 objective = lambda beta: numpy.linalg.norm(y  numpy.dot(X, beta)) 143 # initial guess for the regression parameters 144 beta_init = numpy.random.random(m) 145 # constraints for the regression coefficients 146 cnstr = lambda beta: t  numpy.sum(numpy.abs(beta)) 147 # optimal solution 148 coefficients = scipy.optimize.fmin_cobyla(objective, beta_init, cnstr, iprint=0) 149 150 # set small coefficients to 0 151 def set_2_0(c): return c if abs(c) > self.tol else 0 152 coefficients = numpy.array(map(set_2_0, coefficients)) 153 coefficients /= sigma_x 118 data = self.impute_table(data) 119 domain = data.domain 120 # prepare numpy matrices 121 X, y, _ = data.to_numpy() 122 n, m = X.shape 123 coefficients = zeros(m) 124 std_errors = array([float('nan')] * m) 125 p_vals = array([float('nan')] * m) 126 # standardize y 127 coef0, sigma_y = y.mean(), y.std() + 1e6 128 y = (y  coef0) / sigma_y 129 # standardize X and remove constant vars 130 mu_x = X.mean(axis=0) 131 X = mu_x 132 sigma_x = X.std(axis=0) 133 nz = sigma_x != 0 134 X = ascontiguousarray(X[:, nz]) 135 sigma_x = sigma_x[nz] 136 X /= sigma_x 137 m = sum(nz) 138 139 # run optimization method 140 lipschitz = self.get_lipschitz(X) 141 l = 0.5 * self.lasso_lambda * n / m 142 coefficients[nz] = self.fista(X, y, l, lipschitz) 143 coefficients[nz] *= sigma_y / sigma_x 144 145 d = dict(self.__dict__) 146 d.update({'n_boot': 0, 'n_perm': 0}) 154 147 155 148 # bootstrap estimator of standard error of the coefficient estimators 156 # assumption: fixed t149 # assumption: fixed regularization parameter 157 150 if self.n_boot > 0: 158 151 coeff_b = [] # bootstrapped coefficients 159 152 for i in range(self.n_boot): 160 tmp_ table = get_bootstrap_sample(table)161 l = LassoRegressionLearner( t=t, n_boot=0, n_perm=0)162 c = l(tmp_ table)153 tmp_data = get_bootstrap_sample(data) 154 l = LassoRegressionLearner(**d) 155 c = l(tmp_data) 163 156 coeff_b.append(c.coefficients) 164 std_errors_fixed_t = numpy.std(coeff_b, axis=0) 165 else: 166 std_errors_fixed_t = [float("nan")] * m 167 168 # permutation test to obtain the significance of the regression 169 #coefficients 157 std_errors[nz] = std(coeff_b, axis=0) 158 159 # permutation test to obtain the significance of 160 # the regression coefficients 170 161 if self.n_perm > 0: 171 162 coeff_p = [] 172 163 for i in range(self.n_perm): 173 tmp_ table = permute_responses(table)174 l = LassoRegressionLearner( t=t, n_boot=0, n_perm=0)175 c = l(tmp_ table)164 tmp_data = permute_responses(data) 165 l = LassoRegressionLearner(**d) 166 c = l(tmp_data) 176 167 coeff_p.append(c.coefficients) 177 p_vals = \ 178 numpy.sum(abs(numpy.array(coeff_p)) > \ 179 abs(numpy.array(coefficients)), \ 180 axis=0) / float(self.n_perm) 181 else: 182 p_vals = [float("nan")] * m 168 p_vals[nz] = (abs(coeff_p) > abs(coefficients)).sum(axis=0) 169 p_vals[nz] /= float(self.n_perm) 183 170 184 171 # dictionary of regression coefficients with standard errors 185 172 # and pvalues 186 dict_model = {}173 model = {} 187 174 for i, var in enumerate(domain.attributes): 188 dict_model[var.name] = (coefficients[i], std_errors_fixed_t[i], p_vals[i])175 model[var.name] = (coefficients[i], std_errors[i], p_vals[i]) 189 176 190 177 return LassoRegression(domain=domain, class_var=domain.class_var, 191 coef0=coef0, coefficients=coefficients, 192 std_errors_fixed_t=std_errors_fixed_t, 193 p_vals=p_vals, 194 dict_model=dict_model, 195 mu_x=mu_x) 178 coef0=coef0, coefficients=coefficients, std_errors=std_errors, 179 p_vals=p_vals, model=model, mu_x=mu_x) 196 180 197 181 deprecated_members({"nBoot": "n_boot", … … 201 185 202 186 class LassoRegression(Orange.classification.Classifier): 203 """Lasso regression predicts value of the response variable187 """Lasso regression predicts the value of the response variable 204 188 based on the values of independent variables. 205 189 … … 210 194 .. attribute:: coefficients 211 195 212 Regression coefficients , sotred in list.213 214 .. attribute:: std_errors _fixed_t215 216 Standard errors of the coefficient estimator for thefixed217 tuning parameter t. The standard errors are estimated using218 bootstrapping method.196 Regression coefficients. 197 198 .. attribute:: std_errors 199 200 Standard errors of coefficient estimates for a fixed 201 regularization parameter. The standard errors are estimated 202 using the bootstrapping method. 219 203 220 204 .. attribute:: p_vals 221 205 222 List of pvalues for the null hypothes is that the regression223 coefficients equal 0 based on nonparametric permutation test.224 225 .. attribute:: dict_model226 227 Statistical properties of the model stored in dictionary:206 List of pvalues for the null hypotheses that the regression 207 coefficients equal 0 based on a nonparametric permutation test. 208 209 .. attribute:: model 210 211 Dictionary with the statistical properties of the model: 228 212 Keys  names of the independent variables 229 213 Values  tuples (coefficient, standard error, pvalue) … … 231 215 .. attribute:: mu_x 232 216 233 Sample mean of the allindependent variables.217 Sample mean of independent variables. 234 218 235 219 """ 236 220 def __init__(self, domain=None, class_var=None, coef0=None, 237 coefficients=None, std_errors _fixed_t=None, p_vals=None,238 dict_model=None, mu_x=None):221 coefficients=None, std_errors=None, p_vals=None, 222 model=None, mu_x=None): 239 223 self.domain = domain 240 224 self.class_var = class_var 241 225 self.coef0 = coef0 242 226 self.coefficients = coefficients 243 self.std_errors _fixed_t = std_errors_fixed_t227 self.std_errors = std_errors 244 228 self.p_vals = p_vals 245 self. dict_model = dict_model229 self.model = model 246 230 self.mu_x = mu_x 231 232 def _miss_2_0(self, x): 233 return x if x != '?' else 0 247 234 248 235 @deprecated_keywords({"resultType": "result_type"}) 249 236 def __call__(self, instance, result_type=Orange.core.GetValue): 250 237 """ 251 :param instance: data instance for which the value of the response252 variable will be predicted253 :type instance: 238 :param instance: Data instance for which the value of the response 239 variable will be predicted. 240 :type instance: :obj:`Orange.data.Instance` 254 241 """ 255 242 ins = Orange.data.Instance(self.domain, instance) 256 if "?" in ins: # missing value > corresponding coefficient omitted 257 def miss_2_0(x): return x if x != "?" else 0 258 ins = map(miss_2_0, ins) 259 ins = numpy.array(ins)[:1]  self.mu_x 243 if '?' in ins: # missing value > corresponding coefficient omitted 244 ins = map(self._miss_2_0, ins) 245 ins = array(ins)[:1]  self.mu_x 260 246 else: 261 ins = numpy.array(ins.native())[:1]  self.mu_x262 263 y_hat = numpy.dot(self.coefficients, ins) + self.coef0247 ins = array(ins.native())[:1]  self.mu_x 248 249 y_hat = dot(self.coefficients, ins) + self.coef0 264 250 y_hat = self.class_var(y_hat) 265 251 dist = Orange.statistics.distribution.Continuous(self.class_var) 266 dist[y_hat] = 1. 0252 dist[y_hat] = 1. 267 253 if result_type == Orange.core.GetValue: 268 254 return y_hat … … 274 260 @deprecated_keywords({"skipZero": "skip_zero"}) 275 261 def to_string(self, skip_zero=True): 276 """Prettyprints Lasso regression model,262 """Prettyprints a lasso regression model, 277 263 i.e. estimated regression coefficients with standard errors 278 and significances. Standard errors are obtained using bootstrapping 279 method and significances by the permuation test 280 281 :param skip_zero: if True variables with estimated coefficient equal to 0 282 are omitted 283 :type skip_zero: boolean 284 """ 285 286 from string import join 264 and significances. Standard errors are obtained using the 265 bootstrapping method and significances by a permuation test. 266 267 :param skip_zero: If True, variables with estimated coefficient 268 equal to 0 are omitted. 269 :type skip_zero: bool 270 """ 287 271 labels = ('Variable', 'Coeff Est', 'Std Error', 'p') 288 lines = [ join(['%10s' % l for l in labels], ' ')]289 290 fmt = "%10s " + join(["%10.3f"] * 3, " ") + " %5s"291 fmt1 = "%10s %10.3f"272 lines = [' '.join(['%10s' % l for l in labels])] 273 274 fmt = '%10s ' + ' '.join(['%10.3f'] * 3) + ' %5s' 275 fmt1 = '%10s %10.3f' 292 276 293 277 def get_star(p): 294 if p < 0.001: return "*"*3295 elif p < 0.01: return "*"*2296 elif p < 0.05: return "*"297 elif p < 0.1: return "."298 else: return " "278 if p < 0.001: return '*' * 3 279 elif p < 0.01: return '*' * 2 280 elif p < 0.05: return '*' 281 elif p < 0.1: return '.' 282 else: return ' ' 299 283 300 284 stars = get_star(self.p_vals[0]) … … 307 291 stars = get_star(self.p_vals[i]) 308 292 lines.append(fmt % (self.domain.attributes[i].name, 309 self.coefficients[i], self.std_errors _fixed_t[i],293 self.coefficients[i], self.std_errors[i], 310 294 self.p_vals[i], stars)) 311 lines.append("Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 empty 1") 312 lines.append("\n") 295 lines.append('Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 empty 1\n') 313 296 if skip_zero: 314 297 k = len(skipped) 315 298 if k == 0: 316 lines.append( "All variables have nonzero regression coefficients. ")299 lines.append('All variables have nonzero regression coefficients.') 317 300 else: 318 suff = "s" if k > 1 else "" 319 lines.append("For %d variable%s the regression coefficient equals 0: " \ 320 % (k, suff)) 321 for var in skipped: 322 lines.append(var) 323 return "\n".join(lines) 301 lines.append('For %d variable%s the regression coefficient equals 0:' 302 % (k, 's' if k > 1 else '')) 303 lines.append(', '.join(var for var in skipped)) 304 return '\n'.join(lines) 324 305 325 306 def __str__(self): … … 327 308 328 309 deprecated_members({"muX": "mu_x", 329 "stdErrorsFixedT": "std_errors _fixed_t",310 "stdErrorsFixedT": "std_errors", 330 311 "pVals": "p_vals", 331 "dictModel": " dict_model"},312 "dictModel": "model"}, 332 313 wrap_methods=["__init__"], 333 314 in_place=True)(LassoRegression) 334 315 335 316 if __name__ == "__main__": 336 337 import Orange 338 339 table = Orange.data.Table("housing.tab") 340 341 c = LassoRegressionLearner(table, t=len(table.domain)) 317 data = Orange.data.Table('housing') 318 c = LassoRegressionLearner(data) 342 319 print c
Note: See TracChangeset
for help on using the changeset viewer.