/*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
/*
* NormalMixture.java
* Copyright (C) 2002 Yong Wang
*
*/
package weka.classifiers.functions.pace;
import java.util.Random;
import weka.core.Statistics;
/**
* Class for manipulating normal mixture distributions. <p>
*
* REFERENCES <p>
*
* Wang, Y. (2000). "A new approach to fitting linear models in high
* dimensional spaces." PhD Thesis. Department of Computer Science,
* University of Waikato, New Zealand. <p>
*
* Wang, Y. and Witten, I. H. (2002). "Modeling for optimal probability
* prediction." Proceedings of ICML'2002. Sydney. <p>
*
* @author Yong Wang (yongwang@cs.waikato.ac.nz)
* @version $Revision: 1.1.1.1 $ */
public class NormalMixture extends MixtureDistribution
{
protected double separatingThreshold = 0.05;
protected double trimingThreshold = 0.7;
protected double fittingIntervalLength = 3;
/** Contructs an empty NormalMixture
*/
public NormalMixture() {}
/** Gets the separating threshold value. This value is used by the method
separatable */
public double getSeparatingThreshold(){
return separatingThreshold;
}
/** Sets the separating threshold value
* @param t the threshold value
*/
public void setSeparatingThreshold( double t ){
separatingThreshold = t;
}
/** Gets the triming thresholding value. This value is usef by the
method trim. */
public double getTrimingThreshold(){
return trimingThreshold;
}
/** Sets the triming thresholding value. */
public void setTrimingThreshold( double t ){
trimingThreshold = t;
}
/** Return true if a value can be considered for mixture estimatino
* separately from the data indexed between i0 and i1
* @param data the data supposedly generated from the mixture
* @param i0 the index of the first element in the group
* @param i1 the index of the last element in the group
* @param x the value
*/
public boolean separable( DoubleVector data, int i0, int i1, double x )
{
double p = 0;
for( int i = i0; i <= i1; i++ ) {
p += Maths.pnorm( - Math.abs(x - data.get(i)) );
}
if( p < separatingThreshold ) return true;
else return false;
}
/** Contructs the set of support points for mixture estimation.
* @param data the data supposedly generated from the mixture
* @param ne the number of extra data that are suppposedly discarded
* earlier and not passed into here */
public DoubleVector supportPoints( DoubleVector data, int ne )
{
if( data.size() < 2 )
throw new IllegalArgumentException("data size < 2");
return data.copy();
}
/** Contructs the set of fitting intervals for mixture estimation.
* @param data the data supposedly generated from the mixture
*/
public PaceMatrix fittingIntervals( DoubleVector data )
{
DoubleVector left = data.cat( data.minus( fittingIntervalLength ) );
DoubleVector right = data.plus( fittingIntervalLength ).cat( data );
PaceMatrix a = new PaceMatrix(left.size(), 2);
a.setMatrix(0, left.size()-1, 0, left);
a.setMatrix(0, right.size()-1, 1, right);
return a;
}
/** Contructs the probability matrix for mixture estimation, given a set
* of support points and a set of intervals.
* @param s the set of support points
* @param intervals the intervals */
public PaceMatrix probabilityMatrix( DoubleVector s,
PaceMatrix intervals )
{
int ns = s.size();
int nr = intervals.getRowDimension();
PaceMatrix p = new PaceMatrix(nr, ns);
for( int i = 0; i < nr; i++ ) {
for( int j = 0; j < ns; j++ ) {
p.set( i, j,
Maths.pnorm( intervals.get(i, 1), s.get(j), 1 ) -
Maths.pnorm( intervals.get(i, 0), s.get(j), 1 ) );
}
}
return p;
}
/** Returns the empirical Bayes estimate of a single value.
* @param x the value
*/
public double empiricalBayesEstimate ( double x )
{
if( Math.abs(x) > 10 ) return x; // pratical consideration; modify later
DoubleVector d =
Maths.dnormLog( x, mixingDistribution.getPointValues(), 1 );
d.minusEquals( d.max() );
d = d.map("java.lang.Math", "exp");
d.timesEquals( mixingDistribution.getFunctionValues() );
return mixingDistribution.getPointValues().innerProduct( d ) / d.sum();
}
/** Returns the empirical Bayes estimate of a vector.
* @param x the vector
*/
public DoubleVector empiricalBayesEstimate( DoubleVector x )
{
DoubleVector pred = new DoubleVector( x.size() );
for(int i = 0; i < x.size(); i++ )
pred.set(i, empiricalBayesEstimate(x.get(i)) );
trim( pred );
return pred;
}
/** Returns the optimal nested model estimate of a vector.
* @param x the vector */
public DoubleVector nestedEstimate( DoubleVector x )
{
DoubleVector chf = new DoubleVector( x.size() );
for(int i = 0; i < x.size(); i++ ) chf.set( i, hf( x.get(i) ) );
chf.cumulateInPlace();
int index = chf.indexOfMax();
DoubleVector copy = x.copy();
if( index < x.size()-1 ) copy.set( index + 1, x.size()-1, 0 );
trim( copy );
return copy;
}
/** Returns the estimate of optimal subset selection.
* @param x the vector */
public DoubleVector subsetEstimate( DoubleVector x ) {
DoubleVector h = h( x );
DoubleVector copy = x.copy();
for( int i = 0; i < x.size(); i++ )
if( h.get(i) <= 0 ) copy.set(i, 0);
trim( copy );
return copy;
}
/** Trims the small values of the estaimte
* @param x the estimate vector */
public void trim( DoubleVector x ) {
for(int i = 0; i < x.size(); i++ ) {
if( Math.abs(x.get(i)) <= trimingThreshold ) x.set(i, 0);
}
}
/**
* Computes the value of h(x) / f(x) given the mixture. The
* implementation avoided overflow.
* @param x the value */
public double hf( double x )
{
DoubleVector points = mixingDistribution.getPointValues();
DoubleVector values = mixingDistribution.getFunctionValues();
DoubleVector d = Maths.dnormLog( x, points, 1 );
d.minusEquals( d.max() );
d = (DoubleVector) d.map("java.lang.Math", "exp");
d.timesEquals( values );
return ((DoubleVector) points.times(2*x).minusEquals(x*x))
.innerProduct( d ) / d.sum();
}
/**
* Computes the value of h(x) given the mixture.
* @param x the value */
public double h( double x )
{
DoubleVector points = mixingDistribution.getPointValues();
DoubleVector values = mixingDistribution.getFunctionValues();
DoubleVector d = (DoubleVector) Maths.dnorm( x, points, 1 ).timesEquals( values );
return ((DoubleVector) points.times(2*x).minusEquals(x*x))
.innerProduct( d );
}
/**
* Computes the value of h(x) given the mixture, where x is a vector.
* @param x the vector */
public DoubleVector h( DoubleVector x )
{
DoubleVector h = new DoubleVector( x.size() );
for( int i = 0; i < x.size(); i++ )
h.set( i, h( x.get(i) ) );
return h;
}
/**
* Computes the value of f(x) given the mixture.
* @param x the value */
public double f( double x )
{
DoubleVector points = mixingDistribution.getPointValues();
DoubleVector values = mixingDistribution.getFunctionValues();
return Maths.dchisq( x, points ).timesEquals( values ).sum();
}
/**
* Computes the value of f(x) given the mixture, where x is a vector.
* @param x the vector */
public DoubleVector f( DoubleVector x )
{
DoubleVector f = new DoubleVector( x.size() );
for( int i = 0; i < x.size(); i++ )
f.set( i, h( f.get(i) ) );
return f;
}
/** Converts to a string
*/
public String toString()
{
return mixingDistribution.toString();
}
/** Method to test this class
*/
public static void main(String args[])
{
int n1 = 50;
int n2 = 50;
double mu1 = 0;
double mu2 = 5;
DoubleVector a = Maths.rnorm( n1, mu1, 1, new Random() );
a = a.cat( Maths.rnorm( n2, mu2, 1, new Random() ) );
DoubleVector means = (new DoubleVector( n1, mu1 )).cat(new DoubleVector(n2, mu2));
System.out.println("==========================================================");
System.out.println("This is to test the estimation of the mixing\n" +
"distribution of the mixture of unit variance normal\n" +
"distributions. The example mixture used is of the form: \n\n" +
" 0.5 * N(mu1, 1) + 0.5 * N(mu2, 1)\n" );
System.out.println("It also tests three estimators: the subset\n" +
"selector, the nested model selector, and the empirical Bayes\n" +
"estimator. Quadratic losses of the estimators are given, \n" +
"and are taken as the measure of their performance.");
System.out.println("==========================================================");
System.out.println( "mu1 = " + mu1 + " mu2 = " + mu2 +"\n" );
System.out.println( a.size() + " observations are: \n\n" + a );
System.out.println( "\nQuadratic loss of the raw data (i.e., the MLE) = " +
a.sum2( means ) );
System.out.println("==========================================================");
// find the mixing distribution
NormalMixture d = new NormalMixture();
d.fit( a, NNMMethod );
System.out.println( "The estimated mixing distribution is:\n" + d );
DoubleVector pred = d.nestedEstimate( a.rev() ).rev();
System.out.println( "\nThe Nested Estimate = \n" + pred );
System.out.println( "Quadratic loss = " + pred.sum2( means ) );
pred = d.subsetEstimate( a );
System.out.println( "\nThe Subset Estimate = \n" + pred );
System.out.println( "Quadratic loss = " + pred.sum2( means ) );
pred = d.empiricalBayesEstimate( a );
System.out.println( "\nThe Empirical Bayes Estimate = \n" + pred );
System.out.println( "Quadratic loss = " + pred.sum2( means ) );
}
}