randomvariates.java

一个数据挖掘软件ALPHAMINERR的整个过程的JAVA版源代码

/*
 *    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.
 */

/*
 *    RandomVariates.java
 *    Copyright (C) 2002 Xin Xu
 *
 */

package weka.core;

import java.util.Random;

/**
 * Class implementing some simple random variates generator.
 *
 * @author Xin Xu ([e-mail])
 * @version $Revision$
 */
public final class RandomVariates extends Random{
    
    /** 
     * Simply the constructor of super class
     */
    public RandomVariates(){ super(); }
   
    /** 
     * Simply the constructor of super class
     *
     * @param seed the seed in this random object
     */
    public RandomVariates(long seed){ super(seed); }
    
    /** 
     * Simply use the method of the super class
     *
     * @param bits - random bits
     * @return the next pseudorandom value from this random number 
     * generator's sequence.
     */
    protected int next(int bits) {return super.next(bits);}
    
    /**
     * Generate a value of a variate following standard exponential
     * distribution using simple inverse method.<p>
     *
     * Variates related to standard Exponential can be generated using simple
     * transformations.
     * @return a value of the variate
     */ 
    public double nextExponential(){
	return -Math.log(1.0-super.nextDouble());
    }
    
    /**
     * Generate a value of a variate following standard Erlang distribution.
     * It can be used when the shape parameter is an integer and not too large,
     * say, <100.  When the parameter is not an integer (which is often named
     * Gamma distribution) or is large, use <code>nextGamma(double a)</code>
     * instead.
     *
     * @param a the shape parameter, must be no less than 1
     * @return a value of the variate
     * @exception if parameter less than 1
     */
    public double nextErlang(int a) throws Exception{
	if(a<1)
	    throw new Exception("Shape parameter of Erlang distribution must be greater than 1!");
	
	double product = 1.0;
	for(int i=1; i<=a; i++)
	    product *= super.nextDouble();
	
	return -Math.log(product);
    }
    
    /**
     * Generate a value of a variate following standard Gamma distribution 
     * with shape parameter a.<p>
     * If a>1, it uses a rejection method developed by Minh(1988)"Generating
     * Gamma Variates", ACM Trans. on Math. Software, Vol.14, No.3, pp261-266.
     * <br>
     * If a<1, it uses the algorithm "GS" developed by Ahrens and Dieter(1974)
     * "COMPUTER METHODS FOR SAMPLING FROM GAMMA, BETA, POISSON AND BINOMIAL
     * DISTRIBUTIONS", COMPUTING, 12 (1974), pp223-246, and further implemented
     * in Fortran by Ahrens, Kohrt and Dieter(1983) "Algorithm 599: sampling
     * from Gamma and Poisson distributions", ACM Trans. on Math. Software, 
     * Vol.9 No.2, pp255-257.<p> 
     * 
     * Variates related to standard Gamma can be generated using simple
     * transformations.
     *
     * @param a the shape parameter, must be greater than 1
     * @return a value of the variate
     * @exception if parameter not greater than 1
     */
    public double nextGamma(double a) throws Exception{
	if(a<=0.0)
	    throw new Exception("Shape parameter of Gamma distribution"+
				"must be greater than 0!");
	else if (a==1.0)
	    return nextExponential();
	else if (a<1.0){
	    double b=1.0+Math.exp(-1.0)*a, p,x, condition;
	    do{
		p=b*super.nextDouble();
		if(p<1.0){
		    x = Math.exp(Math.log(p)/a);
		    condition = x;
		}
		else{
		    x = -Math.log((b-p)/a);
		    condition = (1.0-a)*Math.log(x);
		}
	    }
	    while(nextExponential() < condition);
	    return x;	    
	}
	else{ // a>1
	    double b=a-1.0, D=Math.sqrt(b), D1,x1,x2,xl,f1,f2,x4,x5,xr,f4,f5,
		p1,p2,p3,p4;
	    
	    // Initialization
	    if(a<=2.0){
		D1 = b/2.0;
		x1 = 0.0;
		x2 = D1;
		xl = -1.0;
		f1 = 0.0;
	    }
	    else{
		D1 = D-0.5;
		x2 = b-D1;
		x1 = x2-D1;
		xl = 1.0-b/x1;
		f1 = Math.exp(b*Math.log(x1/b)+2.0*D1);
	    }
	    
	    f2=Math.exp(b*Math.log(x2/b)+D1);
	    x4 = b+D;
	    x5 = x4+D;
	    xr = 1.0-b/x5;
	    f4 = Math.exp(b*Math.log(x4/b)-D);
	    f5 = Math.exp(b*Math.log(x5/b)-2.0*D);
	    p1 = 2.0*f4*D;
	    p2 = 2.0*f2*D1+p1;
	    p3 = f5/xr+p2;
	    p4 = -f1/xl+p3;
	    
	    // Generation
	    double u, w=Double.MAX_VALUE, x=b, v, xp;
	    while(Math.log(w) > (b*Math.log(x/b)+b-x)){
		u=super.nextDouble()*p4;
		if(u<=p1){ // step 5-6
		    w = u/D-f4;
		    if(w<=0.0) return (b+u/f4);
		    if(w<=f5)  return (x4+(w*D)/f5);
		    
		    v = super.nextDouble();
		    x=x4+v*D;
		    xp=2.0*x4-x;
		    
		    if(w >= f4+(f4-1)*(x-x4)/(x4-b))
			return xp;
		    if(w <= f4+(b/x4-1)*f4*(x-x4))
			return x;
		    if((w < 2.0*f4-1.0) || 
		       (w < 2.0*f4-Math.exp(b*Math.log(xp/b)+b-xp)))
			continue;
		    return xp;
		}
		else if(u<=p2){ // step 7-8
		    w = (u-p1)/D1-f2;
		    if(w<=0.0) return (b-(u-p1)/f2);
		    if(w<=f1)  return (x1+w*D1/f1);
		    
		    v = super.nextDouble();
		    x=x1+v*D1;
		    xp=2.0*x2-x;
		    
		    if(w >= f2+(f2-1)*(x-x2)/(x2-b))
			return xp;
		    if(w <= f2*(x-x1)/D1)
			return x;
		    if((w < 2.0*f2-1.0) || 
		       (w < 2.0*f2-Math.exp(b*Math.log(xp/b)+b-xp)))
			continue;
		    return xp;
		}
		else if(u<p3){ // step 9
		    w = super.nextDouble();
		    u = (p3-u)/(p3-p2);
		    x = x5-Math.log(u)/xr;
		    if(w <= (xr*(x5-x)+1.0)/u) return x;
		    w = w*f5*u;
		}
		else{ // step 10
		    w = super.nextDouble();
		    u = (p4-u)/(p4-p3);
		    x = x1-Math.log(u)/xl;
		    if(x<0.0) continue;
		    if(w <= (xl*(x1-x)+1.0)/u) return x;
		    w = w*f1*u;
		}
	    }
	    
	    return x;
	}	
    }
    
    
  /**
   * Main method for testing this class.
   *
   * @param ops # of variates/seed, default is 10/45
   */
  public static void main(String[] ops) {

      int n = Integer.parseInt(ops[0]);
      if(n<=0)
	  n=10;
      long seed = Long.parseLong(ops[1]);
      if(seed <= 0)
	  seed = 45;
      RandomVariates var = new RandomVariates(seed);
      double varb[] = new double[n];
      
      try {
	  System.out.println("Generate "+n+" values with std. exp dist:");
	  for(int i=0; i<n; i++){
	      varb[i] = var.nextExponential();
	      System.out.print("["+i+"] "+varb[i]+", ");
	  }

	  System.out.println("\nMean is "+Utils.mean(varb)+
			     ", Variance is "+Utils.variance(varb)+
			     "\n\nGenerate "+n+" values with"+
			     " std. Erlang-5 dist:");

	  for(int i=0; i<n; i++){
	      varb[i] = var.nextErlang(5);
	      System.out.print("["+i+"] "+varb[i]+", ");
	  }

	  System.out.println("\nMean is "+Utils.mean(varb)+
			     ", Variance is "+Utils.variance(varb)+
			     "\n\nGenerate "+n+" values with"+
			     " std. Gamma(4.5) dist:");
	  
	  for(int i=0; i<n; i++){
	      varb[i] = var.nextGamma(4.5);
	      System.out.print("["+i+"] "+varb[i]+", ");
	  }	 
	  
	  System.out.println("\nMean is "+Utils.mean(varb)+
			     ", Variance is "+Utils.variance(varb)+
			     "\n\nGenerate "+n+" values with"+
			     " std. Gamma(0.5) dist:");
	  
	  for(int i=0; i<n; i++){
	      varb[i] = var.nextGamma(0.5);
	      System.out.print("["+i+"] "+varb[i]+", ");
	  }	  	  
	  
	  System.out.println("\nMean is "+Utils.mean(varb)+
			     ", Variance is "+Utils.variance(varb)+
			     "\n\nGenerate "+n+" values with"+
			     " std. Gaussian(5, 2) dist:");
	  
	  for(int i=0; i<n; i++){
	      varb[i] = var.nextGaussian()*2.0+5.0;
	      System.out.print("["+i+"] "+varb[i]+", ");
	  }	  	  
	  System.out.println("\nMean is "+Utils.mean(varb)+
			     ", Variance is "+Utils.variance(varb)+"\n");
	  
      } catch (Exception e) {
	  e.printStackTrace();
      }
  }
}