~ubuntu-branches/ubuntu/precise/weka/precise

« back to all changes in this revision

Viewing changes to weka/classifiers/mi/TLD.java

  • Committer: Bazaar Package Importer
  • Author(s): Soeren Sonnenburg
  • Date: 2008-02-24 09:18:45 UTC
  • Revision ID: james.westby@ubuntu.com-20080224091845-1l8zy6fm6xipbzsr
Tags: upstream-3.5.7+tut1
ImportĀ upstreamĀ versionĀ 3.5.7+tut1

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
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
 * TLD.java
 
19
 * Copyright (C) 2005 University of Waikato, Hamilton, New Zealand
 
20
 *
 
21
 */
 
22
 
 
23
package weka.classifiers.mi;
 
24
 
 
25
import weka.classifiers.RandomizableClassifier;
 
26
import weka.core.Capabilities;
 
27
import weka.core.Instance;
 
28
import weka.core.Instances;
 
29
import weka.core.MultiInstanceCapabilitiesHandler;
 
30
import weka.core.Optimization;
 
31
import weka.core.Option;
 
32
import weka.core.OptionHandler;
 
33
import weka.core.TechnicalInformation;
 
34
import weka.core.TechnicalInformationHandler;
 
35
import weka.core.Utils;
 
36
import weka.core.Capabilities.Capability;
 
37
import weka.core.TechnicalInformation.Field;
 
38
import weka.core.TechnicalInformation.Type;
 
39
 
 
40
import java.util.Enumeration;
 
41
import java.util.Random;
 
42
import java.util.Vector;
 
43
 
 
44
/** 
 
45
 <!-- globalinfo-start -->
 
46
 * Two-Level Distribution approach, changes the starting value of the searching algorithm, supplement the cut-off modification and check missing values.<br/>
 
47
 * <br/>
 
48
 * For more information see:<br/>
 
49
 * <br/>
 
50
 * Xin Xu (2003). Statistical learning in multiple instance problem. Hamilton, NZ.
 
51
 * <p/>
 
52
 <!-- globalinfo-end -->
 
53
 *
 
54
 <!-- technical-bibtex-start -->
 
55
 * BibTeX:
 
56
 * <pre>
 
57
 * &#64;mastersthesis{Xu2003,
 
58
 *    address = {Hamilton, NZ},
 
59
 *    author = {Xin Xu},
 
60
 *    note = {0657.594},
 
61
 *    school = {University of Waikato},
 
62
 *    title = {Statistical learning in multiple instance problem},
 
63
 *    year = {2003}
 
64
 * }
 
65
 * </pre>
 
66
 * <p/>
 
67
 <!-- technical-bibtex-end -->
 
68
 *
 
69
 <!-- options-start -->
 
70
 * Valid options are: <p/>
 
71
 * 
 
72
 * <pre> -C
 
73
 *  Set whether or not use empirical
 
74
 *  log-odds cut-off instead of 0</pre>
 
75
 * 
 
76
 * <pre> -R &lt;numOfRuns&gt;
 
77
 *  Set the number of multiple runs 
 
78
 *  needed for searching the MLE.</pre>
 
79
 * 
 
80
 * <pre> -S &lt;num&gt;
 
81
 *  Random number seed.
 
82
 *  (default 1)</pre>
 
83
 * 
 
84
 * <pre> -D
 
85
 *  If set, classifier is run in debug mode and
 
86
 *  may output additional info to the console</pre>
 
87
 * 
 
88
 <!-- options-end -->
 
89
 *
 
90
 * @author Eibe Frank (eibe@cs.waikato.ac.nz)
 
91
 * @author Xin Xu (xx5@cs.waikato.ac.nz)
 
92
 * @version $Revision: 1.5 $ 
 
93
 */
 
94
public class TLD 
 
95
  extends RandomizableClassifier 
 
96
  implements OptionHandler, MultiInstanceCapabilitiesHandler,
 
97
             TechnicalInformationHandler {
 
98
 
 
99
  /** for serialization */
 
100
  static final long serialVersionUID = 6657315525171152210L;
 
101
  
 
102
  /** The mean for each attribute of each positive exemplar */
 
103
  protected double[][] m_MeanP = null;
 
104
 
 
105
  /** The variance for each attribute of each positive exemplar */
 
106
  protected double[][] m_VarianceP = null;
 
107
 
 
108
  /** The mean for each attribute of each negative exemplar */
 
109
  protected double[][] m_MeanN = null;
 
110
 
 
111
  /** The variance for each attribute of each negative exemplar */
 
112
  protected double[][] m_VarianceN = null;
 
113
 
 
114
  /** The effective sum of weights of each positive exemplar in each dimension*/
 
115
  protected double[][] m_SumP = null;
 
116
 
 
117
  /** The effective sum of weights of each negative exemplar in each dimension*/
 
118
  protected double[][] m_SumN = null;
 
119
 
 
120
  /** The parameters to be estimated for each positive exemplar*/
 
121
  protected double[] m_ParamsP = null;
 
122
 
 
123
  /** The parameters to be estimated for each negative exemplar*/
 
124
  protected double[] m_ParamsN = null;
 
125
 
 
126
  /** The dimension of each exemplar, i.e. (numAttributes-2) */
 
127
  protected int m_Dimension = 0;
 
128
 
 
129
  /** The class label of each exemplar */
 
130
  protected double[] m_Class = null;
 
131
 
 
132
  /** The number of class labels in the data */
 
133
  protected int m_NumClasses = 2;
 
134
 
 
135
  /** The very small number representing zero */
 
136
  static public double ZERO = 1.0e-6;   
 
137
 
 
138
  /** The number of runs to perform */
 
139
  protected int m_Run = 1;
 
140
 
 
141
  protected double m_Cutoff;
 
142
 
 
143
  protected boolean m_UseEmpiricalCutOff = false;   
 
144
 
 
145
  /**
 
146
   * Returns a string describing this filter
 
147
   *
 
148
   * @return a description of the filter suitable for
 
149
   * displaying in the explorer/experimenter gui
 
150
   */
 
151
  public String globalInfo() {
 
152
    return 
 
153
        "Two-Level Distribution approach, changes the starting value of "
 
154
      + "the searching algorithm, supplement the cut-off modification and "
 
155
      + "check missing values.\n\n"
 
156
      + "For more information see:\n\n"
 
157
      + getTechnicalInformation().toString();
 
158
  }
 
159
  
 
160
  /**
 
161
   * Returns an instance of a TechnicalInformation object, containing 
 
162
   * detailed information about the technical background of this class,
 
163
   * e.g., paper reference or book this class is based on.
 
164
   * 
 
165
   * @return the technical information about this class
 
166
   */
 
167
  public TechnicalInformation getTechnicalInformation() {
 
168
    TechnicalInformation        result;
 
169
    
 
170
    result = new TechnicalInformation(Type.MASTERSTHESIS);
 
171
    result.setValue(Field.AUTHOR, "Xin Xu");
 
172
    result.setValue(Field.YEAR, "2003");
 
173
    result.setValue(Field.TITLE, "Statistical learning in multiple instance problem");
 
174
    result.setValue(Field.SCHOOL, "University of Waikato");
 
175
    result.setValue(Field.ADDRESS, "Hamilton, NZ");
 
176
    result.setValue(Field.NOTE, "0657.594");
 
177
    
 
178
    return result;
 
179
  }
 
180
 
 
181
  /**
 
182
   * Returns default capabilities of the classifier.
 
183
   *
 
184
   * @return      the capabilities of this classifier
 
185
   */
 
186
  public Capabilities getCapabilities() {
 
187
    Capabilities result = super.getCapabilities();
 
188
 
 
189
    // attributes
 
190
    result.enable(Capability.NOMINAL_ATTRIBUTES);
 
191
    result.enable(Capability.RELATIONAL_ATTRIBUTES);
 
192
    result.enable(Capability.MISSING_VALUES);
 
193
 
 
194
    // class
 
195
    result.enable(Capability.BINARY_CLASS);
 
196
    result.enable(Capability.MISSING_CLASS_VALUES);
 
197
    
 
198
    // other
 
199
    result.enable(Capability.ONLY_MULTIINSTANCE);
 
200
    
 
201
    return result;
 
202
  }
 
203
 
 
204
  /**
 
205
   * Returns the capabilities of this multi-instance classifier for the
 
206
   * relational data.
 
207
   *
 
208
   * @return            the capabilities of this object
 
209
   * @see               Capabilities
 
210
   */
 
211
  public Capabilities getMultiInstanceCapabilities() {
 
212
    Capabilities result = super.getCapabilities();
 
213
    
 
214
    // attributes
 
215
    result.enable(Capability.NUMERIC_ATTRIBUTES);
 
216
    result.enable(Capability.MISSING_VALUES);
 
217
 
 
218
    // class
 
219
    result.disableAllClasses();
 
220
    result.enable(Capability.NO_CLASS);
 
221
    
 
222
    return result;
 
223
  }
 
224
 
 
225
  /**
 
226
   *
 
227
   * @param exs the training exemplars
 
228
   * @throws Exception if the model cannot be built properly
 
229
   */    
 
230
  public void buildClassifier(Instances exs)throws Exception{
 
231
    // can classifier handle the data?
 
232
    getCapabilities().testWithFail(exs);
 
233
 
 
234
    // remove instances with missing class
 
235
    exs = new Instances(exs);
 
236
    exs.deleteWithMissingClass();
 
237
    
 
238
    int numegs = exs.numInstances();
 
239
    m_Dimension = exs.attribute(1).relation(). numAttributes();
 
240
    Instances pos = new Instances(exs, 0), neg = new Instances(exs, 0);
 
241
 
 
242
    for(int u=0; u<numegs; u++){
 
243
      Instance example = exs.instance(u);
 
244
      if(example.classValue() == 1)
 
245
        pos.add(example);
 
246
      else
 
247
        neg.add(example);
 
248
    }
 
249
 
 
250
    int pnum = pos.numInstances(), nnum = neg.numInstances();   
 
251
 
 
252
    m_MeanP = new double[pnum][m_Dimension];
 
253
    m_VarianceP = new double[pnum][m_Dimension];
 
254
    m_SumP = new double[pnum][m_Dimension];
 
255
    m_MeanN = new double[nnum][m_Dimension];
 
256
    m_VarianceN = new double[nnum][m_Dimension];
 
257
    m_SumN = new double[nnum][m_Dimension];
 
258
    m_ParamsP = new double[4*m_Dimension];
 
259
    m_ParamsN = new double[4*m_Dimension];
 
260
 
 
261
    // Estimation of the parameters: as the start value for search
 
262
    double[] pSumVal=new double[m_Dimension], // for m 
 
263
      nSumVal=new double[m_Dimension]; 
 
264
    double[] maxVarsP=new double[m_Dimension], // for a
 
265
      maxVarsN=new double[m_Dimension]; 
 
266
    // Mean of sample variances: for b, b=a/E(\sigma^2)+2
 
267
    double[] varMeanP = new double[m_Dimension],
 
268
      varMeanN = new double[m_Dimension]; 
 
269
    // Variances of sample means: for w, w=E[var(\mu)]/E[\sigma^2]
 
270
    double[] meanVarP = new double[m_Dimension],
 
271
      meanVarN = new double[m_Dimension];
 
272
    // number of exemplars without all values missing
 
273
    double[] numExsP = new double[m_Dimension],
 
274
      numExsN = new double[m_Dimension];
 
275
 
 
276
    // Extract metadata fro both positive and negative bags
 
277
    for(int v=0; v < pnum; v++){
 
278
      /*Exemplar px = pos.exemplar(v);
 
279
        m_MeanP[v] = px.meanOrMode();
 
280
        m_VarianceP[v] = px.variance();
 
281
        Instances pxi =  px.getInstances();
 
282
        */
 
283
 
 
284
      Instances pxi =  pos.instance(v).relationalValue(1);
 
285
      for (int k=0; k<pxi.numAttributes(); k++) { 
 
286
        m_MeanP[v][k] = pxi.meanOrMode(k);
 
287
        m_VarianceP[v][k] = pxi.variance(k);
 
288
      }
 
289
 
 
290
      for (int w=0,t=0; w < m_Dimension; w++,t++){              
 
291
        //if((t==m_ClassIndex) || (t==m_IdIndex))
 
292
        //  t++;                
 
293
 
 
294
        if(!Double.isNaN(m_MeanP[v][w])){
 
295
          for(int u=0;u<pxi.numInstances();u++){
 
296
            Instance ins = pxi.instance(u);                     
 
297
            if(!ins.isMissing(t))
 
298
              m_SumP[v][w] += ins.weight();                        
 
299
          }   
 
300
          numExsP[w]++;  
 
301
          pSumVal[w] += m_MeanP[v][w];
 
302
          meanVarP[w] += m_MeanP[v][w]*m_MeanP[v][w];    
 
303
          if(maxVarsP[w] < m_VarianceP[v][w])
 
304
            maxVarsP[w] = m_VarianceP[v][w];
 
305
          varMeanP[w] += m_VarianceP[v][w];
 
306
          m_VarianceP[v][w] *= (m_SumP[v][w]-1.0);
 
307
          if(m_VarianceP[v][w] < 0.0)
 
308
            m_VarianceP[v][w] = 0.0;
 
309
        }
 
310
      }
 
311
    }
 
312
 
 
313
    for(int v=0; v < nnum; v++){
 
314
      /*Exemplar nx = neg.exemplar(v);
 
315
        m_MeanN[v] = nx.meanOrMode();
 
316
        m_VarianceN[v] = nx.variance();
 
317
        Instances nxi =  nx.getInstances();
 
318
        */
 
319
      Instances nxi =  neg.instance(v).relationalValue(1);
 
320
      for (int k=0; k<nxi.numAttributes(); k++) {
 
321
        m_MeanN[v][k] = nxi.meanOrMode(k);
 
322
        m_VarianceN[v][k] = nxi.variance(k);
 
323
      }
 
324
 
 
325
      for (int w=0,t=0; w < m_Dimension; w++,t++){              
 
326
        //if((t==m_ClassIndex) || (t==m_IdIndex))
 
327
        //  t++;                
 
328
 
 
329
        if(!Double.isNaN(m_MeanN[v][w])){
 
330
          for(int u=0;u<nxi.numInstances();u++)
 
331
            if(!nxi.instance(u).isMissing(t))
 
332
              m_SumN[v][w] += nxi.instance(u).weight();
 
333
          numExsN[w]++;         
 
334
          nSumVal[w] += m_MeanN[v][w];
 
335
          meanVarN[w] += m_MeanN[v][w]*m_MeanN[v][w]; 
 
336
          if(maxVarsN[w] < m_VarianceN[v][w])
 
337
            maxVarsN[w] = m_VarianceN[v][w];
 
338
          varMeanN[w] += m_VarianceN[v][w];
 
339
          m_VarianceN[v][w] *= (m_SumN[v][w]-1.0);
 
340
          if(m_VarianceN[v][w] < 0.0)
 
341
            m_VarianceN[v][w] = 0.0;
 
342
        }
 
343
      }
 
344
    }
 
345
 
 
346
    for(int w=0; w<m_Dimension; w++){
 
347
      pSumVal[w] /= numExsP[w];
 
348
      nSumVal[w] /= numExsN[w];
 
349
      if(numExsP[w]>1)
 
350
        meanVarP[w] = meanVarP[w]/(numExsP[w]-1.0) 
 
351
          - pSumVal[w]*numExsP[w]/(numExsP[w]-1.0);
 
352
      if(numExsN[w]>1)
 
353
        meanVarN[w] = meanVarN[w]/(numExsN[w]-1.0) 
 
354
          - nSumVal[w]*numExsN[w]/(numExsN[w]-1.0);
 
355
      varMeanP[w] /= numExsP[w];
 
356
      varMeanN[w] /= numExsN[w];
 
357
    }
 
358
 
 
359
    //Bounds and parameter values for each run
 
360
    double[][] bounds = new double[2][4];
 
361
    double[] pThisParam = new double[4], 
 
362
      nThisParam = new double[4];
 
363
 
 
364
    // Initial values for parameters
 
365
    double a, b, w, m;
 
366
 
 
367
    // Optimize for one dimension
 
368
    for (int x=0; x < m_Dimension; x++){
 
369
      if (getDebug())
 
370
        System.err.println("\n\n!!!!!!!!!!!!!!!!!!!!!!???Dimension #"+x);
 
371
 
 
372
      // Positive examplars: first run
 
373
      a = (maxVarsP[x]>ZERO) ? maxVarsP[x]:1.0; 
 
374
      if (varMeanP[x]<=ZERO)   varMeanP[x] = ZERO;  // modified by LinDong (09/2005)
 
375
      b = a/varMeanP[x]+2.0; // a/(b-2) = E(\sigma^2)
 
376
      w = meanVarP[x]/varMeanP[x]; // E[var(\mu)] = w*E[\sigma^2]           
 
377
      if(w<=ZERO)  w=1.0;
 
378
 
 
379
      m = pSumVal[x];     
 
380
      pThisParam[0] = a;    // a
 
381
      pThisParam[1] = b;  // b
 
382
      pThisParam[2] = w;  // w
 
383
      pThisParam[3] = m;  // m
 
384
 
 
385
      // Negative examplars: first run
 
386
      a = (maxVarsN[x]>ZERO) ? maxVarsN[x]:1.0; 
 
387
      if (varMeanN[x]<=ZERO)   varMeanN[x] = ZERO; // modified by LinDong (09/2005)
 
388
      b = a/varMeanN[x]+2.0; // a/(b-2) = E(\sigma^2)
 
389
      w = meanVarN[x]/varMeanN[x]; // E[var(\mu)] = w*E[\sigma^2]           
 
390
      if(w<=ZERO) w=1.0;
 
391
 
 
392
      m = nSumVal[x];     
 
393
      nThisParam[0] = a;    // a
 
394
      nThisParam[1] = b;  // b
 
395
      nThisParam[2] = w;  // w
 
396
      nThisParam[3] = m;  // m
 
397
 
 
398
      // Bound constraints
 
399
      bounds[0][0] = ZERO; // a > 0
 
400
      bounds[0][1] = 2.0+ZERO;  // b > 2 
 
401
      bounds[0][2] = ZERO; // w > 0
 
402
      bounds[0][3] = Double.NaN;
 
403
 
 
404
      for(int t=0; t<4; t++){
 
405
        bounds[1][t] = Double.NaN;
 
406
        m_ParamsP[4*x+t] = pThisParam[t];       
 
407
        m_ParamsN[4*x+t] = nThisParam[t];
 
408
      }
 
409
      double pminVal=Double.MAX_VALUE, nminVal=Double.MAX_VALUE;
 
410
      Random whichEx = new Random(m_Seed); 
 
411
      TLD_Optm pOp=null, nOp=null;      
 
412
      boolean isRunValid = true;
 
413
      double[] sumP=new double[pnum], meanP=new double[pnum],
 
414
        varP=new double[pnum];
 
415
      double[] sumN=new double[nnum], meanN=new double[nnum],
 
416
        varN=new double[nnum];
 
417
 
 
418
      // One dimension
 
419
      for(int p=0; p<pnum; p++){
 
420
        sumP[p] = m_SumP[p][x];
 
421
        meanP[p] = m_MeanP[p][x];
 
422
        varP[p] = m_VarianceP[p][x];
 
423
      }
 
424
      for(int q=0; q<nnum; q++){
 
425
        sumN[q] = m_SumN[q][x];
 
426
        meanN[q] = m_MeanN[q][x];
 
427
        varN[q] = m_VarianceN[q][x];
 
428
      }
 
429
 
 
430
      for(int y=0; y<m_Run;){
 
431
        if (getDebug())
 
432
          System.err.println("\n\n!!!!!!!!!!!!!!!!!!!!!!???Run #"+y);
 
433
        double thisMin;
 
434
 
 
435
        if (getDebug())
 
436
          System.err.println("\nPositive exemplars");
 
437
        pOp = new TLD_Optm();
 
438
        pOp.setNum(sumP);
 
439
        pOp.setSSquare(varP);
 
440
        pOp.setXBar(meanP);
 
441
 
 
442
        pThisParam = pOp.findArgmin(pThisParam, bounds);
 
443
        while(pThisParam==null){
 
444
          pThisParam = pOp.getVarbValues();                 
 
445
          if (getDebug())
 
446
            System.err.println("!!! 200 iterations finished, not enough!");
 
447
          pThisParam = pOp.findArgmin(pThisParam, bounds);
 
448
        }       
 
449
 
 
450
        thisMin = pOp.getMinFunction();
 
451
        if(!Double.isNaN(thisMin) && (thisMin<pminVal)){
 
452
          pminVal = thisMin;
 
453
          for(int z=0; z<4; z++)
 
454
            m_ParamsP[4*x+z] = pThisParam[z];
 
455
        }
 
456
 
 
457
        if(Double.isNaN(thisMin)){
 
458
          pThisParam = new double[4];
 
459
          isRunValid =false;
 
460
        }
 
461
 
 
462
        if (getDebug())
 
463
          System.err.println("\nNegative exemplars");
 
464
        nOp = new TLD_Optm();
 
465
        nOp.setNum(sumN);
 
466
        nOp.setSSquare(varN);
 
467
        nOp.setXBar(meanN);
 
468
 
 
469
        nThisParam = nOp.findArgmin(nThisParam, bounds);
 
470
        while(nThisParam==null){
 
471
          nThisParam = nOp.getVarbValues();
 
472
          if (getDebug())
 
473
            System.err.println("!!! 200 iterations finished, not enough!");
 
474
          nThisParam = nOp.findArgmin(nThisParam, bounds);
 
475
        }       
 
476
        thisMin = nOp.getMinFunction();
 
477
        if(!Double.isNaN(thisMin) && (thisMin<nminVal)){
 
478
          nminVal = thisMin;
 
479
          for(int z=0; z<4; z++)
 
480
            m_ParamsN[4*x+z] = nThisParam[z];     
 
481
        }
 
482
 
 
483
        if(Double.isNaN(thisMin)){
 
484
          nThisParam = new double[4];
 
485
          isRunValid =false;
 
486
        }
 
487
 
 
488
        if(!isRunValid){ y--; isRunValid=true; }                
 
489
 
 
490
        if(++y<m_Run){
 
491
          // Change the initial parameters and restart              
 
492
          int pone = whichEx.nextInt(pnum), // Randomly pick one pos. exmpl.
 
493
              none = whichEx.nextInt(nnum);
 
494
 
 
495
          // Positive exemplars: next run 
 
496
          while((m_SumP[pone][x]<=1.0)||Double.isNaN(m_MeanP[pone][x]))
 
497
            pone = whichEx.nextInt(pnum);
 
498
 
 
499
          a = m_VarianceP[pone][x]/(m_SumP[pone][x]-1.0);               
 
500
          if(a<=ZERO) a=m_ParamsN[4*x]; // Change to negative params
 
501
          m = m_MeanP[pone][x];
 
502
          double sq = (m-m_ParamsP[4*x+3])*(m-m_ParamsP[4*x+3]);
 
503
 
 
504
          b = a*m_ParamsP[4*x+2]/sq+2.0; // b=a/Var+2, assuming Var=Sq/w'
 
505
          if((b<=ZERO) || Double.isNaN(b) || Double.isInfinite(b))
 
506
            b=m_ParamsN[4*x+1];
 
507
 
 
508
          w = sq*(m_ParamsP[4*x+1]-2.0)/m_ParamsP[4*x];//w=Sq/Var, assuming Var=a'/(b'-2)
 
509
          if((w<=ZERO) || Double.isNaN(w) || Double.isInfinite(w))
 
510
            w=m_ParamsN[4*x+2];
 
511
 
 
512
          pThisParam[0] = a;    // a
 
513
          pThisParam[1] = b;  // b
 
514
          pThisParam[2] = w;  // w
 
515
          pThisParam[3] = m;  // m          
 
516
 
 
517
          // Negative exemplars: next run 
 
518
          while((m_SumN[none][x]<=1.0)||Double.isNaN(m_MeanN[none][x]))
 
519
            none = whichEx.nextInt(nnum);           
 
520
 
 
521
          a = m_VarianceN[none][x]/(m_SumN[none][x]-1.0);       
 
522
          if(a<=ZERO) a=m_ParamsP[4*x];       
 
523
          m = m_MeanN[none][x];
 
524
          sq = (m-m_ParamsN[4*x+3])*(m-m_ParamsN[4*x+3]);
 
525
 
 
526
          b = a*m_ParamsN[4*x+2]/sq+2.0; // b=a/Var+2, assuming Var=Sq/w'
 
527
          if((b<=ZERO) || Double.isNaN(b) || Double.isInfinite(b))
 
528
            b=m_ParamsP[4*x+1];
 
529
 
 
530
          w = sq*(m_ParamsN[4*x+1]-2.0)/m_ParamsN[4*x];//w=Sq/Var, assuming Var=a'/(b'-2)
 
531
          if((w<=ZERO) || Double.isNaN(w) || Double.isInfinite(w))
 
532
            w=m_ParamsP[4*x+2];
 
533
 
 
534
          nThisParam[0] = a;    // a
 
535
          nThisParam[1] = b;  // b
 
536
          nThisParam[2] = w;  // w
 
537
          nThisParam[3] = m;  // m                      
 
538
        }
 
539
      }                     
 
540
    }
 
541
 
 
542
    for (int x=0, y=0; x<m_Dimension; x++, y++){
 
543
      //if((x==exs.classIndex()) || (x==exs.idIndex()))
 
544
      //y++;
 
545
      a=m_ParamsP[4*x]; b=m_ParamsP[4*x+1]; 
 
546
      w=m_ParamsP[4*x+2]; m=m_ParamsP[4*x+3];
 
547
      if (getDebug())
 
548
        System.err.println("\n\n???Positive: ( "+exs.attribute(1).relation().attribute(y)+
 
549
          "): a="+a+", b="+b+", w="+w+", m="+m);
 
550
 
 
551
      a=m_ParamsN[4*x]; b=m_ParamsN[4*x+1]; 
 
552
      w=m_ParamsN[4*x+2]; m=m_ParamsN[4*x+3];
 
553
      if (getDebug())
 
554
        System.err.println("???Negative: ("+exs.attribute(1).relation().attribute(y)+
 
555
          "): a="+a+", b="+b+", w="+w+", m="+m);
 
556
    }
 
557
 
 
558
    if(m_UseEmpiricalCutOff){   
 
559
      // Find the empirical cut-off
 
560
      double[] pLogOdds=new double[pnum], nLogOdds=new double[nnum];  
 
561
      for(int p=0; p<pnum; p++)
 
562
        pLogOdds[p] = 
 
563
          likelihoodRatio(m_SumP[p], m_MeanP[p], m_VarianceP[p]);
 
564
 
 
565
      for(int q=0; q<nnum; q++)
 
566
        nLogOdds[q] = 
 
567
          likelihoodRatio(m_SumN[q], m_MeanN[q], m_VarianceN[q]);
 
568
 
 
569
      // Update m_Cutoff
 
570
      findCutOff(pLogOdds, nLogOdds);
 
571
    }
 
572
    else
 
573
      m_Cutoff = -Math.log((double)pnum/(double)nnum);
 
574
 
 
575
    if (getDebug())
 
576
      System.err.println("???Cut-off="+m_Cutoff);
 
577
  }        
 
578
 
 
579
  /**
 
580
   *
 
581
   * @param ex the given test exemplar
 
582
   * @return the classification 
 
583
   * @throws Exception if the exemplar could not be classified
 
584
   * successfully
 
585
   */
 
586
  public double classifyInstance(Instance ex)throws Exception{
 
587
    //Exemplar ex = new Exemplar(e);
 
588
    Instances exi = ex.relationalValue(1);
 
589
    double[] n = new double[m_Dimension];
 
590
    double [] xBar = new double[m_Dimension];
 
591
    double [] sSq = new double[m_Dimension];
 
592
    for (int i=0; i<exi.numAttributes() ; i++){
 
593
      xBar[i] = exi.meanOrMode(i);
 
594
      sSq[i] = exi.variance(i);
 
595
    }
 
596
 
 
597
    for (int w=0, t=0; w < m_Dimension; w++, t++){
 
598
      //if((t==m_ClassIndex) || (t==m_IdIndex))
 
599
      //t++;    
 
600
      for(int u=0;u<exi.numInstances();u++)
 
601
        if(!exi.instance(u).isMissing(t))
 
602
          n[w] += exi.instance(u).weight();
 
603
 
 
604
      sSq[w] = sSq[w]*(n[w]-1.0);
 
605
      if(sSq[w] <= 0.0)
 
606
        sSq[w] = 0.0;
 
607
    }
 
608
 
 
609
    double logOdds = likelihoodRatio(n, xBar, sSq);
 
610
    return (logOdds > m_Cutoff) ? 1 : 0 ;
 
611
  }
 
612
 
 
613
  private double likelihoodRatio(double[] n, double[] xBar, double[] sSq){      
 
614
    double LLP = 0.0, LLN = 0.0;
 
615
 
 
616
    for (int x=0; x<m_Dimension; x++){
 
617
      if(Double.isNaN(xBar[x])) continue; // All missing values
 
618
 
 
619
      int halfN = ((int)n[x])/2;        
 
620
      //Log-likelihood for positive 
 
621
      double a=m_ParamsP[4*x], b=m_ParamsP[4*x+1], 
 
622
             w=m_ParamsP[4*x+2], m=m_ParamsP[4*x+3];
 
623
      LLP += 0.5*b*Math.log(a) + 0.5*(b+n[x]-1.0)*Math.log(1.0+n[x]*w)
 
624
        - 0.5*(b+n[x])*Math.log((1.0+n[x]*w)*(a+sSq[x])+
 
625
            n[x]*(xBar[x]-m)*(xBar[x]-m))
 
626
        - 0.5*n[x]*Math.log(Math.PI);
 
627
      for(int y=1; y<=halfN; y++)
 
628
        LLP += Math.log(b/2.0+n[x]/2.0-(double)y);
 
629
 
 
630
      if(n[x]/2.0 > halfN) // n is odd
 
631
        LLP += TLD_Optm.diffLnGamma(b/2.0);
 
632
 
 
633
      //Log-likelihood for negative 
 
634
      a=m_ParamsN[4*x];
 
635
      b=m_ParamsN[4*x+1]; 
 
636
      w=m_ParamsN[4*x+2];
 
637
      m=m_ParamsN[4*x+3];
 
638
      LLN += 0.5*b*Math.log(a) + 0.5*(b+n[x]-1.0)*Math.log(1.0+n[x]*w)
 
639
        - 0.5*(b+n[x])*Math.log((1.0+n[x]*w)*(a+sSq[x])+
 
640
            n[x]*(xBar[x]-m)*(xBar[x]-m))
 
641
        - 0.5*n[x]*Math.log(Math.PI);
 
642
      for(int y=1; y<=halfN; y++)
 
643
        LLN += Math.log(b/2.0+n[x]/2.0-(double)y);      
 
644
 
 
645
      if(n[x]/2.0 > halfN) // n is odd
 
646
        LLN += TLD_Optm.diffLnGamma(b/2.0);   
 
647
    }
 
648
 
 
649
    return LLP - LLN;
 
650
  }
 
651
 
 
652
  private void findCutOff(double[] pos, double[] neg){
 
653
    int[] pOrder = Utils.sort(pos),
 
654
      nOrder = Utils.sort(neg);
 
655
    /*
 
656
       System.err.println("\n\n???Positive: ");
 
657
       for(int t=0; t<pOrder.length; t++)
 
658
       System.err.print(t+":"+Utils.doubleToString(pos[pOrder[t]],0,2)+" ");
 
659
       System.err.println("\n\n???Negative: ");
 
660
       for(int t=0; t<nOrder.length; t++)
 
661
       System.err.print(t+":"+Utils.doubleToString(neg[nOrder[t]],0,2)+" ");
 
662
       */
 
663
    int pNum = pos.length, nNum = neg.length, count, p=0, n=0;  
 
664
    double fstAccu=0.0, sndAccu=(double)pNum, split; 
 
665
    double maxAccu = 0, minDistTo0 = Double.MAX_VALUE;
 
666
 
 
667
    // Skip continuous negatives        
 
668
    for(;(n<nNum)&&(pos[pOrder[0]]>=neg[nOrder[n]]); n++, fstAccu++);
 
669
 
 
670
    if(n>=nNum){ // totally seperate
 
671
      m_Cutoff = (neg[nOrder[nNum-1]]+pos[pOrder[0]])/2.0;      
 
672
      //m_Cutoff = neg[nOrder[nNum-1]];
 
673
      return;  
 
674
    }   
 
675
 
 
676
    count=n;
 
677
    while((p<pNum)&&(n<nNum)){
 
678
      // Compare the next in the two lists
 
679
      if(pos[pOrder[p]]>=neg[nOrder[n]]){ // Neg has less log-odds
 
680
        fstAccu += 1.0;    
 
681
        split=neg[nOrder[n]];
 
682
        n++;     
 
683
      }
 
684
      else{
 
685
        sndAccu -= 1.0;
 
686
        split=pos[pOrder[p]];
 
687
        p++;
 
688
      }           
 
689
      count++;
 
690
      if((fstAccu+sndAccu > maxAccu) 
 
691
          || ((fstAccu+sndAccu == maxAccu) && (Math.abs(split)<minDistTo0))){
 
692
        maxAccu = fstAccu+sndAccu;
 
693
        m_Cutoff = split;
 
694
        minDistTo0 = Math.abs(split);
 
695
      }     
 
696
    }           
 
697
  }
 
698
 
 
699
  /**
 
700
   * Returns an enumeration describing the available options
 
701
   *
 
702
   * @return an enumeration of all the available options
 
703
   */
 
704
  public Enumeration listOptions() {
 
705
    Vector result = new Vector();
 
706
    
 
707
    result.addElement(new Option(
 
708
          "\tSet whether or not use empirical\n"
 
709
          + "\tlog-odds cut-off instead of 0",
 
710
          "C", 0, "-C"));
 
711
    
 
712
    result.addElement(new Option(
 
713
          "\tSet the number of multiple runs \n"
 
714
          + "\tneeded for searching the MLE.",
 
715
          "R", 1, "-R <numOfRuns>"));
 
716
    
 
717
    Enumeration enu = super.listOptions();
 
718
    while (enu.hasMoreElements()) {
 
719
      result.addElement(enu.nextElement());
 
720
    }
 
721
 
 
722
    return result.elements();
 
723
  }
 
724
 
 
725
  /**
 
726
   * Parses a given list of options. <p/>
 
727
   * 
 
728
   <!-- options-start -->
 
729
   * Valid options are: <p/>
 
730
   * 
 
731
   * <pre> -C
 
732
   *  Set whether or not use empirical
 
733
   *  log-odds cut-off instead of 0</pre>
 
734
   * 
 
735
   * <pre> -R &lt;numOfRuns&gt;
 
736
   *  Set the number of multiple runs 
 
737
   *  needed for searching the MLE.</pre>
 
738
   * 
 
739
   * <pre> -S &lt;num&gt;
 
740
   *  Random number seed.
 
741
   *  (default 1)</pre>
 
742
   * 
 
743
   * <pre> -D
 
744
   *  If set, classifier is run in debug mode and
 
745
   *  may output additional info to the console</pre>
 
746
   * 
 
747
   <!-- options-end -->
 
748
   *
 
749
   * @param options the list of options as an array of strings
 
750
   * @throws Exception if an option is not supported
 
751
   */
 
752
  public void setOptions(String[] options) throws Exception{
 
753
    setDebug(Utils.getFlag('D', options));
 
754
 
 
755
    setUsingCutOff(Utils.getFlag('C', options));
 
756
 
 
757
    String runString = Utils.getOption('R', options);
 
758
    if (runString.length() != 0) 
 
759
      setNumRuns(Integer.parseInt(runString));
 
760
    else 
 
761
      setNumRuns(1);
 
762
 
 
763
    super.setOptions(options);
 
764
  }
 
765
 
 
766
  /**
 
767
   * Gets the current settings of the Classifier.
 
768
   *
 
769
   * @return an array of strings suitable for passing to setOptions
 
770
   */
 
771
  public String[] getOptions() {
 
772
    Vector        result;
 
773
    String[]      options;
 
774
    int           i;
 
775
    
 
776
    result  = new Vector();
 
777
    options = super.getOptions();
 
778
    for (i = 0; i < options.length; i++)
 
779
      result.add(options[i]);
 
780
 
 
781
    if (getDebug())
 
782
      result.add("-D");
 
783
    
 
784
    if (getUsingCutOff())
 
785
      result.add("-C");
 
786
 
 
787
    result.add("-R");
 
788
    result.add("" + getNumRuns());
 
789
 
 
790
    return (String[]) result.toArray(new String[result.size()]);
 
791
  }
 
792
 
 
793
  /**
 
794
   * Returns the tip text for this property
 
795
   *
 
796
   * @return tip text for this property suitable for
 
797
   * displaying in the explorer/experimenter gui
 
798
   */
 
799
  public String numRunsTipText() {
 
800
    return "The number of runs to perform.";
 
801
  }
 
802
 
 
803
  /**
 
804
   * Sets the number of runs to perform.
 
805
   *
 
806
   * @param numRuns   the number of runs to perform
 
807
   */
 
808
  public void setNumRuns(int numRuns) {
 
809
    m_Run = numRuns;
 
810
  }
 
811
 
 
812
  /**
 
813
   * Returns the number of runs to perform.
 
814
   *
 
815
   * @return          the number of runs to perform
 
816
   */
 
817
  public int getNumRuns() {
 
818
    return m_Run;
 
819
  }
 
820
 
 
821
  /**
 
822
   * Returns the tip text for this property
 
823
   *
 
824
   * @return tip text for this property suitable for
 
825
   * displaying in the explorer/experimenter gui
 
826
   */
 
827
  public String usingCutOffTipText() {
 
828
    return "Whether to use an empirical cutoff.";
 
829
  }
 
830
 
 
831
  /**
 
832
   * Sets whether to use an empirical cutoff.
 
833
   *
 
834
   * @param cutOff      whether to use an empirical cutoff
 
835
   */
 
836
  public void setUsingCutOff (boolean cutOff) {
 
837
    m_UseEmpiricalCutOff = cutOff;
 
838
  }
 
839
 
 
840
  /**
 
841
   * Returns whether an empirical cutoff is used
 
842
   *
 
843
   * @return            true if an empirical cutoff is used
 
844
   */
 
845
  public boolean getUsingCutOff() {
 
846
    return m_UseEmpiricalCutOff;
 
847
  }
 
848
 
 
849
  /**
 
850
   * Main method for testing.
 
851
   *
 
852
   * @param args the options for the classifier
 
853
   */
 
854
  public static void main(String[] args) {      
 
855
    runClassifier(new TLD(), args);
 
856
  }
 
857
}
 
858
 
 
859
class TLD_Optm extends Optimization{
 
860
 
 
861
  private double[] num;
 
862
  private double[] sSq;
 
863
  private double[] xBar;
 
864
 
 
865
  public void setNum(double[] n) {num = n;}
 
866
  public void setSSquare(double[] s){sSq = s;}
 
867
  public void setXBar(double[] x){xBar = x;}
 
868
 
 
869
  /**
 
870
   * Compute Ln[Gamma(b+0.5)] - Ln[Gamma(b)]
 
871
   *
 
872
   * @param b the value in the above formula
 
873
   * @return the result
 
874
   */    
 
875
  public static double diffLnGamma(double b){
 
876
    double[] coef= {76.18009172947146, -86.50532032941677,
 
877
      24.01409824083091, -1.231739572450155, 
 
878
      0.1208650973866179e-2, -0.5395239384953e-5};
 
879
    double rt = -0.5;
 
880
    rt += (b+1.0)*Math.log(b+6.0) - (b+0.5)*Math.log(b+5.5);
 
881
    double series1=1.000000000190015, series2=1.000000000190015;
 
882
    for(int i=0; i<6; i++){
 
883
      series1 += coef[i]/(b+1.5+(double)i);
 
884
      series2 += coef[i]/(b+1.0+(double)i);
 
885
    }
 
886
 
 
887
    rt += Math.log(series1*b)-Math.log(series2*(b+0.5));
 
888
    return rt;
 
889
  }
 
890
 
 
891
  /**
 
892
   * Compute dLn[Gamma(x+0.5)]/dx - dLn[Gamma(x)]/dx
 
893
   *
 
894
   * @param x the value in the above formula
 
895
   * @return the result
 
896
   */    
 
897
  protected double diffFstDervLnGamma(double x){
 
898
    double rt=0, series=1.0;// Just make it >0
 
899
    for(int i=0;series>=m_Zero*1e-3;i++){
 
900
      series = 0.5/((x+(double)i)*(x+(double)i+0.5));
 
901
      rt += series;
 
902
    }
 
903
    return rt;
 
904
  }
 
905
 
 
906
  /**
 
907
   * Compute {Ln[Gamma(x+0.5)]}'' - {Ln[Gamma(x)]}''
 
908
   *
 
909
   * @param x the value in the above formula
 
910
   * @return the result
 
911
   */    
 
912
  protected double diffSndDervLnGamma(double x){
 
913
    double rt=0, series=1.0;// Just make it >0
 
914
    for(int i=0;series>=m_Zero*1e-3;i++){
 
915
      series = (x+(double)i+0.25)/
 
916
        ((x+(double)i)*(x+(double)i)*(x+(double)i+0.5)*(x+(double)i+0.5));
 
917
      rt -= series;
 
918
    }
 
919
    return rt;
 
920
  }
 
921
 
 
922
  /**
 
923
   * Implement this procedure to evaluate objective
 
924
   * function to be minimized
 
925
   */
 
926
  protected double objectiveFunction(double[] x){
 
927
    int numExs = num.length;
 
928
    double NLL = 0; // Negative Log-Likelihood
 
929
 
 
930
    double a=x[0], b=x[1], w=x[2], m=x[3];
 
931
    for(int j=0; j < numExs; j++){
 
932
 
 
933
      if(Double.isNaN(xBar[j])) continue; // All missing values
 
934
 
 
935
      NLL += 0.5*(b+num[j])*
 
936
        Math.log((1.0+num[j]*w)*(a+sSq[j]) + 
 
937
            num[j]*(xBar[j]-m)*(xBar[j]-m));        
 
938
 
 
939
      if(Double.isNaN(NLL) && m_Debug){
 
940
        System.err.println("???????????1: "+a+" "+b+" "+w+" "+m
 
941
            +"|x-: "+xBar[j] + 
 
942
            "|n: "+num[j] + "|S^2: "+sSq[j]);
 
943
        System.exit(1);
 
944
      }
 
945
 
 
946
      // Doesn't affect optimization
 
947
      //NLL += 0.5*num[j]*Math.log(Math.PI);            
 
948
 
 
949
      NLL -= 0.5*(b+num[j]-1.0)*Math.log(1.0+num[j]*w);
 
950
 
 
951
 
 
952
      if(Double.isNaN(NLL) && m_Debug){
 
953
        System.err.println("???????????2: "+a+" "+b+" "+w+" "+m
 
954
            +"|x-: "+xBar[j] + 
 
955
            "|n: "+num[j] + "|S^2: "+sSq[j]);
 
956
        System.exit(1);
 
957
      }
 
958
 
 
959
      int halfNum = ((int)num[j])/2;
 
960
      for(int z=1; z<=halfNum; z++)
 
961
        NLL -= Math.log(0.5*b+0.5*num[j]-(double)z);
 
962
 
 
963
      if(0.5*num[j] > halfNum) // num[j] is odd
 
964
        NLL -= diffLnGamma(0.5*b);
 
965
 
 
966
      if(Double.isNaN(NLL) && m_Debug){
 
967
        System.err.println("???????????3: "+a+" "+b+" "+w+" "+m
 
968
            +"|x-: "+xBar[j] + 
 
969
            "|n: "+num[j] + "|S^2: "+sSq[j]);
 
970
        System.exit(1);
 
971
      }                         
 
972
 
 
973
      NLL -= 0.5*Math.log(a)*b;
 
974
      if(Double.isNaN(NLL) && m_Debug){
 
975
        System.err.println("???????????4:"+a+" "+b+" "+w+" "+m);
 
976
        System.exit(1);
 
977
      }     
 
978
    }
 
979
    if(m_Debug)
 
980
      System.err.println("?????????????5: "+NLL);
 
981
    if(Double.isNaN(NLL))          
 
982
      System.exit(1);
 
983
 
 
984
    return NLL;
 
985
  }
 
986
 
 
987
  /**
 
988
   * Subclass should implement this procedure to evaluate gradient
 
989
   * of the objective function
 
990
   */
 
991
  protected double[] evaluateGradient(double[] x){
 
992
    double[] g = new double[x.length];
 
993
    int numExs = num.length;
 
994
 
 
995
    double a=x[0],b=x[1],w=x[2],m=x[3];
 
996
 
 
997
    double da=0.0, db=0.0, dw=0.0, dm=0.0; 
 
998
    for(int j=0; j < numExs; j++){
 
999
 
 
1000
      if(Double.isNaN(xBar[j])) continue; // All missing values
 
1001
 
 
1002
      double denorm = (1.0+num[j]*w)*(a+sSq[j]) + 
 
1003
        num[j]*(xBar[j]-m)*(xBar[j]-m);
 
1004
 
 
1005
      da += 0.5*(b+num[j])*(1.0+num[j]*w)/denorm-0.5*b/a;
 
1006
 
 
1007
      db += 0.5*Math.log(denorm) 
 
1008
        - 0.5*Math.log(1.0+num[j]*w)
 
1009
        - 0.5*Math.log(a);
 
1010
 
 
1011
      int halfNum = ((int)num[j])/2;
 
1012
      for(int z=1; z<=halfNum; z++)
 
1013
        db -= 1.0/(b+num[j]-2.0*(double)z);             
 
1014
      if(num[j]/2.0 > halfNum) // num[j] is odd
 
1015
        db -= 0.5*diffFstDervLnGamma(0.5*b);            
 
1016
 
 
1017
      dw += 0.5*(b+num[j])*(a+sSq[j])*num[j]/denorm -
 
1018
        0.5*(b+num[j]-1.0)*num[j]/(1.0+num[j]*w);
 
1019
 
 
1020
      dm += num[j]*(b+num[j])*(m-xBar[j])/denorm;
 
1021
    }
 
1022
 
 
1023
    g[0] = da;
 
1024
    g[1] = db;
 
1025
    g[2] = dw;
 
1026
    g[3] = dm;
 
1027
    return g;
 
1028
  }
 
1029
 
 
1030
  /**
 
1031
   * Subclass should implement this procedure to evaluate second-order
 
1032
   * gradient of the objective function
 
1033
   */
 
1034
  protected double[] evaluateHessian(double[] x, int index){
 
1035
    double[] h = new double[x.length];
 
1036
 
 
1037
    // # of exemplars, # of dimensions
 
1038
    // which dimension and which variable for 'index'
 
1039
    int numExs = num.length;
 
1040
    double a,b,w,m;
 
1041
    // Take the 2nd-order derivative
 
1042
    switch(index){
 
1043
      case 0:  // a        
 
1044
        a=x[0];b=x[1];w=x[2];m=x[3];
 
1045
 
 
1046
        for(int j=0; j < numExs; j++){
 
1047
          if(Double.isNaN(xBar[j])) continue; //All missing values
 
1048
          double denorm = (1.0+num[j]*w)*(a+sSq[j]) + 
 
1049
            num[j]*(xBar[j]-m)*(xBar[j]-m);
 
1050
 
 
1051
          h[0] += 0.5*b/(a*a) 
 
1052
            - 0.5*(b+num[j])*(1.0+num[j]*w)*(1.0+num[j]*w)
 
1053
            /(denorm*denorm);
 
1054
 
 
1055
          h[1] += 0.5*(1.0+num[j]*w)/denorm - 0.5/a;
 
1056
 
 
1057
          h[2] += 0.5*num[j]*num[j]*(b+num[j])*
 
1058
            (xBar[j]-m)*(xBar[j]-m)/(denorm*denorm);
 
1059
 
 
1060
          h[3] -= num[j]*(b+num[j])*(m-xBar[j])
 
1061
            *(1.0+num[j]*w)/(denorm*denorm);
 
1062
        }
 
1063
        break;
 
1064
 
 
1065
      case 1: // b      
 
1066
        a=x[0];b=x[1];w=x[2];m=x[3];
 
1067
 
 
1068
        for(int j=0; j < numExs; j++){
 
1069
          if(Double.isNaN(xBar[j])) continue; //All missing values
 
1070
          double denorm = (1.0+num[j]*w)*(a+sSq[j]) + 
 
1071
            num[j]*(xBar[j]-m)*(xBar[j]-m);
 
1072
 
 
1073
          h[0] += 0.5*(1.0+num[j]*w)/denorm - 0.5/a;
 
1074
 
 
1075
          int halfNum = ((int)num[j])/2;
 
1076
          for(int z=1; z<=halfNum; z++)
 
1077
            h[1] += 
 
1078
              1.0/((b+num[j]-2.0*(double)z)*(b+num[j]-2.0*(double)z));
 
1079
          if(num[j]/2.0 > halfNum) // num[j] is odd
 
1080
            h[1] -= 0.25*diffSndDervLnGamma(0.5*b); 
 
1081
 
 
1082
          h[2] += 0.5*(a+sSq[j])*num[j]/denorm -
 
1083
            0.5*num[j]/(1.0+num[j]*w);
 
1084
 
 
1085
          h[3] += num[j]*(m-xBar[j])/denorm;
 
1086
        }
 
1087
        break;
 
1088
 
 
1089
      case 2: // w   
 
1090
        a=x[0];b=x[1];w=x[2];m=x[3];
 
1091
 
 
1092
        for(int j=0; j < numExs; j++){
 
1093
          if(Double.isNaN(xBar[j])) continue; //All missing values
 
1094
          double denorm = (1.0+num[j]*w)*(a+sSq[j]) + 
 
1095
            num[j]*(xBar[j]-m)*(xBar[j]-m);
 
1096
 
 
1097
          h[0] += 0.5*num[j]*num[j]*(b+num[j])*
 
1098
            (xBar[j]-m)*(xBar[j]-m)/(denorm*denorm);
 
1099
 
 
1100
          h[1] += 0.5*(a+sSq[j])*num[j]/denorm -
 
1101
            0.5*num[j]/(1.0+num[j]*w);
 
1102
 
 
1103
          h[2] += 0.5*(b+num[j]-1.0)*num[j]*num[j]/
 
1104
            ((1.0+num[j]*w)*(1.0+num[j]*w)) -
 
1105
            0.5*(b+num[j])*(a+sSq[j])*(a+sSq[j])*
 
1106
            num[j]*num[j]/(denorm*denorm);
 
1107
 
 
1108
          h[3] -= num[j]*num[j]*(b+num[j])*
 
1109
            (m-xBar[j])*(a+sSq[j])/(denorm*denorm);
 
1110
        }
 
1111
        break;
 
1112
 
 
1113
      case 3: // m
 
1114
        a=x[0];b=x[1];w=x[2];m=x[3];
 
1115
 
 
1116
        for(int j=0; j < numExs; j++){
 
1117
          if(Double.isNaN(xBar[j])) continue; //All missing values
 
1118
          double denorm = (1.0+num[j]*w)*(a+sSq[j]) + 
 
1119
            num[j]*(xBar[j]-m)*(xBar[j]-m);
 
1120
 
 
1121
          h[0] -= num[j]*(b+num[j])*(m-xBar[j])
 
1122
            *(1.0+num[j]*w)/(denorm*denorm);
 
1123
 
 
1124
          h[1] += num[j]*(m-xBar[j])/denorm;
 
1125
 
 
1126
          h[2] -= num[j]*num[j]*(b+num[j])*
 
1127
            (m-xBar[j])*(a+sSq[j])/(denorm*denorm);
 
1128
 
 
1129
          h[3] += num[j]*(b+num[j])*
 
1130
            ((1.0+num[j]*w)*(a+sSq[j])-
 
1131
             num[j]*(m-xBar[j])*(m-xBar[j]))
 
1132
            /(denorm*denorm);
 
1133
        }
 
1134
    }
 
1135
 
 
1136
    return h;
 
1137
  }
 
1138
}