[4] | 1 | /* |
---|
| 2 | * This program is free software; you can redistribute it and/or modify |
---|
| 3 | * it under the terms of the GNU General Public License as published by |
---|
| 4 | * the Free Software Foundation; either version 2 of the License, or |
---|
| 5 | * (at your option) any later version. |
---|
| 6 | * |
---|
| 7 | * This program is distributed in the hope that it will be useful, |
---|
| 8 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
| 9 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
| 10 | * GNU General Public License for more details. |
---|
| 11 | * |
---|
| 12 | * You should have received a copy of the GNU General Public License |
---|
| 13 | * along with this program; if not, write to the Free Software |
---|
| 14 | * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. |
---|
| 15 | */ |
---|
| 16 | |
---|
| 17 | /* |
---|
| 18 | * Logistic.java |
---|
| 19 | * Copyright (C) 2003 University of Waikato, Hamilton, New Zealand |
---|
| 20 | * |
---|
| 21 | */ |
---|
| 22 | |
---|
| 23 | package weka.classifiers.functions; |
---|
| 24 | |
---|
| 25 | import weka.classifiers.Classifier; |
---|
| 26 | import weka.classifiers.AbstractClassifier; |
---|
| 27 | import weka.core.Capabilities; |
---|
| 28 | import weka.core.Instance; |
---|
| 29 | import weka.core.Instances; |
---|
| 30 | import weka.core.Optimization; |
---|
| 31 | import weka.core.Option; |
---|
| 32 | import weka.core.OptionHandler; |
---|
| 33 | import weka.core.RevisionUtils; |
---|
| 34 | import weka.core.TechnicalInformation; |
---|
| 35 | import weka.core.TechnicalInformationHandler; |
---|
| 36 | import weka.core.Utils; |
---|
| 37 | import weka.core.WeightedInstancesHandler; |
---|
| 38 | import weka.core.Capabilities.Capability; |
---|
| 39 | import weka.core.TechnicalInformation.Field; |
---|
| 40 | import weka.core.TechnicalInformation.Type; |
---|
| 41 | import weka.filters.Filter; |
---|
| 42 | import weka.filters.unsupervised.attribute.NominalToBinary; |
---|
| 43 | import weka.filters.unsupervised.attribute.RemoveUseless; |
---|
| 44 | import weka.filters.unsupervised.attribute.ReplaceMissingValues; |
---|
| 45 | |
---|
| 46 | import java.util.Enumeration; |
---|
| 47 | import java.util.Vector; |
---|
| 48 | |
---|
| 49 | /** |
---|
| 50 | <!-- globalinfo-start --> |
---|
| 51 | * Class for building and using a multinomial logistic regression model with a ridge estimator.<br/> |
---|
| 52 | * <br/> |
---|
| 53 | * There are some modifications, however, compared to the paper of leCessie and van Houwelingen(1992): <br/> |
---|
| 54 | * <br/> |
---|
| 55 | * If there are k classes for n instances with m attributes, the parameter matrix B to be calculated will be an m*(k-1) matrix.<br/> |
---|
| 56 | * <br/> |
---|
| 57 | * The probability for class j with the exception of the last class is<br/> |
---|
| 58 | * <br/> |
---|
| 59 | * Pj(Xi) = exp(XiBj)/((sum[j=1..(k-1)]exp(Xi*Bj))+1) <br/> |
---|
| 60 | * <br/> |
---|
| 61 | * The last class has probability<br/> |
---|
| 62 | * <br/> |
---|
| 63 | * 1-(sum[j=1..(k-1)]Pj(Xi)) <br/> |
---|
| 64 | * = 1/((sum[j=1..(k-1)]exp(Xi*Bj))+1)<br/> |
---|
| 65 | * <br/> |
---|
| 66 | * The (negative) multinomial log-likelihood is thus: <br/> |
---|
| 67 | * <br/> |
---|
| 68 | * L = -sum[i=1..n]{<br/> |
---|
| 69 | * sum[j=1..(k-1)](Yij * ln(Pj(Xi)))<br/> |
---|
| 70 | * +(1 - (sum[j=1..(k-1)]Yij)) <br/> |
---|
| 71 | * * ln(1 - sum[j=1..(k-1)]Pj(Xi))<br/> |
---|
| 72 | * } + ridge * (B^2)<br/> |
---|
| 73 | * <br/> |
---|
| 74 | * In order to find the matrix B for which L is minimised, a Quasi-Newton Method is used to search for the optimized values of the m*(k-1) variables. Note that before we use the optimization procedure, we 'squeeze' the matrix B into a m*(k-1) vector. For details of the optimization procedure, please check weka.core.Optimization class.<br/> |
---|
| 75 | * <br/> |
---|
| 76 | * Although original Logistic Regression does not deal with instance weights, we modify the algorithm a little bit to handle the instance weights.<br/> |
---|
| 77 | * <br/> |
---|
| 78 | * For more information see:<br/> |
---|
| 79 | * <br/> |
---|
| 80 | * le Cessie, S., van Houwelingen, J.C. (1992). Ridge Estimators in Logistic Regression. Applied Statistics. 41(1):191-201.<br/> |
---|
| 81 | * <br/> |
---|
| 82 | * Note: Missing values are replaced using a ReplaceMissingValuesFilter, and nominal attributes are transformed into numeric attributes using a NominalToBinaryFilter. |
---|
| 83 | * <p/> |
---|
| 84 | <!-- globalinfo-end --> |
---|
| 85 | * |
---|
| 86 | <!-- technical-bibtex-start --> |
---|
| 87 | * BibTeX: |
---|
| 88 | * <pre> |
---|
| 89 | * @article{leCessie1992, |
---|
| 90 | * author = {le Cessie, S. and van Houwelingen, J.C.}, |
---|
| 91 | * journal = {Applied Statistics}, |
---|
| 92 | * number = {1}, |
---|
| 93 | * pages = {191-201}, |
---|
| 94 | * title = {Ridge Estimators in Logistic Regression}, |
---|
| 95 | * volume = {41}, |
---|
| 96 | * year = {1992} |
---|
| 97 | * } |
---|
| 98 | * </pre> |
---|
| 99 | * <p/> |
---|
| 100 | <!-- technical-bibtex-end --> |
---|
| 101 | * |
---|
| 102 | <!-- options-start --> |
---|
| 103 | * Valid options are: <p/> |
---|
| 104 | * |
---|
| 105 | * <pre> -D |
---|
| 106 | * Turn on debugging output.</pre> |
---|
| 107 | * |
---|
| 108 | * <pre> -R <ridge> |
---|
| 109 | * Set the ridge in the log-likelihood.</pre> |
---|
| 110 | * |
---|
| 111 | * <pre> -M <number> |
---|
| 112 | * Set the maximum number of iterations (default -1, until convergence).</pre> |
---|
| 113 | * |
---|
| 114 | <!-- options-end --> |
---|
| 115 | * |
---|
| 116 | * @author Xin Xu (xx5@cs.waikato.ac.nz) |
---|
| 117 | * @version $Revision: 5928 $ |
---|
| 118 | */ |
---|
| 119 | public class Logistic extends AbstractClassifier |
---|
| 120 | implements OptionHandler, WeightedInstancesHandler, TechnicalInformationHandler { |
---|
| 121 | |
---|
| 122 | /** for serialization */ |
---|
| 123 | static final long serialVersionUID = 3932117032546553727L; |
---|
| 124 | |
---|
| 125 | /** The coefficients (optimized parameters) of the model */ |
---|
| 126 | protected double [][] m_Par; |
---|
| 127 | |
---|
| 128 | /** The data saved as a matrix */ |
---|
| 129 | protected double [][] m_Data; |
---|
| 130 | |
---|
| 131 | /** The number of attributes in the model */ |
---|
| 132 | protected int m_NumPredictors; |
---|
| 133 | |
---|
| 134 | /** The index of the class attribute */ |
---|
| 135 | protected int m_ClassIndex; |
---|
| 136 | |
---|
| 137 | /** The number of the class labels */ |
---|
| 138 | protected int m_NumClasses; |
---|
| 139 | |
---|
| 140 | /** The ridge parameter. */ |
---|
| 141 | protected double m_Ridge = 1e-8; |
---|
| 142 | |
---|
| 143 | /** An attribute filter */ |
---|
| 144 | private RemoveUseless m_AttFilter; |
---|
| 145 | |
---|
| 146 | /** The filter used to make attributes numeric. */ |
---|
| 147 | private NominalToBinary m_NominalToBinary; |
---|
| 148 | |
---|
| 149 | /** The filter used to get rid of missing values. */ |
---|
| 150 | private ReplaceMissingValues m_ReplaceMissingValues; |
---|
| 151 | |
---|
| 152 | /** Debugging output */ |
---|
| 153 | protected boolean m_Debug; |
---|
| 154 | |
---|
| 155 | /** Log-likelihood of the searched model */ |
---|
| 156 | protected double m_LL; |
---|
| 157 | |
---|
| 158 | /** The maximum number of iterations. */ |
---|
| 159 | private int m_MaxIts = -1; |
---|
| 160 | |
---|
| 161 | private Instances m_structure; |
---|
| 162 | |
---|
| 163 | /** |
---|
| 164 | * Returns a string describing this classifier |
---|
| 165 | * @return a description of the classifier suitable for |
---|
| 166 | * displaying in the explorer/experimenter gui |
---|
| 167 | */ |
---|
| 168 | public String globalInfo() { |
---|
| 169 | return "Class for building and using a multinomial logistic " |
---|
| 170 | +"regression model with a ridge estimator.\n\n" |
---|
| 171 | +"There are some modifications, however, compared to the paper of " |
---|
| 172 | +"leCessie and van Houwelingen(1992): \n\n" |
---|
| 173 | +"If there are k classes for n instances with m attributes, the " |
---|
| 174 | +"parameter matrix B to be calculated will be an m*(k-1) matrix.\n\n" |
---|
| 175 | +"The probability for class j with the exception of the last class is\n\n" |
---|
| 176 | +"Pj(Xi) = exp(XiBj)/((sum[j=1..(k-1)]exp(Xi*Bj))+1) \n\n" |
---|
| 177 | +"The last class has probability\n\n" |
---|
| 178 | +"1-(sum[j=1..(k-1)]Pj(Xi)) \n\t= 1/((sum[j=1..(k-1)]exp(Xi*Bj))+1)\n\n" |
---|
| 179 | +"The (negative) multinomial log-likelihood is thus: \n\n" |
---|
| 180 | +"L = -sum[i=1..n]{\n\tsum[j=1..(k-1)](Yij * ln(Pj(Xi)))" |
---|
| 181 | +"\n\t+(1 - (sum[j=1..(k-1)]Yij)) \n\t* ln(1 - sum[j=1..(k-1)]Pj(Xi))" |
---|
| 182 | +"\n\t} + ridge * (B^2)\n\n" |
---|
| 183 | +"In order to find the matrix B for which L is minimised, a " |
---|
| 184 | +"Quasi-Newton Method is used to search for the optimized values of " |
---|
| 185 | +"the m*(k-1) variables. Note that before we use the optimization " |
---|
| 186 | +"procedure, we 'squeeze' the matrix B into a m*(k-1) vector. For " |
---|
| 187 | +"details of the optimization procedure, please check " |
---|
| 188 | +"weka.core.Optimization class.\n\n" |
---|
| 189 | +"Although original Logistic Regression does not deal with instance " |
---|
| 190 | +"weights, we modify the algorithm a little bit to handle the " |
---|
| 191 | +"instance weights.\n\n" |
---|
| 192 | +"For more information see:\n\n" |
---|
| 193 | + getTechnicalInformation().toString() + "\n\n" |
---|
| 194 | +"Note: Missing values are replaced using a ReplaceMissingValuesFilter, and " |
---|
| 195 | +"nominal attributes are transformed into numeric attributes using a " |
---|
| 196 | +"NominalToBinaryFilter."; |
---|
| 197 | } |
---|
| 198 | |
---|
| 199 | /** |
---|
| 200 | * Returns an instance of a TechnicalInformation object, containing |
---|
| 201 | * detailed information about the technical background of this class, |
---|
| 202 | * e.g., paper reference or book this class is based on. |
---|
| 203 | * |
---|
| 204 | * @return the technical information about this class |
---|
| 205 | */ |
---|
| 206 | public TechnicalInformation getTechnicalInformation() { |
---|
| 207 | TechnicalInformation result; |
---|
| 208 | |
---|
| 209 | result = new TechnicalInformation(Type.ARTICLE); |
---|
| 210 | result.setValue(Field.AUTHOR, "le Cessie, S. and van Houwelingen, J.C."); |
---|
| 211 | result.setValue(Field.YEAR, "1992"); |
---|
| 212 | result.setValue(Field.TITLE, "Ridge Estimators in Logistic Regression"); |
---|
| 213 | result.setValue(Field.JOURNAL, "Applied Statistics"); |
---|
| 214 | result.setValue(Field.VOLUME, "41"); |
---|
| 215 | result.setValue(Field.NUMBER, "1"); |
---|
| 216 | result.setValue(Field.PAGES, "191-201"); |
---|
| 217 | |
---|
| 218 | return result; |
---|
| 219 | } |
---|
| 220 | |
---|
| 221 | /** |
---|
| 222 | * Returns an enumeration describing the available options |
---|
| 223 | * |
---|
| 224 | * @return an enumeration of all the available options |
---|
| 225 | */ |
---|
| 226 | public Enumeration listOptions() { |
---|
| 227 | Vector newVector = new Vector(3); |
---|
| 228 | newVector.addElement(new Option("\tTurn on debugging output.", |
---|
| 229 | "D", 0, "-D")); |
---|
| 230 | newVector.addElement(new Option("\tSet the ridge in the log-likelihood.", |
---|
| 231 | "R", 1, "-R <ridge>")); |
---|
| 232 | newVector.addElement(new Option("\tSet the maximum number of iterations"+ |
---|
| 233 | " (default -1, until convergence).", |
---|
| 234 | "M", 1, "-M <number>")); |
---|
| 235 | return newVector.elements(); |
---|
| 236 | } |
---|
| 237 | |
---|
| 238 | /** |
---|
| 239 | * Parses a given list of options. <p/> |
---|
| 240 | * |
---|
| 241 | <!-- options-start --> |
---|
| 242 | * Valid options are: <p/> |
---|
| 243 | * |
---|
| 244 | * <pre> -D |
---|
| 245 | * Turn on debugging output.</pre> |
---|
| 246 | * |
---|
| 247 | * <pre> -R <ridge> |
---|
| 248 | * Set the ridge in the log-likelihood.</pre> |
---|
| 249 | * |
---|
| 250 | * <pre> -M <number> |
---|
| 251 | * Set the maximum number of iterations (default -1, until convergence).</pre> |
---|
| 252 | * |
---|
| 253 | <!-- options-end --> |
---|
| 254 | * |
---|
| 255 | * @param options the list of options as an array of strings |
---|
| 256 | * @throws Exception if an option is not supported |
---|
| 257 | */ |
---|
| 258 | public void setOptions(String[] options) throws Exception { |
---|
| 259 | setDebug(Utils.getFlag('D', options)); |
---|
| 260 | |
---|
| 261 | String ridgeString = Utils.getOption('R', options); |
---|
| 262 | if (ridgeString.length() != 0) |
---|
| 263 | m_Ridge = Double.parseDouble(ridgeString); |
---|
| 264 | else |
---|
| 265 | m_Ridge = 1.0e-8; |
---|
| 266 | |
---|
| 267 | String maxItsString = Utils.getOption('M', options); |
---|
| 268 | if (maxItsString.length() != 0) |
---|
| 269 | m_MaxIts = Integer.parseInt(maxItsString); |
---|
| 270 | else |
---|
| 271 | m_MaxIts = -1; |
---|
| 272 | } |
---|
| 273 | |
---|
| 274 | /** |
---|
| 275 | * Gets the current settings of the classifier. |
---|
| 276 | * |
---|
| 277 | * @return an array of strings suitable for passing to setOptions |
---|
| 278 | */ |
---|
| 279 | public String [] getOptions() { |
---|
| 280 | |
---|
| 281 | String [] options = new String [5]; |
---|
| 282 | int current = 0; |
---|
| 283 | |
---|
| 284 | if (getDebug()) |
---|
| 285 | options[current++] = "-D"; |
---|
| 286 | options[current++] = "-R"; |
---|
| 287 | options[current++] = ""+m_Ridge; |
---|
| 288 | options[current++] = "-M"; |
---|
| 289 | options[current++] = ""+m_MaxIts; |
---|
| 290 | while (current < options.length) |
---|
| 291 | options[current++] = ""; |
---|
| 292 | return options; |
---|
| 293 | } |
---|
| 294 | |
---|
| 295 | /** |
---|
| 296 | * Returns the tip text for this property |
---|
| 297 | * @return tip text for this property suitable for |
---|
| 298 | * displaying in the explorer/experimenter gui |
---|
| 299 | */ |
---|
| 300 | public String debugTipText() { |
---|
| 301 | return "Output debug information to the console."; |
---|
| 302 | } |
---|
| 303 | |
---|
| 304 | /** |
---|
| 305 | * Sets whether debugging output will be printed. |
---|
| 306 | * |
---|
| 307 | * @param debug true if debugging output should be printed |
---|
| 308 | */ |
---|
| 309 | public void setDebug(boolean debug) { |
---|
| 310 | m_Debug = debug; |
---|
| 311 | } |
---|
| 312 | |
---|
| 313 | /** |
---|
| 314 | * Gets whether debugging output will be printed. |
---|
| 315 | * |
---|
| 316 | * @return true if debugging output will be printed |
---|
| 317 | */ |
---|
| 318 | public boolean getDebug() { |
---|
| 319 | return m_Debug; |
---|
| 320 | } |
---|
| 321 | |
---|
| 322 | /** |
---|
| 323 | * Returns the tip text for this property |
---|
| 324 | * @return tip text for this property suitable for |
---|
| 325 | * displaying in the explorer/experimenter gui |
---|
| 326 | */ |
---|
| 327 | public String ridgeTipText() { |
---|
| 328 | return "Set the Ridge value in the log-likelihood."; |
---|
| 329 | } |
---|
| 330 | |
---|
| 331 | /** |
---|
| 332 | * Sets the ridge in the log-likelihood. |
---|
| 333 | * |
---|
| 334 | * @param ridge the ridge |
---|
| 335 | */ |
---|
| 336 | public void setRidge(double ridge) { |
---|
| 337 | m_Ridge = ridge; |
---|
| 338 | } |
---|
| 339 | |
---|
| 340 | /** |
---|
| 341 | * Gets the ridge in the log-likelihood. |
---|
| 342 | * |
---|
| 343 | * @return the ridge |
---|
| 344 | */ |
---|
| 345 | public double getRidge() { |
---|
| 346 | return m_Ridge; |
---|
| 347 | } |
---|
| 348 | |
---|
| 349 | /** |
---|
| 350 | * Returns the tip text for this property |
---|
| 351 | * @return tip text for this property suitable for |
---|
| 352 | * displaying in the explorer/experimenter gui |
---|
| 353 | */ |
---|
| 354 | public String maxItsTipText() { |
---|
| 355 | return "Maximum number of iterations to perform."; |
---|
| 356 | } |
---|
| 357 | |
---|
| 358 | /** |
---|
| 359 | * Get the value of MaxIts. |
---|
| 360 | * |
---|
| 361 | * @return Value of MaxIts. |
---|
| 362 | */ |
---|
| 363 | public int getMaxIts() { |
---|
| 364 | |
---|
| 365 | return m_MaxIts; |
---|
| 366 | } |
---|
| 367 | |
---|
| 368 | /** |
---|
| 369 | * Set the value of MaxIts. |
---|
| 370 | * |
---|
| 371 | * @param newMaxIts Value to assign to MaxIts. |
---|
| 372 | */ |
---|
| 373 | public void setMaxIts(int newMaxIts) { |
---|
| 374 | |
---|
| 375 | m_MaxIts = newMaxIts; |
---|
| 376 | } |
---|
| 377 | |
---|
| 378 | private class OptEng extends Optimization{ |
---|
| 379 | /** Weights of instances in the data */ |
---|
| 380 | private double[] weights; |
---|
| 381 | |
---|
| 382 | /** Class labels of instances */ |
---|
| 383 | private int[] cls; |
---|
| 384 | |
---|
| 385 | /** |
---|
| 386 | * Set the weights of instances |
---|
| 387 | * @param w the weights to be set |
---|
| 388 | */ |
---|
| 389 | public void setWeights(double[] w) { |
---|
| 390 | weights = w; |
---|
| 391 | } |
---|
| 392 | |
---|
| 393 | /** |
---|
| 394 | * Set the class labels of instances |
---|
| 395 | * @param c the class labels to be set |
---|
| 396 | */ |
---|
| 397 | public void setClassLabels(int[] c) { |
---|
| 398 | cls = c; |
---|
| 399 | } |
---|
| 400 | |
---|
| 401 | /** |
---|
| 402 | * Evaluate objective function |
---|
| 403 | * @param x the current values of variables |
---|
| 404 | * @return the value of the objective function |
---|
| 405 | */ |
---|
| 406 | protected double objectiveFunction(double[] x){ |
---|
| 407 | double nll = 0; // -LogLikelihood |
---|
| 408 | int dim = m_NumPredictors+1; // Number of variables per class |
---|
| 409 | |
---|
| 410 | for(int i=0; i<cls.length; i++){ // ith instance |
---|
| 411 | |
---|
| 412 | double[] exp = new double[m_NumClasses-1]; |
---|
| 413 | int index; |
---|
| 414 | for(int offset=0; offset<m_NumClasses-1; offset++){ |
---|
| 415 | index = offset * dim; |
---|
| 416 | for(int j=0; j<dim; j++) |
---|
| 417 | exp[offset] += m_Data[i][j]*x[index + j]; |
---|
| 418 | } |
---|
| 419 | double max = exp[Utils.maxIndex(exp)]; |
---|
| 420 | double denom = Math.exp(-max); |
---|
| 421 | double num; |
---|
| 422 | if (cls[i] == m_NumClasses - 1) { // Class of this instance |
---|
| 423 | num = -max; |
---|
| 424 | } else { |
---|
| 425 | num = exp[cls[i]] - max; |
---|
| 426 | } |
---|
| 427 | for(int offset=0; offset<m_NumClasses-1; offset++){ |
---|
| 428 | denom += Math.exp(exp[offset] - max); |
---|
| 429 | } |
---|
| 430 | |
---|
| 431 | nll -= weights[i]*(num - Math.log(denom)); // Weighted NLL |
---|
| 432 | } |
---|
| 433 | |
---|
| 434 | // Ridge: note that intercepts NOT included |
---|
| 435 | for(int offset=0; offset<m_NumClasses-1; offset++){ |
---|
| 436 | for(int r=1; r<dim; r++) |
---|
| 437 | nll += m_Ridge*x[offset*dim+r]*x[offset*dim+r]; |
---|
| 438 | } |
---|
| 439 | |
---|
| 440 | return nll; |
---|
| 441 | } |
---|
| 442 | |
---|
| 443 | /** |
---|
| 444 | * Evaluate Jacobian vector |
---|
| 445 | * @param x the current values of variables |
---|
| 446 | * @return the gradient vector |
---|
| 447 | */ |
---|
| 448 | protected double[] evaluateGradient(double[] x){ |
---|
| 449 | double[] grad = new double[x.length]; |
---|
| 450 | int dim = m_NumPredictors+1; // Number of variables per class |
---|
| 451 | |
---|
| 452 | for(int i=0; i<cls.length; i++){ // ith instance |
---|
| 453 | double[] num=new double[m_NumClasses-1]; // numerator of [-log(1+sum(exp))]' |
---|
| 454 | int index; |
---|
| 455 | for(int offset=0; offset<m_NumClasses-1; offset++){ // Which part of x |
---|
| 456 | double exp=0.0; |
---|
| 457 | index = offset * dim; |
---|
| 458 | for(int j=0; j<dim; j++) |
---|
| 459 | exp += m_Data[i][j]*x[index + j]; |
---|
| 460 | num[offset] = exp; |
---|
| 461 | } |
---|
| 462 | |
---|
| 463 | double max = num[Utils.maxIndex(num)]; |
---|
| 464 | double denom = Math.exp(-max); // Denominator of [-log(1+sum(exp))]' |
---|
| 465 | for(int offset=0; offset<m_NumClasses-1; offset++){ |
---|
| 466 | num[offset] = Math.exp(num[offset] - max); |
---|
| 467 | denom += num[offset]; |
---|
| 468 | } |
---|
| 469 | Utils.normalize(num, denom); |
---|
| 470 | |
---|
| 471 | // Update denominator of the gradient of -log(Posterior) |
---|
| 472 | double firstTerm; |
---|
| 473 | for(int offset=0; offset<m_NumClasses-1; offset++){ // Which part of x |
---|
| 474 | index = offset * dim; |
---|
| 475 | firstTerm = weights[i] * num[offset]; |
---|
| 476 | for(int q=0; q<dim; q++){ |
---|
| 477 | grad[index + q] += firstTerm * m_Data[i][q]; |
---|
| 478 | } |
---|
| 479 | } |
---|
| 480 | |
---|
| 481 | if(cls[i] != m_NumClasses-1){ // Not the last class |
---|
| 482 | for(int p=0; p<dim; p++){ |
---|
| 483 | grad[cls[i]*dim+p] -= weights[i]*m_Data[i][p]; |
---|
| 484 | } |
---|
| 485 | } |
---|
| 486 | } |
---|
| 487 | |
---|
| 488 | // Ridge: note that intercepts NOT included |
---|
| 489 | for(int offset=0; offset<m_NumClasses-1; offset++){ |
---|
| 490 | for(int r=1; r<dim; r++) |
---|
| 491 | grad[offset*dim+r] += 2*m_Ridge*x[offset*dim+r]; |
---|
| 492 | } |
---|
| 493 | |
---|
| 494 | return grad; |
---|
| 495 | } |
---|
| 496 | |
---|
| 497 | /** |
---|
| 498 | * Returns the revision string. |
---|
| 499 | * |
---|
| 500 | * @return the revision |
---|
| 501 | */ |
---|
| 502 | public String getRevision() { |
---|
| 503 | return RevisionUtils.extract("$Revision: 5928 $"); |
---|
| 504 | } |
---|
| 505 | } |
---|
| 506 | |
---|
| 507 | /** |
---|
| 508 | * Returns default capabilities of the classifier. |
---|
| 509 | * |
---|
| 510 | * @return the capabilities of this classifier |
---|
| 511 | */ |
---|
| 512 | public Capabilities getCapabilities() { |
---|
| 513 | Capabilities result = super.getCapabilities(); |
---|
| 514 | result.disableAll(); |
---|
| 515 | |
---|
| 516 | // attributes |
---|
| 517 | result.enable(Capability.NOMINAL_ATTRIBUTES); |
---|
| 518 | result.enable(Capability.NUMERIC_ATTRIBUTES); |
---|
| 519 | result.enable(Capability.DATE_ATTRIBUTES); |
---|
| 520 | result.enable(Capability.MISSING_VALUES); |
---|
| 521 | |
---|
| 522 | // class |
---|
| 523 | result.enable(Capability.NOMINAL_CLASS); |
---|
| 524 | result.enable(Capability.MISSING_CLASS_VALUES); |
---|
| 525 | |
---|
| 526 | return result; |
---|
| 527 | } |
---|
| 528 | |
---|
| 529 | /** |
---|
| 530 | * Builds the classifier |
---|
| 531 | * |
---|
| 532 | * @param train the training data to be used for generating the |
---|
| 533 | * boosted classifier. |
---|
| 534 | * @throws Exception if the classifier could not be built successfully |
---|
| 535 | */ |
---|
| 536 | public void buildClassifier(Instances train) throws Exception { |
---|
| 537 | // can classifier handle the data? |
---|
| 538 | getCapabilities().testWithFail(train); |
---|
| 539 | |
---|
| 540 | // remove instances with missing class |
---|
| 541 | train = new Instances(train); |
---|
| 542 | train.deleteWithMissingClass(); |
---|
| 543 | |
---|
| 544 | // Replace missing values |
---|
| 545 | m_ReplaceMissingValues = new ReplaceMissingValues(); |
---|
| 546 | m_ReplaceMissingValues.setInputFormat(train); |
---|
| 547 | train = Filter.useFilter(train, m_ReplaceMissingValues); |
---|
| 548 | |
---|
| 549 | // Remove useless attributes |
---|
| 550 | m_AttFilter = new RemoveUseless(); |
---|
| 551 | m_AttFilter.setInputFormat(train); |
---|
| 552 | train = Filter.useFilter(train, m_AttFilter); |
---|
| 553 | |
---|
| 554 | // Transform attributes |
---|
| 555 | m_NominalToBinary = new NominalToBinary(); |
---|
| 556 | m_NominalToBinary.setInputFormat(train); |
---|
| 557 | train = Filter.useFilter(train, m_NominalToBinary); |
---|
| 558 | |
---|
| 559 | // Save the structure for printing the model |
---|
| 560 | m_structure = new Instances(train, 0); |
---|
| 561 | |
---|
| 562 | // Extract data |
---|
| 563 | m_ClassIndex = train.classIndex(); |
---|
| 564 | m_NumClasses = train.numClasses(); |
---|
| 565 | |
---|
| 566 | int nK = m_NumClasses - 1; // Only K-1 class labels needed |
---|
| 567 | int nR = m_NumPredictors = train.numAttributes() - 1; |
---|
| 568 | int nC = train.numInstances(); |
---|
| 569 | |
---|
| 570 | m_Data = new double[nC][nR + 1]; // Data values |
---|
| 571 | int [] Y = new int[nC]; // Class labels |
---|
| 572 | double [] xMean= new double[nR + 1]; // Attribute means |
---|
| 573 | double [] xSD = new double[nR + 1]; // Attribute stddev's |
---|
| 574 | double [] sY = new double[nK + 1]; // Number of classes |
---|
| 575 | double [] weights = new double[nC]; // Weights of instances |
---|
| 576 | double totWeights = 0; // Total weights of the instances |
---|
| 577 | m_Par = new double[nR + 1][nK]; // Optimized parameter values |
---|
| 578 | |
---|
| 579 | if (m_Debug) { |
---|
| 580 | System.out.println("Extracting data..."); |
---|
| 581 | } |
---|
| 582 | |
---|
| 583 | for (int i = 0; i < nC; i++) { |
---|
| 584 | // initialize X[][] |
---|
| 585 | Instance current = train.instance(i); |
---|
| 586 | Y[i] = (int)current.classValue(); // Class value starts from 0 |
---|
| 587 | weights[i] = current.weight(); // Dealing with weights |
---|
| 588 | totWeights += weights[i]; |
---|
| 589 | |
---|
| 590 | m_Data[i][0] = 1; |
---|
| 591 | int j = 1; |
---|
| 592 | for (int k = 0; k <= nR; k++) { |
---|
| 593 | if (k != m_ClassIndex) { |
---|
| 594 | double x = current.value(k); |
---|
| 595 | m_Data[i][j] = x; |
---|
| 596 | xMean[j] += weights[i]*x; |
---|
| 597 | xSD[j] += weights[i]*x*x; |
---|
| 598 | j++; |
---|
| 599 | } |
---|
| 600 | } |
---|
| 601 | |
---|
| 602 | // Class count |
---|
| 603 | sY[Y[i]]++; |
---|
| 604 | } |
---|
| 605 | |
---|
| 606 | if((totWeights <= 1) && (nC > 1)) |
---|
| 607 | throw new Exception("Sum of weights of instances less than 1, please reweight!"); |
---|
| 608 | |
---|
| 609 | xMean[0] = 0; xSD[0] = 1; |
---|
| 610 | for (int j = 1; j <= nR; j++) { |
---|
| 611 | xMean[j] = xMean[j] / totWeights; |
---|
| 612 | if(totWeights > 1) |
---|
| 613 | xSD[j] = Math.sqrt(Math.abs(xSD[j] - totWeights*xMean[j]*xMean[j])/(totWeights-1)); |
---|
| 614 | else |
---|
| 615 | xSD[j] = 0; |
---|
| 616 | } |
---|
| 617 | |
---|
| 618 | if (m_Debug) { |
---|
| 619 | // Output stats about input data |
---|
| 620 | System.out.println("Descriptives..."); |
---|
| 621 | for (int m = 0; m <= nK; m++) |
---|
| 622 | System.out.println(sY[m] + " cases have class " + m); |
---|
| 623 | System.out.println("\n Variable Avg SD "); |
---|
| 624 | for (int j = 1; j <= nR; j++) |
---|
| 625 | System.out.println(Utils.doubleToString(j,8,4) |
---|
| 626 | + Utils.doubleToString(xMean[j], 10, 4) |
---|
| 627 | + Utils.doubleToString(xSD[j], 10, 4) |
---|
| 628 | ); |
---|
| 629 | } |
---|
| 630 | |
---|
| 631 | // Normalise input data |
---|
| 632 | for (int i = 0; i < nC; i++) { |
---|
| 633 | for (int j = 0; j <= nR; j++) { |
---|
| 634 | if (xSD[j] != 0) { |
---|
| 635 | m_Data[i][j] = (m_Data[i][j] - xMean[j]) / xSD[j]; |
---|
| 636 | } |
---|
| 637 | } |
---|
| 638 | } |
---|
| 639 | |
---|
| 640 | if (m_Debug) { |
---|
| 641 | System.out.println("\nIteration History..." ); |
---|
| 642 | } |
---|
| 643 | |
---|
| 644 | double x[] = new double[(nR+1)*nK]; |
---|
| 645 | double[][] b = new double[2][x.length]; // Boundary constraints, N/A here |
---|
| 646 | |
---|
| 647 | // Initialize |
---|
| 648 | for(int p=0; p<nK; p++){ |
---|
| 649 | int offset=p*(nR+1); |
---|
| 650 | x[offset] = Math.log(sY[p]+1.0) - Math.log(sY[nK]+1.0); // Null model |
---|
| 651 | b[0][offset] = Double.NaN; |
---|
| 652 | b[1][offset] = Double.NaN; |
---|
| 653 | for (int q=1; q <= nR; q++){ |
---|
| 654 | x[offset+q] = 0.0; |
---|
| 655 | b[0][offset+q] = Double.NaN; |
---|
| 656 | b[1][offset+q] = Double.NaN; |
---|
| 657 | } |
---|
| 658 | } |
---|
| 659 | |
---|
| 660 | OptEng opt = new OptEng(); |
---|
| 661 | opt.setDebug(m_Debug); |
---|
| 662 | opt.setWeights(weights); |
---|
| 663 | opt.setClassLabels(Y); |
---|
| 664 | |
---|
| 665 | if(m_MaxIts == -1){ // Search until convergence |
---|
| 666 | x = opt.findArgmin(x, b); |
---|
| 667 | while(x==null){ |
---|
| 668 | x = opt.getVarbValues(); |
---|
| 669 | if (m_Debug) |
---|
| 670 | System.out.println("200 iterations finished, not enough!"); |
---|
| 671 | x = opt.findArgmin(x, b); |
---|
| 672 | } |
---|
| 673 | if (m_Debug) |
---|
| 674 | System.out.println(" -------------<Converged>--------------"); |
---|
| 675 | } |
---|
| 676 | else{ |
---|
| 677 | opt.setMaxIteration(m_MaxIts); |
---|
| 678 | x = opt.findArgmin(x, b); |
---|
| 679 | if(x==null) // Not enough, but use the current value |
---|
| 680 | x = opt.getVarbValues(); |
---|
| 681 | } |
---|
| 682 | |
---|
| 683 | m_LL = -opt.getMinFunction(); // Log-likelihood |
---|
| 684 | |
---|
| 685 | // Don't need data matrix anymore |
---|
| 686 | m_Data = null; |
---|
| 687 | |
---|
| 688 | // Convert coefficients back to non-normalized attribute units |
---|
| 689 | for(int i=0; i < nK; i++){ |
---|
| 690 | m_Par[0][i] = x[i*(nR+1)]; |
---|
| 691 | for(int j = 1; j <= nR; j++) { |
---|
| 692 | m_Par[j][i] = x[i*(nR+1)+j]; |
---|
| 693 | if (xSD[j] != 0) { |
---|
| 694 | m_Par[j][i] /= xSD[j]; |
---|
| 695 | m_Par[0][i] -= m_Par[j][i] * xMean[j]; |
---|
| 696 | } |
---|
| 697 | } |
---|
| 698 | } |
---|
| 699 | } |
---|
| 700 | |
---|
| 701 | /** |
---|
| 702 | * Computes the distribution for a given instance |
---|
| 703 | * |
---|
| 704 | * @param instance the instance for which distribution is computed |
---|
| 705 | * @return the distribution |
---|
| 706 | * @throws Exception if the distribution can't be computed successfully |
---|
| 707 | */ |
---|
| 708 | public double [] distributionForInstance(Instance instance) |
---|
| 709 | throws Exception { |
---|
| 710 | |
---|
| 711 | m_ReplaceMissingValues.input(instance); |
---|
| 712 | instance = m_ReplaceMissingValues.output(); |
---|
| 713 | m_AttFilter.input(instance); |
---|
| 714 | instance = m_AttFilter.output(); |
---|
| 715 | m_NominalToBinary.input(instance); |
---|
| 716 | instance = m_NominalToBinary.output(); |
---|
| 717 | |
---|
| 718 | // Extract the predictor columns into an array |
---|
| 719 | double [] instDat = new double [m_NumPredictors + 1]; |
---|
| 720 | int j = 1; |
---|
| 721 | instDat[0] = 1; |
---|
| 722 | for (int k = 0; k <= m_NumPredictors; k++) { |
---|
| 723 | if (k != m_ClassIndex) { |
---|
| 724 | instDat[j++] = instance.value(k); |
---|
| 725 | } |
---|
| 726 | } |
---|
| 727 | |
---|
| 728 | double [] distribution = evaluateProbability(instDat); |
---|
| 729 | return distribution; |
---|
| 730 | } |
---|
| 731 | |
---|
| 732 | /** |
---|
| 733 | * Compute the posterior distribution using optimized parameter values |
---|
| 734 | * and the testing instance. |
---|
| 735 | * @param data the testing instance |
---|
| 736 | * @return the posterior probability distribution |
---|
| 737 | */ |
---|
| 738 | private double[] evaluateProbability(double[] data){ |
---|
| 739 | double[] prob = new double[m_NumClasses], |
---|
| 740 | v = new double[m_NumClasses]; |
---|
| 741 | |
---|
| 742 | // Log-posterior before normalizing |
---|
| 743 | for(int j = 0; j < m_NumClasses-1; j++){ |
---|
| 744 | for(int k = 0; k <= m_NumPredictors; k++){ |
---|
| 745 | v[j] += m_Par[k][j] * data[k]; |
---|
| 746 | } |
---|
| 747 | } |
---|
| 748 | v[m_NumClasses-1] = 0; |
---|
| 749 | |
---|
| 750 | // Do so to avoid scaling problems |
---|
| 751 | for(int m=0; m < m_NumClasses; m++){ |
---|
| 752 | double sum = 0; |
---|
| 753 | for(int n=0; n < m_NumClasses-1; n++) |
---|
| 754 | sum += Math.exp(v[n] - v[m]); |
---|
| 755 | prob[m] = 1 / (sum + Math.exp(-v[m])); |
---|
| 756 | } |
---|
| 757 | |
---|
| 758 | return prob; |
---|
| 759 | } |
---|
| 760 | |
---|
| 761 | /** |
---|
| 762 | * Returns the coefficients for this logistic model. |
---|
| 763 | * The first dimension indexes the attributes, and |
---|
| 764 | * the second the classes. |
---|
| 765 | * |
---|
| 766 | * @return the coefficients for this logistic model |
---|
| 767 | */ |
---|
| 768 | public double [][] coefficients() { |
---|
| 769 | return m_Par; |
---|
| 770 | } |
---|
| 771 | |
---|
| 772 | /** |
---|
| 773 | * Gets a string describing the classifier. |
---|
| 774 | * |
---|
| 775 | * @return a string describing the classifer built. |
---|
| 776 | */ |
---|
| 777 | public String toString() { |
---|
| 778 | StringBuffer temp = new StringBuffer(); |
---|
| 779 | |
---|
| 780 | String result = ""; |
---|
| 781 | temp.append("Logistic Regression with ridge parameter of " + m_Ridge); |
---|
| 782 | if (m_Par == null) { |
---|
| 783 | return result + ": No model built yet."; |
---|
| 784 | } |
---|
| 785 | |
---|
| 786 | // find longest attribute name |
---|
| 787 | int attLength = 0; |
---|
| 788 | for (int i = 0; i < m_structure.numAttributes(); i++) { |
---|
| 789 | if (i != m_structure.classIndex() && |
---|
| 790 | m_structure.attribute(i).name().length() > attLength) { |
---|
| 791 | attLength = m_structure.attribute(i).name().length(); |
---|
| 792 | } |
---|
| 793 | } |
---|
| 794 | |
---|
| 795 | if ("Intercept".length() > attLength) { |
---|
| 796 | attLength = "Intercept".length(); |
---|
| 797 | } |
---|
| 798 | |
---|
| 799 | if ("Variable".length() > attLength) { |
---|
| 800 | attLength = "Variable".length(); |
---|
| 801 | } |
---|
| 802 | attLength += 2; |
---|
| 803 | |
---|
| 804 | int colWidth = 0; |
---|
| 805 | // check length of class names |
---|
| 806 | for (int i = 0; i < m_structure.classAttribute().numValues() - 1; i++) { |
---|
| 807 | if (m_structure.classAttribute().value(i).length() > colWidth) { |
---|
| 808 | colWidth = m_structure.classAttribute().value(i).length(); |
---|
| 809 | } |
---|
| 810 | } |
---|
| 811 | |
---|
| 812 | // check against coefficients and odds ratios |
---|
| 813 | for (int j = 1; j <= m_NumPredictors; j++) { |
---|
| 814 | for (int k = 0; k < m_NumClasses - 1; k++) { |
---|
| 815 | if (Utils.doubleToString(m_Par[j][k], 12, 4).trim().length() > colWidth) { |
---|
| 816 | colWidth = Utils.doubleToString(m_Par[j][k], 12, 4).trim().length(); |
---|
| 817 | } |
---|
| 818 | double ORc = Math.exp(m_Par[j][k]); |
---|
| 819 | String t = " " + ((ORc > 1e10) ? "" + ORc : Utils.doubleToString(ORc, 12, 4)); |
---|
| 820 | if (t.trim().length() > colWidth) { |
---|
| 821 | colWidth = t.trim().length(); |
---|
| 822 | } |
---|
| 823 | } |
---|
| 824 | } |
---|
| 825 | |
---|
| 826 | if ("Class".length() > colWidth) { |
---|
| 827 | colWidth = "Class".length(); |
---|
| 828 | } |
---|
| 829 | colWidth += 2; |
---|
| 830 | |
---|
| 831 | |
---|
| 832 | temp.append("\nCoefficients...\n"); |
---|
| 833 | temp.append(Utils.padLeft(" ", attLength) + Utils.padLeft("Class", colWidth) + "\n"); |
---|
| 834 | temp.append(Utils.padRight("Variable", attLength)); |
---|
| 835 | |
---|
| 836 | for (int i = 0; i < m_NumClasses - 1; i++) { |
---|
| 837 | String className = m_structure.classAttribute().value(i); |
---|
| 838 | temp.append(Utils.padLeft(className, colWidth)); |
---|
| 839 | } |
---|
| 840 | temp.append("\n"); |
---|
| 841 | int separatorL = attLength + ((m_NumClasses - 1) * colWidth); |
---|
| 842 | for (int i = 0; i < separatorL; i++) { |
---|
| 843 | temp.append("="); |
---|
| 844 | } |
---|
| 845 | temp.append("\n"); |
---|
| 846 | |
---|
| 847 | int j = 1; |
---|
| 848 | for (int i = 0; i < m_structure.numAttributes(); i++) { |
---|
| 849 | if (i != m_structure.classIndex()) { |
---|
| 850 | temp.append(Utils.padRight(m_structure.attribute(i).name(), attLength)); |
---|
| 851 | for (int k = 0; k < m_NumClasses-1; k++) { |
---|
| 852 | temp.append(Utils.padLeft(Utils.doubleToString(m_Par[j][k], 12, 4).trim(), colWidth)); |
---|
| 853 | } |
---|
| 854 | temp.append("\n"); |
---|
| 855 | j++; |
---|
| 856 | } |
---|
| 857 | } |
---|
| 858 | |
---|
| 859 | temp.append(Utils.padRight("Intercept", attLength)); |
---|
| 860 | for (int k = 0; k < m_NumClasses-1; k++) { |
---|
| 861 | temp.append(Utils.padLeft(Utils.doubleToString(m_Par[0][k], 10, 4).trim(), colWidth)); |
---|
| 862 | } |
---|
| 863 | temp.append("\n"); |
---|
| 864 | |
---|
| 865 | temp.append("\n\nOdds Ratios...\n"); |
---|
| 866 | temp.append(Utils.padLeft(" ", attLength) + Utils.padLeft("Class", colWidth) + "\n"); |
---|
| 867 | temp.append(Utils.padRight("Variable", attLength)); |
---|
| 868 | |
---|
| 869 | for (int i = 0; i < m_NumClasses - 1; i++) { |
---|
| 870 | String className = m_structure.classAttribute().value(i); |
---|
| 871 | temp.append(Utils.padLeft(className, colWidth)); |
---|
| 872 | } |
---|
| 873 | temp.append("\n"); |
---|
| 874 | for (int i = 0; i < separatorL; i++) { |
---|
| 875 | temp.append("="); |
---|
| 876 | } |
---|
| 877 | temp.append("\n"); |
---|
| 878 | |
---|
| 879 | j = 1; |
---|
| 880 | for (int i = 0; i < m_structure.numAttributes(); i++) { |
---|
| 881 | if (i != m_structure.classIndex()) { |
---|
| 882 | temp.append(Utils.padRight(m_structure.attribute(i).name(), attLength)); |
---|
| 883 | for (int k = 0; k < m_NumClasses-1; k++) { |
---|
| 884 | double ORc = Math.exp(m_Par[j][k]); |
---|
| 885 | String ORs = " " + ((ORc > 1e10) ? "" + ORc : Utils.doubleToString(ORc, 12, 4)); |
---|
| 886 | temp.append(Utils.padLeft(ORs.trim(), colWidth)); |
---|
| 887 | } |
---|
| 888 | temp.append("\n"); |
---|
| 889 | j++; |
---|
| 890 | } |
---|
| 891 | } |
---|
| 892 | |
---|
| 893 | return temp.toString(); |
---|
| 894 | } |
---|
| 895 | |
---|
| 896 | /** |
---|
| 897 | * Returns the revision string. |
---|
| 898 | * |
---|
| 899 | * @return the revision |
---|
| 900 | */ |
---|
| 901 | public String getRevision() { |
---|
| 902 | return RevisionUtils.extract("$Revision: 5928 $"); |
---|
| 903 | } |
---|
| 904 | |
---|
| 905 | /** |
---|
| 906 | * Main method for testing this class. |
---|
| 907 | * |
---|
| 908 | * @param argv should contain the command line arguments to the |
---|
| 909 | * scheme (see Evaluation) |
---|
| 910 | */ |
---|
| 911 | public static void main(String [] argv) { |
---|
| 912 | runClassifier(new Logistic(), argv); |
---|
| 913 | } |
---|
| 914 | } |
---|