fastmath.java

用于multivariate时间序列分类

	if (x < (a+1.0))
	    return Math.log(1-gammaSeriesExpansion(a,x));
	else 
	    return lnGammaFraction(a,x);
    }
    /**
     * @author Jaco van Kooten
     */
    private static double gammaSeriesExpansion(double a, double x) {
	double ap=a;
	double del=1.0/a;
	double sum = del;
	for (int n=1; n < MAX_ITERATIONS; n++) {
	    ++ap;
	    del *= x/ap;
	    sum += del;
	    if (del < sum*PRECISION)
		return sum*Math.exp(-x + a*Math.log(x) - logGamma(a));
	}
	return 0.0;
    }
    /**
     * @author Jaco van Kooten
     */
    private static double gammaFraction(double a, double x) {
	double b=x+1.0-a;
	double c=1.0/XMININ;
	double d=1.0/b;
	double h=d;
	double del=0.0;
	double an; 
	for (int i=1; i<MAX_ITERATIONS && Math.abs(del-1.0)>PRECISION; i++) {
	    an = -i*(i-a);
	    b += 2.0;
	    d=an*d+b;
	    c=b+an/c;
	    if (Math.abs(c) < XMININ)
		c=XMININ;
	    if (Math.abs(d) < XMININ)
		c=XMININ;
	    d=1.0/d;
	    del=d*c;
	    h *= del;
	}
	return Math.exp(-x + a*Math.log(x) - logGamma(a))*h;
    }
    private static double lnGammaFraction(double a, double x) {
	double b=x+1.0-a;
	double c=1.0/XMININ;
	double d=1.0/b;
	double h=d;
	double del=0.0;
	double an; 
	for (int i=1; i<MAX_ITERATIONS && Math.abs(del-1.0)>PRECISION; i++) {
	    an = -i*(i-a);
	    b += 2.0;
	    d=an*d+b;
	    c=b+an/c;
	    if (Math.abs(c) < XMININ)
		c=XMININ;
	    if (Math.abs(d) < XMININ)
		c=XMININ;
	    d=1.0/d;
	    del=d*c;
	    h *= del;
	}
	return -x + a*Math.log(x) - logGamma(a) + Math.log(h);
    }

    /** 
     * Returns the log of n choose k
     */
    
    public static float log2Choose(int n, int k){
        if(n==k) return 0; // log of (n choose n) = 0; 
        if(n == 0) return 0; 
        if(k == 0) return 0; 
        float retval = (float) logGamma(n) - (float) logGamma(k) - (float) logGamma(n-k); 
        retval = retval/ln2; 
        return retval; 
    }
    
    /**
     * Returns the normalCDF for the value; accurate to 10^-6 or so. 
     */
    public static float normalCDF(float val) {
        if(val < min){
            return (float) Statistics.normalProbability(val); 
        }
        if(val > max){
            return (float) Statistics.normalProbability(val); 
        }
        else {
            //Find array position: 
            float scaledVal = (val-min)/interval;
            int arrayIndex = (int) scaledVal; 
            float remainder = scaledVal - arrayIndex; 
            float lowerBnd = nCDFArray[arrayIndex]; 
            float upperBnd = nCDFArray[arrayIndex+1]; 
            return(lowerBnd + (upperBnd - lowerBnd)*remainder); 
        }
    }
    
    public static float logNormalCDF(float val) {
        if(val < min){
            return (float) Math.log(Statistics.normalProbability(val)); 
        }
        if(val > max){
            return (float) Math.log(Statistics.normalProbability(val)); 
        }
        else {
            //Find array position: 
            float scaledVal = (val-min)/interval;
            int arrayIndex = (int) scaledVal; 
            float remainder = scaledVal - arrayIndex; 
            float lowerBnd = nLCDFArray[arrayIndex]; 
            float upperBnd = nLCDFArray[arrayIndex+1]; 
            return(lowerBnd + (upperBnd - lowerBnd)*remainder); 
        }
    }
    
  
    public static void main(String args[]){
        for(float i=-4.955f; i < 5.1; i+= 0.1){
            System.out.println("In: " + i + " N: " + normalCDF(i) + " LN: " + logNormalCDF(i) + " Ex: " + Statistics.normalProbability(i) + " Err: " + Math.abs(normalCDF(i)-Statistics.normalProbability(i))); 
        }
        float maxErr = 0; 
        float maxErrPoint = 0; 
        for(float i = -4.995f; i < 5.1; i += 0.01){
            if(Math.abs(normalCDF(i)- Statistics.normalProbability(i)) > maxErr){
                maxErr = (float) Math.abs(normalCDF(i)- Statistics.normalProbability(i)); 
                maxErrPoint = i; 
            }
        }
        System.err.println("Max Error Point = " + maxErrPoint + " Max Error = " + maxErr); 
    }
    
    static final float nCDFArray[] = {
        2.866515718791934E-7f, 3.018964621715206E-7f, 3.179213654509154E-7f,
        3.3476450711587856E-7f, 3.5246589655419515E-7f, 3.710674058325394E-7f,
        3.906128516317548E-7f, 4.111480805510372E-7f, 4.3272105790841255E-7f,
        4.55381960169229E-7f, 4.791832711388164E-7f, 5.041798820600425E-7f,
        5.304291957611717E-7f, 5.579912350042749E-7f, 5.869287551894048E-7f,
        6.173073615748139E-7f, 6.49195631178801E-7f, 6.826652395340734E-7f,
        7.177910924711068E-7f, 7.546514631126327E-7f, 7.933281342672059E-7f,
        8.339065464158082E-7f, 8.764759514915567E-7f, 9.211295726589157E-7f,
        9.679647703052444E-7f, 1.0170832144641394E-6f, 1.0685910638968524E-6f,
        1.1225991520649828E-6f, 1.1792231802347991E-6f, 1.238583917960896E-6f,
        1.300807411204274E-6f, 1.3660251983477276E-6f, 1.4343745343791672E-6f,
        1.5059986235216042E-6f, 1.5810468605967448E-6f,
        1.6596750814174736E-6f, 1.7420458225132655E-6f,
        1.8283285905011901E-6f, 1.9187001414242804E-6f,
        2.0133447703882077E-6f, 2.112454611836565E-6f, 2.2162299508147E-6f,
        2.324879545581791E-6f, 2.438620961940885E-6f, 2.557680919666854E-6f,
        2.682295651422565E-6f, 2.812711274564298E-6f, 2.9491841762481335E-6f,
        3.091981412260259E-6f, 3.2413811200052373E-6f, 3.397672946097838E-6f,
        3.5611584890157243E-6f, 3.7321517572821486E-6f, 3.910979643659972E-6f,
        4.0979824158505285E-6f, 4.293514224203553E-6f, 4.497943626957144E-6f,
        4.711654133539426E-6f, 4.93504476647713E-6f, 5.1685306424694686E-6f,
        5.412543573199252E-6f, 5.667532686467238E-6f, 5.933965068249263E-6f,
        6.212326426290541E-6f, 6.5031217758652945E-6f, 6.806876148344881E-6f,
        7.124135323232339E-6f, 7.455466584336005E-6f, 7.801459500770616E-6f,
        8.162726733488658E-6f, 8.539904868061665E-6f, 8.93365527444557E-6f,
        9.344664994480821E-6f, 9.773647657893986E-6f, 1.0221344427583445E-5f,
        1.068852497498845E-5f, 1.1175988486357363E-5f, 1.1684564700747021E-5f,
        1.2215114980602643E-5f, 1.2768533415784225E-5f,
        1.3345747961922201E-5f, 1.3947721614002726E-5f,
        1.4575453616099798E-5f, 1.5229980708188631E-5f,
        1.5912378410992766E-5f, 1.662376234983364E-5f, 1.7365289618469893E-5f,
        1.8138160183931137E-5f, 1.894361833336725E-5f, 1.9782954163953555E-5f,
        2.065750511690864E-5f, 2.1568657556699246E-5f, 2.251784839652505E-5f,
        2.3506566771191853E-5f, 2.4536355758501522E-5f,
        2.5608814150302184E-5f, 2.672559827436046E-5f, 2.7888423868234892E-5f,
        2.909906800634576E-5f, 3.035937108145441E-5f, 3.167123884178135E-5f,
        3.3036644485008066E-5f, 3.44576308104261E-5f, 3.593631243051017E-5f,
        3.747487804320984E-5f, 3.9075592766269305E-5f, 4.074080053489903E-5f,
        4.247292656413995E-5f, 4.427447987727305E-5f, 4.614805590164292E-5f,
        4.809633913327755E-5f, 5.0122105871698614E-5f, 5.2228227026330666E-5f,
        5.441767099592988E-5f, 5.6693506622463423E-5f, 5.905890622088378E-5f,
        6.151714868625241E-5f, 6.407162267967406E-5f, 6.672582989451952E-5f,
        6.948338840441554E-5f, 7.234803609449349E-5f, 7.532363417739445E-5f,
        7.841417079553409E-5f, 8.162376471113726E-5f, 8.495666908555678E-5f,
        8.841727534939478E-5f, 9.201011716494704E-5f, 9.573987448249548E-5f,
        9.961137769197046E-5f, 1.0362961187151002E-4f, 1.077997211344361E-4f,
        1.1212701307617247E-4f, 1.1661696332261562E-4f,
        1.2127522018147847E-4f, 1.261076093981082E-4f, 1.311201390172805E-4f,
        1.3631900435245994E-4f, 1.4171059306400797E-4f, 1.473014903478087E-4f,
        1.5309848423576832E-4f, 1.5910857100963132E-4f, 1.653389607295399E-4f,