import cern.jet.random.Gamma; //导入依赖的package包/类
public static void dirichlet_rnd(double[] x, double[] alpha, int vSize) {
	double sum = 0.0;
	double xt = 0.0;
	int i = 0;
	
	for(i=0;i<vSize;i++) {
		if (alpha[i] <= 0.0) {
			xt = 0.0;
		} else {
			xt = Gamma.staticNextDouble(alpha[i],1.0);
		}
		if (Double.isNaN(xt)) {
			xt = 0.0;
		}
		x[i] = xt;
		sum = sum + x[i];
	}
	
	for(i=0;i<vSize;i++) {
		x[i] = x[i]/sum;
	}
}
 

import cern.jet.random.Gamma; //导入依赖的package包/类
public static double logDirichletPDF(double[] x, double[] alpha, int vSize) {
	double alpha_sum = 0.0;
	int i = 0;
	double lik = 0.0;
	
	for (i=0;i<vSize;i++) {
		if (alpha[i] > 0.0) {
			if (x[i] > 0.0) {
				lik = lik + (alpha[i] - 1.0)*Math.log(x[i]);
			}
			lik = lik - cern.jet.stat.Gamma.logGamma(alpha[i]);
			alpha_sum = alpha_sum + alpha[i];
		}
	}
	
	lik = lik+cern.jet.stat.Gamma.logGamma(alpha_sum);
	
	return lik;
}
 

import cern.jet.random.Gamma; //导入依赖的package包/类
public static double logDensity(double[] alpha, double[] x) {
	if(StatUtils.sum(x)!=1) return 0;
	double A = org.apache.commons.math3.special.
			Gamma.logGamma(StatUtils.sum(alpha));
	double B = 0;
	for(int i=0; i<alpha.length; i++) 
		B += org.apache.commons.math3.special.
			Gamma.logGamma(alpha[i]);
	double C = 0;
	for(int i=0; i<x.length; i++)
		C += (alpha[i]-1)*Math.log(x[i]);
	
	System.out.println(A);
	System.out.println(B);
	System.out.println(C);
	return A-B+C;
}
 

import cern.jet.random.Gamma; //导入依赖的package包/类
/**
 * This resets the generators
 * @param seed
 */
public static void setSeed(int seed) {
    uniform = new Random(seed);
    poisson = new Poisson(1.0, new MersenneTwister(uniform.nextInt()));
    gamma = new Gamma(1.0, 5.0, new MersenneTwister(uniform.nextInt()));
    gaussian = new Normal(0.0, 1.0, new MersenneTwister(uniform.nextInt()));
}
 

import cern.jet.random.Gamma; //导入依赖的package包/类
/**
 * NegativeBinomial (r p parameterization)
 * Parameterized as the distribution of the number of failures (X) before the rth success in independent trials
 * @param r : Number of successes 
 * @param p : Probability of success in each trial
 */
public NegativeBinomialDistrib(double r, double p){
	this.r = r;
	this.p = p;
	gamma = new Gamma(1, 1.0, new DRand());
	poisson = new Poisson(10, new DRand());
	nb = new NegativeBinomialDist(r, p);
}
 

import cern.jet.random.Gamma; //导入依赖的package包/类
public void set(double u){
	if(u==Double.POSITIVE_INFINITY){
		// throw new RuntimeException("should not set u to pos inf");
		// return;
	}
	for(int i=0; i<g.length; i++){
		if(dist[i]>0){
			if(g[i]!=null) g[i].setState(dist[i]*u, 1);
			g[i] = new Gamma(dist[i]*u, 1, re) ;
		}
	}
}
 

import cern.jet.random.Gamma; //导入依赖的package包/类
/**
 * @param k The shape parameter of the distribution, not negative or zero
 * @param theta The scale parameter of the distribution, not negative or zero
 * @param engine A uniform random number generator, not null
 */
public GammaDistribution(final double k, final double theta, final RandomEngine engine) {
  Validate.isTrue(k > 0, "k must be > 0");
  Validate.isTrue(theta > 0, "theta must be > 0");
  Validate.notNull(engine);
  _gamma = new Gamma(k, 1. / theta, engine);
  _k = k;
  _theta = theta;
}
 

import cern.jet.random.Gamma; //导入依赖的package包/类
/**
 * Constructor
 */
public GammaDistribution(double shape, double scale) {
    this.shape = shape;
    this.scale = scale;
    this.samples = 0;

    if (TRY_COLT) {
        randomEngine = new MersenneTwister(MathUtils.nextInt());
        System.out.println("Colt Gamma(" + shape + "," + scale + ")");
        coltGamma = new Gamma(shape, 1.0/scale, randomEngine);
    }
}
 

import cern.jet.random.Gamma; //导入依赖的package包/类
/**
 * NegativeBinomial (r p parameterization)
 * @param r : Number of successes 
 * @param p : Probability of success in each trial
 */
public NegativeBinomialDistrib(double r, double p){
	this.r = r;
	this.p = p;
	gamma = new Gamma(1, 1.0, new DRand());
	poisson = new Poisson(10, new DRand());
	nb = new NegativeBinomialDist(r, p);
}
 

import cern.jet.random.Gamma; //导入依赖的package包/类
public Dirichlet(double[] dist, double u){
	super(dist, u);
	this.g = new Gamma[dist.length];
	set(u);
}
 

import cern.jet.random.Gamma; //导入依赖的package包/类
public void deconv() {
    if (nogamma) {
        /* deconvolve the actual data: convolution of two gammas is a gamma of half the size
         */
        deconv = new double[smoothed.length / 2 + 1];
        for (int i = 0; i < deconv.length; i++) {
            deconv[i] = 0;
        }
        for (int i = 0; i < smoothed.length; i++) {
            deconv[i / 2] += smoothed[i];
        }
        double max = 0;
        for (int i = 0; i < deconv.length; i++) {
            if (deconv[i] > max) {
                max = deconv[i];
            }
        }
        for (int i = 0; i < deconv.length; i++) {
            deconv[i] /= max;
        }
    } else {
        double mean, var;
        int min = 0;
        if (!Double.isNaN(alpha) && !Double.isNaN(beta)) {
            mean = alpha * beta;
            var = alpha * beta * beta;
        } else {
            mean = 0;
            /* The gamma distribution isn't a good fit if it doesn't
               start at ~0.  So find the first index at which there's much
               probability mass and then effectively shift smoothed[] by
               that much
            */
            for (int i = 0; i < smoothed.length; i++) {
                if (smoothed[i] < .00000001) {
                    smoothed[i] = 0;
                    min = i;
                } else {
                    break;
                }
            }
            for (int i = min; i < smoothed.length; i++) {
                System.err.println(String.format("smoothed[%d] = %.20f", i, smoothed[i]));
                mean += smoothed[i] * (i - min);
            }
            var = 0;
            for (int i = min; i < smoothed.length; i++) {
                var += smoothed[i] * (i - min - mean) * (i - min - mean);
            }
            beta = var / mean;
            alpha = mean / beta;
        }
        System.err.println(String.format("min %d, mean %f, var %f, alpha %f, beta %f",
                                         min, mean, var, alpha, beta));
        alpha = alpha / 2;  // do the deconvolution
        min /= 2;
        deconv = new double[maxdist];
        for (int i = 0; i < min; i++) {
            deconv[i] = 0;
        }
        RandomEngine engine = RandomEngine.makeDefault();
        Gamma gamma = new Gamma(alpha, beta, engine);
        for (int i = min; i < deconv.length; i++) {
            deconv[i] = gamma.pdf(i - min);
        }            
    }
}
 

import cern.jet.random.Gamma; //导入依赖的package包/类
/**
 * Returns the hypergeometric density probability of number <tt>x</tt> 
 * when the sample size is <tt>n</tt>, the size of the positive set <tt>s</tt>
 * and the population size <tt>N</tt>.  <br>
 * It works with the gamma function and in log space in an attempt to avoid
 * numerical problems.
 * @param x # observed successes in the sample
 * @param N population size
 * @param s # of successes in population (e.g., size of positive set in the population)
 * @param n sample size
 * @return
 */
public static double hyperGeometricPDF(int x, int N, int s, int n) {
	double xd = (double) x;
	double Nd = (double) N;
	double sd = (double) s;
	double nd = (double) n;
	double kx = cern.jet.stat.Gamma.logGamma(sd+1)-cern.jet.stat.Gamma.logGamma(x+1)-cern.jet.stat.Gamma.logGamma(sd-x+1);
	double mn = cern.jet.stat.Gamma.logGamma(Nd+1)-cern.jet.stat.Gamma.logGamma(nd+1)-cern.jet.stat.Gamma.logGamma(Nd-nd+1);
	double mknx = cern.jet.stat.Gamma.logGamma(Nd-sd+1)-cern.jet.stat.Gamma.logGamma(nd-x+1)-cern.jet.stat.Gamma.logGamma(Nd-sd-(nd-x)+1);
	return Math.exp(kx + mknx - mn);
}
 

import cern.jet.random.Gamma; //导入依赖的package包/类
/**
 *
 * @return Gamma distribution RNG
 */
public static Gamma getGammaGenerator() {
    return gamma;
}