binormal.java

包含了模式识别中常用的一些分类器设计算法

/*
* file: BiNormal.java
*
* Last edited: Ryan Irwin
*
*/

/*
// These imports are not needed - Phil T. 6-23-03
// import java.awt.*;
// import java.applet.*;
// import javax.swing.*;
// import java.lang.*;
// import java.io.*;
*/
import java.util.*;
/**
 * Handles random generation of Gaussian uses
 */
public class BiNormal 
{
    
    // declare global variables
    //
    int snormi = 0;
    double snorm = 0.0;
    double snormu = 0.0;
    double snorms = 0.0;
    double ustar = 0.0;
    double snormaa = 0.0;
    double snormw = 0.0;
    double snormy = 0.0;
    double snormtt = 0.0;

    int mltmod = 0;
    int mltmoda0 = 0;
    int mltmoda1 = 0;
    int mltmodk = 0;
    int mltmodp = 0;
    int mltmodq = 0;
    int mltmodqh = 0;
    int mltmodrh = 0;
    int mltmodh = 32768;
    
    int Xm1 = 0;
    int Xm2 = 0;
    int Xa1 = 0;
    int Xa2 = 0;
    int Xa1w = 0;
    int Xa2w = 0;
    int Xa1vw = 0;
    int Xa2vw = 0;

    int Xcg1[] = new int[32];
    int Xcg2[] = new int[32];
    int Xig1[] = new int[32];
    int Xig2[] = new int[32];
    int Xlg1[] = new int[32];
    int Xlg2[] = new int[32];
    int Xqanti[] = new int[32];

    Integer info = new Integer(0);
    Integer qinit = new Integer(0);
    Integer curntg = new Integer(1);
    Integer qstate = new Integer(0);
    Random random = new Random();

     double[] a = {
	0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,
	0.2776904,0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,
	0.5791322,0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,
	0.9467818,1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,
	1.534121,1.67594,1.862732,2.153875
    };
     double[] d = {
	0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
	0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
	0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
	0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
    };
     double[] t = {
	7.673828e-4,2.30687e-3,3.860618e-3,5.438454e-3,7.0507e-3,8.708396e-3,
	1.042357e-2,1.220953e-2,1.408125e-2,1.605579e-2,1.81529e-2,2.039573e-2,
	2.281177e-2,2.543407e-2,2.830296e-2,3.146822e-2,3.499233e-2,
	3.895483e-2,4.345878e-2,4.864035e-2,5.468334e-2,6.184222e-2,
	7.047983e-2,8.113195e-2,9.462444e-2,0.1123001,0.136498,0.1716886,
	0.2276241,0.330498,0.5847031
    };
     double[] h = {
	3.920617e-2,3.932705e-2,3.951e-2,3.975703e-2,4.007093e-2,4.045533e-2,
	4.091481e-2,4.145507e-2,4.208311e-2,4.280748e-2,4.363863e-2,
	4.458932e-2,4.567523e-2,4.691571e-2,4.833487e-2,4.996298e-2,
	5.183859e-2,5.401138e-2,5.654656e-2,5.95313e-2,6.308489e-2,
	6.737503e-2,7.264544e-2,7.926471e-2,8.781922e-2,9.930398e-2,
	0.11556,0.1404344,0.1836142,0.2790016,0.7010474
    };

    /**
     * methods sets up the parameters needed to generate the multivariate 
     * normal deviates form the inputs given
     *
     *
     * @param  meanv mean vector of multivariate normal distribution
     * @param  covm covariance matrix of multivariate 
     *         normal distribution
     * @param  p dimension of the normal
     * @param  parm array of parameters needed to generate 
     *                   multivariate normal deviates
     *
     */
    public  void setgmn(double[] meanv, double[] covm, int p, 
			      double[] parm) 
    {

	int i = 0;
	int icount = 0;
	int j = 0;
	int D2 = 0;
	int D3 = 0;
	int D4 = 0;
	int D5 = 0;

	// test the input
	//
	if (p <= 0) {
	    // System.out.println("P nonpositive in SETGMN");
	    return;
	}

	parm[0] = p;


	// put p and meanv into param
	//
	for(i=2, D2=1, D3=(p+1-i+D2)/D2; D3>0; D3--, i+=D2) {
	    parm[i-1] = meanv[i-2];
	}

	// Cholesky decomposition to find A s.t. trans(A)*(A) = COVM
	//
	spofa(covm, p, p);
	if(info.intValue() != 0) {
	    // System.out.println(" COVM not positive definite in SETGMN");
	    return;
	}

	icount = p+1;

	// put upper half of a, which is now the cholesky factor, into parm
	// covm(1,1) = parm(p+2)
	// covm(1,2) = parm(p+3)
	// :
	// covm(1,p) = parm(2p+1)
	// covm(2,2) = parm(2p+2)  ...
	//
	for(i=1, D4=1, D5=(p-i+D4)/D4; D5>0; D5--, i+=D4) {
	    for(j=i-1; j<p; j++) {
		icount += 1;
		parm[icount-1] = covm[i-1+j*p];
	    }
	}
    }
    
    /**
     * methods generates the multivariate normal deviates using the procedure:
     * 1) generate p independent standard normal deviates - ei ~ n(0,1)
     * 2) using cholesky decomposition find a s.t. trans(a)*a = covm
     * 3) trans(a)e + meanv ~ n(meanv,covm)
     *
     * @param parm array of parameters needed to generate 
     *             multivariate normal deviates
     * @param x vector deviate generated
     * @param work scratch array
     *
     */
    public  void genmn(double[] parm, double[] x, double[] work) 
    {

	int i = 0;
	int j = 0;
	int p = 0;
	int D1 = 0;
	int D2 = 0;
	int D3 = 0;
	int D4 = 0;
	int icount = 0;
	double ae = 0.0;
	
	p = (int)parm[0];
	
	//  generate p independent normal deviates - work ~ n(0,1)
	//
	for(i=1; i<=p; i++) 
	{
	    work[i-1] = snorm();
	}

	for(i=1, D3=1, D4=(p-i+D3)/D3; D4>0; D4--, i+=D3) 
	{

	    // parm (p+2 : p*(p+3)/2 + 1) contains a, the cholesky
	    // decomposition of the desired covariance matrix.
	    // trans(a)(1,1) = parm(p+2)
	    // trans(a)(2,1) = parm(p+3)
	    // trans(a)(2,2) = parm(p+2+p)
	    // trans(a)(3,1) = parm(p+4)
	    // trans(a)(3,2) = parm(p+3+p)
	    // trans(a)(3,3) = parm(p+2-1+2p)  ...
	    // trans(a)*work + meanv ~ n(meanv,covm)
	    //
	    icount = 0;
	    ae = 0.0;

	    for(j=1, D1=1, D2=(i-j+D1)/D1; D2>0; D2--, j+=D1) {
		icount += (j-1);
		ae += (parm[i+(j-1)*p-icount+p]*work[j-1]);
	    }
	    x[i-1] = ae+parm[i];
	}
    }

    /**
     * generates a uniform distribution over 0 - 1
     *
     * @return  random floating point number from a uniform distribution
     *            over 0 - 1 using the current generator
     *
     */
    public double ranf() 
    {
	return random.nextDouble();
    }
    
    /**
     *
     * linpack.  this version dated 08/14/78
     * cleve moler, university of new mexico, argonne national lab
     *
     */
    public void snorm40() {
	if(ustar <= t[snormi-1]) 
	{
	    snorm60();
	    return;
	}
	snormw = (ustar-t[snormi-1])*h[snormi-1];
	snorm50();
    }

    /**
     *
     * linpack.  this version dated 08/14/78
     * cleve moler, university of new mexico, argonne national lab
     *
     */
    public void snorm50() 
    {
	snormy = snormaa+snormw;
	snorm = snormy;
	if(snorms == 1.0) snorm = -snormy;
    }

    /**
     *
     * linpack.  this version dated 08/14/78
     * cleve moler, university of new mexico, argonne national lab
     *
     */
    public void snorm60() 
    {
	snormu = ranf();
	snormw = snormu*(a[snormi]-snormaa);
	snormtt = (0.5*snormw+snormaa)*snormw;
	snorm80();
    }

    /**
     *
     * linpack.  this version dated 08/14/78
     * cleve moler, university of new mexico, argonne national lab
     *
     */
    public void snorm70() 
    {
	snormtt = snormu;
	ustar = ranf();
	snorm80();
    }

    /**
     *
     * linpack.  this version dated 08/14/78
     * cleve moler, university of new mexico, argonne national lab
     *
     */
    public void snorm80() 
    {
	if(ustar > snormtt) 
	{
	    snorm50();
	    return;
	}
	snormu = ranf();
	if(ustar >= snormu) 
	{
	    snorm70();
	    return;
	}
	ustar = ranf();
	snorm40();
    }

    /**
     *
     * linpack.  this version dated 08/14/78
     * cleve moler, university of new mexico, argonne national lab
     *
     */
    public void snorm100() 
    {
	snormi = 6;
	snormaa = a[31];
	snorm120();
    }

    /**
     *
     * linpack.  this version dated 08/14/78
     * cleve moler, university of new mexico, argonne national lab
     *
     */
    public void snorm110() 
    {
	snormaa += d[snormi-1];
	snormi += 1;
	snorm120();
    }

    /**
     *
     * linpack.  this version dated 08/14/78
     * cleve moler, university of new mexico, argonne national lab
     *
     */
    public void snorm120() 
    {
	snormu += snormu;
	if(snormu < 1.0) 
	{
	    snorm110();
	    return;
	}
	snormu -= 1.0;
	snorm140();
    }
    
    /**
     *
     * linpack.  this version dated 08/14/78
     * cleve moler, university of new mexico, argonne national lab
     *
     */
    public void snorm140() 
    {
	snormw = snormu*d[snormi-1];
	snormtt = (0.5*snormw+snormaa)*snormw;
	snorm160();
    }

    /**
     *
     * linpack.  this version dated 08/14/78
     * cleve moler, university of new mexico, argonne national lab
     *
     */
    public void snorm150() 
    {
	snormtt = snormu;
	snorm160();
    }

    /**
     *
     * linpack.  this version dated 08/14/78
     * cleve moler, university of new mexico, argonne national lab
     *
     */
    public void snorm160() 
    {
	ustar = ranf();
	if(ustar > snormtt) 
	{
	    snorm50();
	    return;
	}
	snormu = ranf();
	if(ustar >= snormu) 
	{
	    snorm150();
	    return;
	}
	snormu = ranf();
	snorm140();
    }
    
    /**
     * ahrens, j.h. and dieter, u.                            
     * extensions of forsythe's method for random             
     * sampling from the normal distribution.                 
     * math. comput., 27,124 (oct. 1973), 927 - 937.    
     *
     * @return  standard normal distribution
     *
     *
     */
    public double snorm() 
    {

	snormu = ranf();
	snorms = 0.0;
	
	if(snormu > 0.5) 
	{
	    snorms = 1.0;
	}
	
	snormu += (snormu-snorms);
	snormu = 32.0*snormu;
	snormi = (int)snormu;
	
	if(snormi == 32) 
	{
	    snormi = 31;
	}
	if(snormi == 0) 
	{
	    snorm100();
	    return snorm;
	}
	
	// start center
	//
	ustar = snormu-(double)snormi;
	snormaa = a[snormi-1];
	
	// start center
	//
	snorm40();	
	return snorm;
    }

    /**
     *
     * linpack.  this version dated 08/14/78
     * cleve moler, university of new mexico, argonne national lab
     *
     * @param n    integer
     * @param sx   array of doubles
     * @param dx   integer
     * @param incx integer
     * @param sy   array of doubles
     * @param dy   integer
     * @param incy integer
     *
     * @return double value of sdot
     */
    public  double sdot(int n, double[] sx, int dx, int incx, 
			      double[] sy, int dy, int incy) 
    {

	int i = 0;
	int ix = 0;
	int iy = 0;
	int m = 0;
	int mp1 = 0;
	double sdot = 0.0;
	double stemp = 0.0;

	stemp = sdot = 0.0;
	if(n <= 0) return sdot;

	if(incx != 1 && incy != 1) 
	{
	    
	    ix = iy = 1;
	    if(incx < 0) ix = (-n+1) * incx + 1;
	    if(incy < 0) iy = (-n+1) * incy + 1;
	    
	    for(i=1; i<=n; i++) 
	    {
		stemp += (sx[dx+ix-1] * sy[dy+iy-1]);
		ix += incx;
		iy += incy;
	    }
	    
	    sdot = stemp;
	    return sdot;
	}
	
	m = n % 5;
	
	if(m == 0) 
	{
	    mp1 = m+1;
	    for(i=mp1; i<=n; i+=5) 
	    {
		stemp += (sx[dx+i-1] * sy[dy+i-1] +
			  sx[dx+i]   * sy[dy+i]   +
			  sx[dx+i+1] * sy[dy+i+1] +
			  sx[dx+i+2] * sy[dy+i+2] +
			  sx[dx+i+3] * sy[dy+i+3]);
	    }
	    sdot = stemp;
	    return sdot;
	} 
	else 
	{
	    
	    for(i=0; i<m; i++) 
	    {
		stemp += (sx[dx+i] * sy[dy+i]);
	    }
	    if(n < 5) 
	    {
		sdot = stemp;
		return sdot;
	    } 
	    else 
	    {
		mp1 = m+1;
		for(i=mp1; i<=n; i+=5) 
		{
		    stemp += (sx[dx+i-1] * sy[dy+i-1] +
			      sx[dx+i]   * sy[dy+i]   +
			      sx[dx+i+1] * sy[dy+i+1] +
			      sx[dx+i+2] * sy[dy+i+2] +
			      sx[dx+i+3] * sy[dy+i+3]);
		}
		sdot = stemp;
		return sdot;
	    }
	}
    }
    
    /**
     * linpack.  this version dated 08/14/78
     * cleve moler, university of new mexico, argonne national lab
     *
     * @param a array of doubles
     * @param lda integer
     * @param n integer
     */
    public  void spofa(double[] a, int lda, int n) 
    {
    
	
	int j = 0;
	int jm1 = 0;
	int k = 0;
	double t = 0.0;
	double s = 0.0;
	int flag = 0;

    S40:
	for(j=1; j<=n; j++) 
	{
	    
	    info = new Integer(j);
	    s = 0.0;
	    jm1 = j-1;
	    
	    if(jm1 >= 1) 
	    {
		for(k=0; k<jm1; k++) 
		{
		    t = a[k+(j-1)*lda] - 
			sdot(k, a, k*lda, 1, a, (j-1)*lda, 1);
		    t /=  a[k+k*lda];
		    a[k+(j-1)*lda] = t;
		    s += (t*t);
		}
	    }
	    
	    s = a[j-1+(j-1)*lda] - s;
	    
	    if(s <= 0.0) 
	    {
		flag = 1;
		break S40;
	    }
	    a[j-1+(j-1)*lda] = Math.sqrt(s);
	}
	if (flag == 0) info = new Integer(0);
	
	return;
    }
    
    /**
     * Generates binormal gaussian random deviates
     *
     * @param   max number of deviates
     * @param   mx mean, X-COORDINATE
     * @param   my mean, y-coordinate
     * @param   c11 a11 element of the covariance matrix
     * @param   c12 a12 element of the covariance matrix
     * @param   c21 a21 element of the covariance matrix
     * @param   c22 a22 element of the covariance matrix
     * @param   xval x-vector of the gaussian random deviates
     * @param   yval y-vector of the gaussian random deviates
     *
     *
     */
    public  void gaussian(int max, double mx, double my, 
			  double[] xval,double[] yval, double c11,
				double c12, double c21, double c22) 
    {

	//System.out.println(c11 + "   " + c12); // 0.05 0.0
	//System.out.println(c21 + "   " + c22); // 0.0  0.05
	// declare local variables
	//
	double[] val = null;
	double[] work = null;
	double[] mean = null;
	double[] covm = null;
	double[] param = null;

	// allocate memory 
	//
	mean = new double[2];
	covm = new double[4];
	val = new double[2];
	work = new double[max];
	param = new double[max];

	// set the mean for the binormal gaussian generator
	//
	mean[0] = mx;
	mean[1] = my;

	// set the covariance matrix for the binormal gaussian generator
	//
	covm[0] = c11;
	covm[1] = c12;
	covm[2] = c21;
	covm[3] = c22;

	// initialize the binormal gaussian random generator
	//
	setgmn(mean, covm, 2, param);

	// generate the binormal gaussian random deviates
	//
	for(int i = 0; i < max; i++) 
	{
	    genmn(param, val, work);
	    xval[i] = val[0]; 
	    yval[i] = val[1];
	}
    }
}