You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@mahout.apache.org by sr...@apache.org on 2009/11/25 04:31:49 UTC
svn commit: r883972 [4/10] - in /lucene/mahout/trunk:
core/src/main/java/org/apache/mahout/fpm/pfpgrowth/
matrix/src/main/java/org/apache/mahout/jet/math/
matrix/src/main/java/org/apache/mahout/jet/random/
matrix/src/main/java/org/apache/mahout/jet/ran...
Modified: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/EmpiricalWalker.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/EmpiricalWalker.java?rev=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/EmpiricalWalker.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/EmpiricalWalker.java Wed Nov 25 03:31:47 2009
@@ -37,19 +37,17 @@
* Computer Programming, Volume 2 (Seminumerical algorithms), 3rd
* edition, Addison-Wesley (1997), p120.
*
- * @author wolfgang.hoschek@cern.ch
- * @version 1.0, 09/24/99
*/
/**
* @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
*/
@Deprecated
public class EmpiricalWalker extends AbstractDiscreteDistribution {
- protected int K;
- protected int[] A;
- protected double[] F;
+ protected int K;
+ protected int[] A;
+ protected double[] F;
- protected double[] cdf; // cumulative distribution function
+ protected double[] cdf; // cumulative distribution function
/*
* James Theiler, jt@lanl.gov, the author of the GSL routine this port is based on, describes his approach as follows:
*
@@ -169,17 +167,17 @@
* @throws IllegalArgumentException if at least one of the three conditions above is violated.
*/
public EmpiricalWalker(double[] pdf, int interpolationType, RandomEngine randomGenerator) {
- setRandomGenerator(randomGenerator);
- setState(pdf,interpolationType);
- setState2(pdf);
+ setRandomGenerator(randomGenerator);
+ setState(pdf,interpolationType);
+ setState2(pdf);
}
/**
* Returns the cumulative distribution function.
*/
public double cdf(int k) {
- if (k < 0) return 0.0;
- if (k >= cdf.length-1) return 1.0;
- return cdf[k];
+ if (k < 0) return 0.0;
+ if (k >= cdf.length-1) return 1.0;
+ return cdf[k];
}
/**
* Returns a deep copy of the receiver; the copy will produce identical sequences.
@@ -188,42 +186,42 @@
* @return a copy of the receiver.
*/
public Object clone() {
- EmpiricalWalker copy = (EmpiricalWalker) super.clone();
- if (this.cdf != null) copy.cdf = (double[]) this.cdf.clone();
- if (this.A != null) copy.A = (int[]) this.A.clone();
- if (this.F != null) copy.F = (double[]) this.F.clone();
- return copy;
+ EmpiricalWalker copy = (EmpiricalWalker) super.clone();
+ if (this.cdf != null) copy.cdf = (double[]) this.cdf.clone();
+ if (this.A != null) copy.A = (int[]) this.A.clone();
+ if (this.F != null) copy.F = (double[]) this.F.clone();
+ return copy;
}
/**
* Returns a random integer <tt>k</tt> with probability <tt>pdf(k)</tt>.
*/
public int nextInt() {
- int c=0;
- double u,f;
- u = this.randomGenerator.raw();
+ int c=0;
+ double u,f;
+ u = this.randomGenerator.raw();
//#if KNUTH_CONVENTION
// c = (int)(u*(g->K));
//#else
- u *= this.K;
- c = (int)u;
- u -= c;
+ u *= this.K;
+ c = (int)u;
+ u -= c;
//#endif
- f = this.F[c];
- // fprintf(stderr,"c,f,u: %d %.4f %f\n",c,f,u);
- if (f == 1.0) return c;
- if (u < f) {
- return c;
- }
- else {
- return this.A[c];
- }
+ f = this.F[c];
+ // fprintf(stderr,"c,f,u: %d %.4f %f\n",c,f,u);
+ if (f == 1.0) return c;
+ if (u < f) {
+ return c;
+ }
+ else {
+ return this.A[c];
+ }
}
/**
* Returns the probability distribution function.
*/
public double pdf(int k) {
- if (k < 0 || k >= cdf.length-1) return 0.0;
- return cdf[k-1] - cdf[k];
+ if (k < 0 || k >= cdf.length-1) return 0.0;
+ return cdf[k-1] - cdf[k];
}
/**
* Sets the distribution parameters.
@@ -237,25 +235,25 @@
* @throws IllegalArgumentException if at least one of the three conditions above is violated.
*/
public void setState(double[] pdf, int interpolationType) {
- if (pdf==null || pdf.length==0) {
- throw new IllegalArgumentException("Non-existing pdf");
- }
-
- // compute cumulative distribution function (cdf) from probability distribution function (pdf)
- int nBins = pdf.length;
- this.cdf = new double[nBins+1];
-
- cdf[0] = 0;
- for (int i = 0; i<nBins; i++ ) {
- if (pdf[i] < 0.0) throw new IllegalArgumentException("Negative probability");
- cdf[i+1] = cdf[i] + pdf[i];
- }
- if (cdf[nBins] <= 0.0) throw new IllegalArgumentException("At leat one probability must be > 0.0");
- // now normalize to 1 (relative probabilities).
- for (int i = 0; i < nBins+1; i++ ) {
- cdf[i] /= cdf[nBins];
- }
- // cdf is now cached...
+ if (pdf==null || pdf.length==0) {
+ throw new IllegalArgumentException("Non-existing pdf");
+ }
+
+ // compute cumulative distribution function (cdf) from probability distribution function (pdf)
+ int nBins = pdf.length;
+ this.cdf = new double[nBins+1];
+
+ cdf[0] = 0;
+ for (int i = 0; i<nBins; i++ ) {
+ if (pdf[i] < 0.0) throw new IllegalArgumentException("Negative probability");
+ cdf[i+1] = cdf[i] + pdf[i];
+ }
+ if (cdf[nBins] <= 0.0) throw new IllegalArgumentException("At leat one probability must be > 0.0");
+ // now normalize to 1 (relative probabilities).
+ for (int i = 0; i < nBins+1; i++ ) {
+ cdf[i] /= cdf[nBins];
+ }
+ // cdf is now cached...
}
/**
* Sets the distribution parameters.
@@ -268,136 +266,136 @@
* @throws IllegalArgumentException if at least one of the three conditions above is violated.
*/
public void setState2(double[] pdf) {
- int size = pdf.length;
- int k,s,b;
- int nBigs, nSmalls;
- Stack Bigs;
- Stack Smalls;
- double[] E;
- double pTotal=0;
- double mean,d;
-
- //if (size < 1) {
- // throw new IllegalArgumentException("must have size greater than zero");
- //}
- /* Make sure elements of ProbArray[] are positive.
- * Won't enforce that sum is unity; instead will just normalize
- */
- for (k=0; k<size; ++k) {
- //if (pdf[k] < 0) {
- //throw new IllegalArgumentException("Probabilities must be >= 0: "+pdf[k]);
- //}
- pTotal += pdf[k];
- }
-
- /* Begin setting up the internal state */
- this.K = size;
- this.F = new double[size];
- this.A = new int[size];
-
- // normalize to relative probability
- E = new double[size];
- for (k=0; k<size; ++k) {
- E[k] = pdf[k]/pTotal;
- }
-
- /* Now create the Bigs and the Smalls */
- mean = 1.0/size;
- nSmalls=0;
- nBigs=0;
- for (k=0; k<size; ++k) {
- if (E[k] < mean) ++nSmalls;
- else ++nBigs;
- }
- Bigs = new Stack(nBigs);
- Smalls = new Stack(nSmalls);
- for (k=0; k<size; ++k) {
- if (E[k] < mean) {
- Smalls.push(k);
- }
- else {
- Bigs.push(k);
- }
- }
- /* Now work through the smalls */
- while (Smalls.size() > 0) {
- s = Smalls.pop();
- if (Bigs.size() == 0) {
- /* Then we are on our last value */
- this.A[s]=s;
- this.F[s]=1.0;
- break;
- }
- b = Bigs.pop();
- this.A[s]=b;
- this.F[s]=size*E[s];
+ int size = pdf.length;
+ int k,s,b;
+ int nBigs, nSmalls;
+ Stack Bigs;
+ Stack Smalls;
+ double[] E;
+ double pTotal=0;
+ double mean,d;
+
+ //if (size < 1) {
+ // throw new IllegalArgumentException("must have size greater than zero");
+ //}
+ /* Make sure elements of ProbArray[] are positive.
+ * Won't enforce that sum is unity; instead will just normalize
+ */
+ for (k=0; k<size; ++k) {
+ //if (pdf[k] < 0) {
+ //throw new IllegalArgumentException("Probabilities must be >= 0: "+pdf[k]);
+ //}
+ pTotal += pdf[k];
+ }
+
+ /* Begin setting up the internal state */
+ this.K = size;
+ this.F = new double[size];
+ this.A = new int[size];
+
+ // normalize to relative probability
+ E = new double[size];
+ for (k=0; k<size; ++k) {
+ E[k] = pdf[k]/pTotal;
+ }
+
+ /* Now create the Bigs and the Smalls */
+ mean = 1.0/size;
+ nSmalls=0;
+ nBigs=0;
+ for (k=0; k<size; ++k) {
+ if (E[k] < mean) ++nSmalls;
+ else ++nBigs;
+ }
+ Bigs = new Stack(nBigs);
+ Smalls = new Stack(nSmalls);
+ for (k=0; k<size; ++k) {
+ if (E[k] < mean) {
+ Smalls.push(k);
+ }
+ else {
+ Bigs.push(k);
+ }
+ }
+ /* Now work through the smalls */
+ while (Smalls.size() > 0) {
+ s = Smalls.pop();
+ if (Bigs.size() == 0) {
+ /* Then we are on our last value */
+ this.A[s]=s;
+ this.F[s]=1.0;
+ break;
+ }
+ b = Bigs.pop();
+ this.A[s]=b;
+ this.F[s]=size*E[s];
/*
#if DEBUG
- fprintf(stderr,"s=%2d, A=%2d, F=%.4f\n",s,(g->A)[s],(g->F)[s]);
+ fprintf(stderr,"s=%2d, A=%2d, F=%.4f\n",s,(g->A)[s],(g->F)[s]);
#endif
*/
- d = mean - E[s];
- E[s] += d; /* now E[s] == mean */
- E[b] -= d;
- if (E[b] < mean) {
- Smalls.push(b); /* no longer big, join ranks of the small */
- }
- else if (E[b] > mean) {
- Bigs.push(b); /* still big, put it back where you found it */
- }
- else {
- /* E[b]==mean implies it is finished too */
- this.A[b]=b;
- this.F[b]=1.0;
- }
- }
- while (Bigs.size() > 0) {
- b = Bigs.pop();
- this.A[b]=b;
- this.F[b]=1.0;
- }
- /* Stacks have been emptied, and A and F have been filled */
+ d = mean - E[s];
+ E[s] += d; /* now E[s] == mean */
+ E[b] -= d;
+ if (E[b] < mean) {
+ Smalls.push(b); /* no longer big, join ranks of the small */
+ }
+ else if (E[b] > mean) {
+ Bigs.push(b); /* still big, put it back where you found it */
+ }
+ else {
+ /* E[b]==mean implies it is finished too */
+ this.A[b]=b;
+ this.F[b]=1.0;
+ }
+ }
+ while (Bigs.size() > 0) {
+ b = Bigs.pop();
+ this.A[b]=b;
+ this.F[b]=1.0;
+ }
+ /* Stacks have been emptied, and A and F have been filled */
-
+
//#if 0
- /* if 1, then artificially set all F[k]'s to unity. This will
- * give wrong answers, but you'll get them faster. But, not
- * that much faster (I get maybe 20%); that's an upper bound
- * on what the optimal preprocessing would give.
- */
+ /* if 1, then artificially set all F[k]'s to unity. This will
+ * give wrong answers, but you'll get them faster. But, not
+ * that much faster (I get maybe 20%); that's an upper bound
+ * on what the optimal preprocessing would give.
+ */
/*
- for (k=0; k<size; ++k) {
- F[k] = 1.0;
- }
+ for (k=0; k<size; ++k) {
+ F[k] = 1.0;
+ }
//#endif
*/
//#if KNUTH_CONVENTION
- /* For convenience, set F'[k]=(k+F[k])/K */
- /* This saves some arithmetic in gsl_ran_discrete(); I find that
- * it doesn't actually make much difference.
- */
- /*
- for (k=0; k<size; ++k) {
- F[k] += k;
- F[k] /= size;
- }
+ /* For convenience, set F'[k]=(k+F[k])/K */
+ /* This saves some arithmetic in gsl_ran_discrete(); I find that
+ * it doesn't actually make much difference.
+ */
+ /*
+ for (k=0; k<size; ++k) {
+ F[k] += k;
+ F[k] /= size;
+ }
#endif
*/
- /*
- free_stack(Bigs);
- free_stack(Smalls);
- free((char *)E);
+ /*
+ free_stack(Bigs);
+ free_stack(Smalls);
+ free((char *)E);
- return g;
- */
+ return g;
+ */
}
/**
* Returns a String representation of the receiver.
*/
public String toString() {
- String interpolation = null;
- return this.getClass().getName()+"("+ ((cdf!=null) ? cdf.length : 0)+")";
+ String interpolation = null;
+ return this.getClass().getName()+"("+ ((cdf!=null) ? cdf.length : 0)+")";
}
}
Modified: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Exponential.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Exponential.java?rev=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Exponential.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Exponential.java Wed Nov 25 03:31:47 2009
@@ -20,78 +20,76 @@
* Static methods operate on a default uniform random number generator; they are synchronized.
* <p>
*
- * @author wolfgang.hoschek@cern.ch
- * @version 1.0, 09/24/99
*/
/**
* @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
*/
@Deprecated
public class Exponential extends AbstractContinousDistribution {
- protected double lambda;
+ protected double lambda;
- // The uniform random number generated shared by all <b>static</b> methods.
- protected static Exponential shared = new Exponential(1.0,makeDefaultGenerator());
+ // The uniform random number generated shared by all <b>static</b> methods.
+ protected static Exponential shared = new Exponential(1.0,makeDefaultGenerator());
/**
* Constructs a Negative Exponential distribution.
*/
public Exponential(double lambda, RandomEngine randomGenerator) {
- setRandomGenerator(randomGenerator);
- setState(lambda);
+ setRandomGenerator(randomGenerator);
+ setState(lambda);
}
/**
* Returns the cumulative distribution function.
*/
public double cdf(double x) {
- if (x <= 0.0) return 0.0;
- return 1.0 - Math.exp(-x * lambda);
+ if (x <= 0.0) return 0.0;
+ return 1.0 - Math.exp(-x * lambda);
}
/**
* Returns a random number from the distribution.
*/
public double nextDouble() {
- return nextDouble(lambda);
+ return nextDouble(lambda);
}
/**
* Returns a random number from the distribution; bypasses the internal state.
*/
public double nextDouble(double lambda) {
- return - Math.log(randomGenerator.raw()) / lambda;
+ return - Math.log(randomGenerator.raw()) / lambda;
}
/**
* Returns the probability distribution function.
*/
public double pdf(double x) {
- if (x < 0.0) return 0.0;
- return lambda*Math.exp(-x*lambda);
+ if (x < 0.0) return 0.0;
+ return lambda*Math.exp(-x*lambda);
}
/**
* Sets the mean.
*/
public void setState(double lambda) {
- this.lambda = lambda;
+ this.lambda = lambda;
}
/**
* Returns a random number from the distribution with the given lambda.
*/
public static double staticNextDouble(double lambda) {
- synchronized (shared) {
- return shared.nextDouble(lambda);
- }
+ synchronized (shared) {
+ return shared.nextDouble(lambda);
+ }
}
/**
* Returns a String representation of the receiver.
*/
public String toString() {
- return this.getClass().getName()+"("+lambda+")";
+ return this.getClass().getName()+"("+lambda+")";
}
/**
* Sets the uniform random number generated shared by all <b>static</b> methods.
* @param randomGenerator the new uniform random number generator to be shared.
*/
private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
- synchronized (shared) {
- shared.setRandomGenerator(randomGenerator);
- }
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
}
}
Modified: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/ExponentialPower.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/ExponentialPower.java?rev=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/ExponentialPower.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/ExponentialPower.java Wed Nov 25 03:31:47 2009
@@ -26,106 +26,104 @@
* L. Devroye (1986): Non-Uniform Random Variate Generation , Springer Verlag, New York.
* <p>
*
- * @author wolfgang.hoschek@cern.ch
- * @version 1.0, 09/24/99
*/
/**
* @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
*/
@Deprecated
public class ExponentialPower extends AbstractContinousDistribution {
- protected double tau;
+ protected double tau;
- // cached vars for method nextDouble(tau) (for performance only)
- private double s,sm1,tau_set = -1.0;
+ // cached vars for method nextDouble(tau) (for performance only)
+ private double s,sm1,tau_set = -1.0;
- // The uniform random number generated shared by all <b>static</b> methods.
- protected static ExponentialPower shared = new ExponentialPower(1.0,makeDefaultGenerator());
+ // The uniform random number generated shared by all <b>static</b> methods.
+ protected static ExponentialPower shared = new ExponentialPower(1.0,makeDefaultGenerator());
/**
* Constructs an Exponential Power distribution.
* Example: tau=1.0.
* @throws IllegalArgumentException if <tt>tau < 1.0</tt>.
*/
public ExponentialPower(double tau, RandomEngine randomGenerator) {
- setRandomGenerator(randomGenerator);
- setState(tau);
+ setRandomGenerator(randomGenerator);
+ setState(tau);
}
/**
* Returns a random number from the distribution.
*/
public double nextDouble() {
- return nextDouble(this.tau);
+ return nextDouble(this.tau);
}
/**
* Returns a random number from the distribution; bypasses the internal state.
* @throws IllegalArgumentException if <tt>tau < 1.0</tt>.
*/
public double nextDouble(double tau) {
- double u,u1,v,x,y;
+ double u,u1,v,x,y;
- if (tau != tau_set) { // SET-UP
- s = 1.0/tau;
- sm1 = 1.0 - s;
-
- tau_set = tau;
- }
-
- // GENERATOR
- do {
- u = randomGenerator.raw(); // U(0/1)
- u = (2.0*u) - 1.0; // U(-1.0/1.0)
- u1 = Math.abs(u); // u1=|u|
- v = randomGenerator.raw(); // U(0/1)
-
- if (u1 <= sm1) { // Uniform hat-function for x <= (1-1/tau)
- x = u1;
- }
- else { // Exponential hat-function for x > (1-1/tau)
- y = tau*(1.0 - u1); // U(0/1)
- x = sm1 - s*Math.log(y);
- v = v*y;
- }
- }
-
- // Acceptance/Rejection
- while (Math.log(v) > -Math.exp(Math.log(x)*tau));
-
- // Random sign
- if (u < 0.0)
- return x;
- else
- return -x;
+ if (tau != tau_set) { // SET-UP
+ s = 1.0/tau;
+ sm1 = 1.0 - s;
+
+ tau_set = tau;
+ }
+
+ // GENERATOR
+ do {
+ u = randomGenerator.raw(); // U(0/1)
+ u = (2.0*u) - 1.0; // U(-1.0/1.0)
+ u1 = Math.abs(u); // u1=|u|
+ v = randomGenerator.raw(); // U(0/1)
+
+ if (u1 <= sm1) { // Uniform hat-function for x <= (1-1/tau)
+ x = u1;
+ }
+ else { // Exponential hat-function for x > (1-1/tau)
+ y = tau*(1.0 - u1); // U(0/1)
+ x = sm1 - s*Math.log(y);
+ v = v*y;
+ }
+ }
+
+ // Acceptance/Rejection
+ while (Math.log(v) > -Math.exp(Math.log(x)*tau));
+
+ // Random sign
+ if (u < 0.0)
+ return x;
+ else
+ return -x;
}
/**
* Sets the distribution parameter.
* @throws IllegalArgumentException if <tt>tau < 1.0</tt>.
*/
public void setState(double tau) {
- if (tau<1.0) throw new IllegalArgumentException();
- this.tau = tau;
+ if (tau<1.0) throw new IllegalArgumentException();
+ this.tau = tau;
}
/**
* Returns a random number from the distribution.
* @throws IllegalArgumentException if <tt>tau < 1.0</tt>.
*/
public static double staticNextDouble(double tau) {
- synchronized (shared) {
- return shared.nextDouble(tau);
- }
+ synchronized (shared) {
+ return shared.nextDouble(tau);
+ }
}
/**
* Returns a String representation of the receiver.
*/
public String toString() {
- return this.getClass().getName()+"("+tau+")";
+ return this.getClass().getName()+"("+tau+")";
}
/**
* Sets the uniform random number generated shared by all <b>static</b> methods.
* @param randomGenerator the new uniform random number generator to be shared.
*/
private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
- synchronized (shared) {
- shared.setRandomGenerator(randomGenerator);
- }
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
}
}
Modified: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Fun.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Fun.java?rev=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Fun.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Fun.java Wed Nov 25 03:31:47 2009
@@ -14,289 +14,287 @@
* <b>Implementation:</b> High performance implementation.
* <dt>This is a port of <tt>gen_fun.cpp</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
*
- * @author wolfgang.hoschek@cern.ch
- * @version 1.0, 09/24/99
*/
class Fun {
/**
* Makes this class non instantiable, but still let's others inherit from it.
*/
protected Fun() {
- throw new RuntimeException("Non instantiable");
+ throw new RuntimeException("Non instantiable");
}
private static double _fkt_value(double lambda, double z1, double z2, double x_value) {
- double y_value;
+ double y_value;
- y_value = Math.cos(z1*x_value)/(Math.pow((x_value*x_value + z2*z2),(lambda + 0.5)));
- return y_value;
+ y_value = Math.cos(z1*x_value)/(Math.pow((x_value*x_value + z2*z2),(lambda + 0.5)));
+ return y_value;
}
public static double bessel2_fkt(double lambda, double beta) {
- final double pi = Math.PI;
- double sum,x,step,x1,first_value,new_value;
- double epsilon = 0.01;
- double y,fx,z1,z2,erg;
- double period,border,first_sum,second_sum;
- double my,c,prod=0.0,diff,value;
- int i,j,nr_per;
-
- final double b0[] = { -1.5787132, -0.6130827, 0.1735823, 1.4793411,
- 2.6667307, 4.9086836, 8.1355339,
- };
- final double b05[] = { -1.9694802, -0.7642538, 0.0826017, 1.4276355,
- 2.6303682, 4.8857787, 8.1207968,
- };
- final double b1[] = { -2.9807345, -1.1969943, -0.1843161, 1.2739241,
- 2.5218256, 4.8172216, 8.0765633,
- };
- final double b2[] = { -5.9889676, -2.7145389, -1.1781269, 0.6782201,
- 2.0954009, 4.5452152, 7.9003173,
- };
- final double b3[] = { -9.6803440, -4.8211925, -2.6533185, -0.2583337,
- 1.4091915, 4.0993448, 7.6088310,
- };
- final double b5[] = {-18.1567152,-10.0939408, -6.5819139, -2.9371545,
- -0.6289005, 2.7270412, 6.6936799,
- };
- final double b8[] = {-32.4910195,-19.6065943,-14.0347298, -8.3839439,
- -4.9679730, -0.3567823, 4.5589697,
- };
-
- if (lambda == 0.0) {
- if (beta == 0.1) return(b0[0]);
- if (beta == 0.5) return(b0[1]);
- if (beta == 1.0) return(b0[2]);
- if (beta == 2.0) return(b0[3]);
- if (beta == 3.0) return(b0[4]);
- if (beta == 5.0) return(b0[5]);
- if (beta == 8.0) return(b0[6]);
- }
-
- if (lambda == 0.5) {
- if (beta == 0.1) return(b05[0]);
- if (beta == 0.5) return(b05[1]);
- if (beta == 1.0) return(b05[2]);
- if (beta == 2.0) return(b05[3]);
- if (beta == 3.0) return(b05[4]);
- if (beta == 5.0) return(b05[5]);
- if (beta == 8.0) return(b05[6]);
- }
-
- if (lambda == 1.0) {
- if (beta == 0.1) return(b1[0]);
- if (beta == 0.5) return(b1[1]);
- if (beta == 1.0) return(b1[2]);
- if (beta == 2.0) return(b1[3]);
- if (beta == 3.0) return(b1[4]);
- if (beta == 5.0) return(b1[5]);
- if (beta == 8.0) return(b1[6]);
- }
-
- if (lambda == 2.0) {
- if (beta == 0.1) return(b2[0]);
- if (beta == 0.5) return(b2[1]);
- if (beta == 1.0) return(b2[2]);
- if (beta == 2.0) return(b2[3]);
- if (beta == 3.0) return(b2[4]);
- if (beta == 5.0) return(b2[5]);
- if (beta == 8.0) return(b2[6]);
- }
-
- if (lambda == 3.0) {
- if (beta == 0.1) return(b3[0]);
- if (beta == 0.5) return(b3[1]);
- if (beta == 1.0) return(b3[2]);
- if (beta == 2.0) return(b3[3]);
- if (beta == 3.0) return(b3[4]);
- if (beta == 5.0) return(b3[5]);
- if (beta == 8.0) return(b3[6]);
- }
-
- if (lambda == 5.0) {
- if (beta == 0.1) return(b5[0]);
- if (beta == 0.5) return(b5[1]);
- if (beta == 1.0) return(b5[2]);
- if (beta == 2.0) return(b5[3]);
- if (beta == 3.0) return(b5[4]);
- if (beta == 5.0) return(b5[5]);
- if (beta == 8.0) return(b5[6]);
- }
-
- if (lambda == 8.0) {
- if (beta == 0.1) return(b8[0]);
- if (beta == 0.5) return(b8[1]);
- if (beta == 1.0) return(b8[2]);
- if (beta == 2.0) return(b8[3]);
- if (beta == 3.0) return(b8[4]);
- if (beta == 5.0) return(b8[5]);
- if (beta == 8.0) return(b8[6]);
- }
-
-
- if ((beta - 5.0*lambda - 8.0) >= 0.0) {
- my = 4.0*lambda*lambda;
- c = -0.9189385 + 0.5*Math.log(beta) + beta;
- sum = 1.0;
- value = 1.0;
- diff = 8.0;
- i = 1;
- for (;;) { //while (!NULL) {
- if ((factorial(i)*Math.pow((8.0*beta),i)) > 1.0e250) break;
- if (i > 10) break;
- if (i == 1) prod = my - 1.0;
- else {
- value += diff;
- prod = prod*(my - value);
- diff *= 2.0;
- }
- sum =sum + prod/(factorial(i)*Math.pow((8.0*beta),i));
- i++;
- }
- erg = c - Math.log(sum);
- return(erg);
- }
-
- if ((lambda > 0.0) && ((beta - 0.04*lambda) <= 0.0)) {
- if (lambda < 11.5) {
- erg = -Math.log(gamma(lambda)) - lambda*Math.log(2.0) + lambda*Math.log(beta);
- return(erg);
- }
- else {
- erg = -(lambda + 1.0)*Math.log(2.0) - (lambda - 0.5)*Math.log(lambda) + lambda + lambda*Math.log(beta) - 0.5*Math.log(0.5*pi);
- return(erg);
- }
- }
-
-
- // otherwise numerical integration of the function defined above
-
- x = 0.0;
-
- if (beta < 1.57) {
- fx = (fkt2_value(lambda,beta,x))*0.01;
- y = 0.0;
- for (;;) { //while (!NULL) {
- y += 0.1;
- if ((fkt2_value(lambda,beta,y)) < fx) break;
- }
- step = y*0.001;
- x1 = step;
- sum = (0.5*(10.0*step + fkt2_value(lambda,beta,x1)))*step;
- first_value = sum;
- for (;;) { //while (!NULL) {
- x = x1;
- x1 += step;
- new_value = (0.5*(fkt2_value(lambda,beta,x) + fkt2_value(lambda,beta,x1)))*step;
- sum += new_value;
- if ((new_value/first_value) < epsilon) break;
- }
- erg = -Math.log(2.0*sum);
- return(erg);
- }
- else {
- z2 = 1.57;
- z1 = beta/1.57;
- sum = 0.0;
- period = pi/z1;
- step = 0.1*period;
- border = 100.0/((lambda + 0.1)*(lambda + 0.1));
- nr_per = (int) Math.ceil((border/period)) + 20;
- x1 = step;
- for (i = 1; i <= nr_per; i++) {
- for (j = 1; j <= 10; j++) {
- new_value = (0.5*(_fkt_value(lambda,z1,z2,x) + _fkt_value(lambda,z1,z2,x1)))*step;
- sum += new_value;
- x = x1;
- x1 += step;
- }
- }
- for (j = 1; j <= 5; j++) {
- new_value = (0.5*(_fkt_value(lambda,z1,z2,x) + _fkt_value(lambda,z1,z2,x1)))*step;
- sum += new_value;
- x = x1;
- x1 += step;
- }
- first_sum = sum;
- for (j = 1; j <= 10; j++) {
- new_value = (0.5*(_fkt_value(lambda,z1,z2,x) + _fkt_value(lambda,z1,z2,x1)))*step;
- sum += new_value;
- x = x1;
- x1 += step;
- }
- second_sum = sum;
- sum = 0.5*(first_sum + second_sum);
- erg = gamma(lambda + 0.5)*Math.pow((2.0*z2),lambda)/(Math.sqrt(pi)*Math.pow(z1,lambda))*sum;
- erg = -Math.log(2.0*erg);
- return(erg);
- }
+ final double pi = Math.PI;
+ double sum,x,step,x1,first_value,new_value;
+ double epsilon = 0.01;
+ double y,fx,z1,z2,erg;
+ double period,border,first_sum,second_sum;
+ double my,c,prod=0.0,diff,value;
+ int i,j,nr_per;
+
+ final double b0[] = { -1.5787132, -0.6130827, 0.1735823, 1.4793411,
+ 2.6667307, 4.9086836, 8.1355339,
+ };
+ final double b05[] = { -1.9694802, -0.7642538, 0.0826017, 1.4276355,
+ 2.6303682, 4.8857787, 8.1207968,
+ };
+ final double b1[] = { -2.9807345, -1.1969943, -0.1843161, 1.2739241,
+ 2.5218256, 4.8172216, 8.0765633,
+ };
+ final double b2[] = { -5.9889676, -2.7145389, -1.1781269, 0.6782201,
+ 2.0954009, 4.5452152, 7.9003173,
+ };
+ final double b3[] = { -9.6803440, -4.8211925, -2.6533185, -0.2583337,
+ 1.4091915, 4.0993448, 7.6088310,
+ };
+ final double b5[] = {-18.1567152,-10.0939408, -6.5819139, -2.9371545,
+ -0.6289005, 2.7270412, 6.6936799,
+ };
+ final double b8[] = {-32.4910195,-19.6065943,-14.0347298, -8.3839439,
+ -4.9679730, -0.3567823, 4.5589697,
+ };
+
+ if (lambda == 0.0) {
+ if (beta == 0.1) return(b0[0]);
+ if (beta == 0.5) return(b0[1]);
+ if (beta == 1.0) return(b0[2]);
+ if (beta == 2.0) return(b0[3]);
+ if (beta == 3.0) return(b0[4]);
+ if (beta == 5.0) return(b0[5]);
+ if (beta == 8.0) return(b0[6]);
+ }
+
+ if (lambda == 0.5) {
+ if (beta == 0.1) return(b05[0]);
+ if (beta == 0.5) return(b05[1]);
+ if (beta == 1.0) return(b05[2]);
+ if (beta == 2.0) return(b05[3]);
+ if (beta == 3.0) return(b05[4]);
+ if (beta == 5.0) return(b05[5]);
+ if (beta == 8.0) return(b05[6]);
+ }
+
+ if (lambda == 1.0) {
+ if (beta == 0.1) return(b1[0]);
+ if (beta == 0.5) return(b1[1]);
+ if (beta == 1.0) return(b1[2]);
+ if (beta == 2.0) return(b1[3]);
+ if (beta == 3.0) return(b1[4]);
+ if (beta == 5.0) return(b1[5]);
+ if (beta == 8.0) return(b1[6]);
+ }
+
+ if (lambda == 2.0) {
+ if (beta == 0.1) return(b2[0]);
+ if (beta == 0.5) return(b2[1]);
+ if (beta == 1.0) return(b2[2]);
+ if (beta == 2.0) return(b2[3]);
+ if (beta == 3.0) return(b2[4]);
+ if (beta == 5.0) return(b2[5]);
+ if (beta == 8.0) return(b2[6]);
+ }
+
+ if (lambda == 3.0) {
+ if (beta == 0.1) return(b3[0]);
+ if (beta == 0.5) return(b3[1]);
+ if (beta == 1.0) return(b3[2]);
+ if (beta == 2.0) return(b3[3]);
+ if (beta == 3.0) return(b3[4]);
+ if (beta == 5.0) return(b3[5]);
+ if (beta == 8.0) return(b3[6]);
+ }
+
+ if (lambda == 5.0) {
+ if (beta == 0.1) return(b5[0]);
+ if (beta == 0.5) return(b5[1]);
+ if (beta == 1.0) return(b5[2]);
+ if (beta == 2.0) return(b5[3]);
+ if (beta == 3.0) return(b5[4]);
+ if (beta == 5.0) return(b5[5]);
+ if (beta == 8.0) return(b5[6]);
+ }
+
+ if (lambda == 8.0) {
+ if (beta == 0.1) return(b8[0]);
+ if (beta == 0.5) return(b8[1]);
+ if (beta == 1.0) return(b8[2]);
+ if (beta == 2.0) return(b8[3]);
+ if (beta == 3.0) return(b8[4]);
+ if (beta == 5.0) return(b8[5]);
+ if (beta == 8.0) return(b8[6]);
+ }
+
+
+ if ((beta - 5.0*lambda - 8.0) >= 0.0) {
+ my = 4.0*lambda*lambda;
+ c = -0.9189385 + 0.5*Math.log(beta) + beta;
+ sum = 1.0;
+ value = 1.0;
+ diff = 8.0;
+ i = 1;
+ for (;;) { //while (!NULL) {
+ if ((factorial(i)*Math.pow((8.0*beta),i)) > 1.0e250) break;
+ if (i > 10) break;
+ if (i == 1) prod = my - 1.0;
+ else {
+ value += diff;
+ prod = prod*(my - value);
+ diff *= 2.0;
+ }
+ sum =sum + prod/(factorial(i)*Math.pow((8.0*beta),i));
+ i++;
+ }
+ erg = c - Math.log(sum);
+ return(erg);
+ }
+
+ if ((lambda > 0.0) && ((beta - 0.04*lambda) <= 0.0)) {
+ if (lambda < 11.5) {
+ erg = -Math.log(gamma(lambda)) - lambda*Math.log(2.0) + lambda*Math.log(beta);
+ return(erg);
+ }
+ else {
+ erg = -(lambda + 1.0)*Math.log(2.0) - (lambda - 0.5)*Math.log(lambda) + lambda + lambda*Math.log(beta) - 0.5*Math.log(0.5*pi);
+ return(erg);
+ }
+ }
+
+
+ // otherwise numerical integration of the function defined above
+
+ x = 0.0;
+
+ if (beta < 1.57) {
+ fx = (fkt2_value(lambda,beta,x))*0.01;
+ y = 0.0;
+ for (;;) { //while (!NULL) {
+ y += 0.1;
+ if ((fkt2_value(lambda,beta,y)) < fx) break;
+ }
+ step = y*0.001;
+ x1 = step;
+ sum = (0.5*(10.0*step + fkt2_value(lambda,beta,x1)))*step;
+ first_value = sum;
+ for (;;) { //while (!NULL) {
+ x = x1;
+ x1 += step;
+ new_value = (0.5*(fkt2_value(lambda,beta,x) + fkt2_value(lambda,beta,x1)))*step;
+ sum += new_value;
+ if ((new_value/first_value) < epsilon) break;
+ }
+ erg = -Math.log(2.0*sum);
+ return(erg);
+ }
+ else {
+ z2 = 1.57;
+ z1 = beta/1.57;
+ sum = 0.0;
+ period = pi/z1;
+ step = 0.1*period;
+ border = 100.0/((lambda + 0.1)*(lambda + 0.1));
+ nr_per = (int) Math.ceil((border/period)) + 20;
+ x1 = step;
+ for (i = 1; i <= nr_per; i++) {
+ for (j = 1; j <= 10; j++) {
+ new_value = (0.5*(_fkt_value(lambda,z1,z2,x) + _fkt_value(lambda,z1,z2,x1)))*step;
+ sum += new_value;
+ x = x1;
+ x1 += step;
+ }
+ }
+ for (j = 1; j <= 5; j++) {
+ new_value = (0.5*(_fkt_value(lambda,z1,z2,x) + _fkt_value(lambda,z1,z2,x1)))*step;
+ sum += new_value;
+ x = x1;
+ x1 += step;
+ }
+ first_sum = sum;
+ for (j = 1; j <= 10; j++) {
+ new_value = (0.5*(_fkt_value(lambda,z1,z2,x) + _fkt_value(lambda,z1,z2,x1)))*step;
+ sum += new_value;
+ x = x1;
+ x1 += step;
+ }
+ second_sum = sum;
+ sum = 0.5*(first_sum + second_sum);
+ erg = gamma(lambda + 0.5)*Math.pow((2.0*z2),lambda)/(Math.sqrt(pi)*Math.pow(z1,lambda))*sum;
+ erg = -Math.log(2.0*erg);
+ return(erg);
+ }
}
/**
* Modified Bessel Functions of First Kind - Order 0.
*/
public static double bessi0(double x) {
- double ax,ans;
- double y;
+ double ax,ans;
+ double y;
- if ((ax=Math.abs(x)) < 3.75) {
- y=x/3.75;
- y*=y;
- ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
- +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
- }
- else {
- y=3.75/ax;
- ans=(Math.exp(ax)/Math.sqrt(ax))*(0.39894228+y*(0.1328592e-1
- +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
- +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
- +y*0.392377e-2))))))));
- }
- return ans;
+ if ((ax=Math.abs(x)) < 3.75) {
+ y=x/3.75;
+ y*=y;
+ ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
+ +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
+ }
+ else {
+ y=3.75/ax;
+ ans=(Math.exp(ax)/Math.sqrt(ax))*(0.39894228+y*(0.1328592e-1
+ +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
+ +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
+ +y*0.392377e-2))))))));
+ }
+ return ans;
}
/**
* Modified Bessel Functions of First Kind - Order 1.
*/
public static double bessi1(double x) {
- double ax,ans;
- double y;
+ double ax,ans;
+ double y;
- if ((ax=Math.abs(x)) < 3.75) {
- y=x/3.75;
- y*=y;
- ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934
- +y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
- }
- else {
- y=3.75/ax;
- ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1
- -y*0.420059e-2));
- ans=0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2
- +y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
- ans *= (Math.exp(ax)/Math.sqrt(ax));
- }
- return x < 0.0 ? -ans : ans;
+ if ((ax=Math.abs(x)) < 3.75) {
+ y=x/3.75;
+ y*=y;
+ ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934
+ +y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
+ }
+ else {
+ y=3.75/ax;
+ ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1
+ -y*0.420059e-2));
+ ans=0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2
+ +y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
+ ans *= (Math.exp(ax)/Math.sqrt(ax));
+ }
+ return x < 0.0 ? -ans : ans;
}
/**
* Returns <tt>n!</tt>.
*/
public static long factorial(int n) {
- return org.apache.mahout.jet.math.Arithmetic.longFactorial(n);
- /*
- long i,prod;
-
- prod = 1;
- if (n != 0) {
- for (i = 2; i <= n; i++) prod *= i;
- }
- return prod;
- */
+ return org.apache.mahout.jet.math.Arithmetic.longFactorial(n);
+ /*
+ long i,prod;
+
+ prod = 1;
+ if (n != 0) {
+ for (i = 2; i <= n; i++) prod *= i;
+ }
+ return prod;
+ */
}
private static double fkt2_value(double lambda, double beta, double x_value) {
- double y_value;
+ double y_value;
- y_value = cosh(lambda*x_value)*Math.exp(-beta*cosh(x_value));
- return y_value;
+ y_value = cosh(lambda*x_value)*Math.exp(-beta*cosh(x_value));
+ return y_value;
}
private static double cosh(double x) {
- return (Math.exp(x) + Math.exp(-x)) / 2.0;
+ return (Math.exp(x) + Math.exp(-x)) / 2.0;
}
@@ -304,28 +302,28 @@
* Returns the gamma function <tt>gamma(x)</tt>.
*/
public static double gamma(double x) {
- x = logGamma(x);
- //if (x > Math.log(Double.MAX_VALUE)) return Double.MAX_VALUE;
- return Math.exp(x);
+ x = logGamma(x);
+ //if (x > Math.log(Double.MAX_VALUE)) return Double.MAX_VALUE;
+ return Math.exp(x);
}
/**
* Returns a quick approximation of <tt>log(gamma(x))</tt>.
*/
public static double logGamma(double x) {
- final double c0 = 9.1893853320467274e-01, c1 = 8.3333333333333333e-02,
- c2 =-2.7777777777777777e-03, c3 = 7.9365079365079365e-04,
- c4 =-5.9523809523809524e-04, c5 = 8.4175084175084175e-04,
- c6 =-1.9175269175269175e-03;
- double g,r,z;
-
- if (x <= 0.0 /* || x > 1.3e19 */ ) return -999;
-
- for (z = 1.0; x < 11.0; x++) z *= x;
-
- r = 1.0 / (x * x);
- g = c1 + r * (c2 + r * (c3 + r * (c4 + r * (c5 + r + c6))));
- g = (x - 0.5) * Math.log(x) - x + c0 + g / x;
- if (z == 1.0) return g;
- return g - Math.log(z);
+ final double c0 = 9.1893853320467274e-01, c1 = 8.3333333333333333e-02,
+ c2 =-2.7777777777777777e-03, c3 = 7.9365079365079365e-04,
+ c4 =-5.9523809523809524e-04, c5 = 8.4175084175084175e-04,
+ c6 =-1.9175269175269175e-03;
+ double g,r,z;
+
+ if (x <= 0.0 /* || x > 1.3e19 */ ) return -999;
+
+ for (z = 1.0; x < 11.0; x++) z *= x;
+
+ r = 1.0 / (x * x);
+ g = c1 + r * (c2 + r * (c3 + r * (c4 + r * (c5 + r + c6))));
+ g = (x - 0.5) * Math.log(x) - x + c0 + g / x;
+ if (z == 1.0) return g;
+ return g - Math.log(z);
}
}
Modified: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Gamma.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Gamma.java?rev=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Gamma.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Gamma.java Wed Nov 25 03:31:47 2009
@@ -44,39 +44,37 @@
* J.H. Ahrens, U. Dieter (1982): Generating gamma variates by a modified rejection technique,
* Communications of the ACM 25, 47-54.
*
- * @author wolfgang.hoschek@cern.ch
- * @version 1.0, 09/24/99
*/
/**
* @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
*/
@Deprecated
public class Gamma extends AbstractContinousDistribution {
- protected double alpha;
- protected double lambda;
+ protected double alpha;
+ protected double lambda;
- // The uniform random number generated shared by all <b>static</b> methods.
- protected static Gamma shared = new Gamma(1.0,1.0,makeDefaultGenerator());
+ // The uniform random number generated shared by all <b>static</b> methods.
+ protected static Gamma shared = new Gamma(1.0,1.0,makeDefaultGenerator());
/**
* Constructs a Gamma distribution.
* Example: alpha=1.0, lambda=1.0.
* @throws IllegalArgumentException if <tt>alpha <= 0.0 || lambda <= 0.0</tt>.
*/
public Gamma(double alpha, double lambda, RandomEngine randomGenerator) {
- setRandomGenerator(randomGenerator);
- setState(alpha,lambda);
+ setRandomGenerator(randomGenerator);
+ setState(alpha,lambda);
}
/**
* Returns the cumulative distribution function.
*/
public double cdf(double x) {
- return Probability.gamma(alpha,lambda,x);
+ return Probability.gamma(alpha,lambda,x);
}
/**
* Returns a random number from the distribution.
*/
public double nextDouble() {
- return nextDouble(alpha, lambda);
+ return nextDouble(alpha, lambda);
}
/**
* Returns a random number from the distribution; bypasses the internal state.
@@ -104,174 +102,174 @@
* - NORMAL(seed) ... Normal generator N(0,1). *
* *
******************************************************************/
- double a = alpha;
- double aa = -1.0, aaa = -1.0,
- b=0.0, c=0.0, d=0.0, e, r, s=0.0, si=0.0, ss=0.0, q0=0.0,
- q1 = 0.0416666664, q2 = 0.0208333723, q3 = 0.0079849875,
- q4 = 0.0015746717, q5 = -0.0003349403, q6 = 0.0003340332,
- q7 = 0.0006053049, q8 = -0.0004701849, q9 = 0.0001710320,
- a1 = 0.333333333, a2 = -0.249999949, a3 = 0.199999867,
- a4 =-0.166677482, a5 = 0.142873973, a6 =-0.124385581,
- a7 = 0.110368310, a8 = -0.112750886, a9 = 0.104089866,
- e1 = 1.000000000, e2 = 0.499999994, e3 = 0.166666848,
- e4 = 0.041664508, e5 = 0.008345522, e6 = 0.001353826,
- e7 = 0.000247453;
-
- double gds,p,q,t,sign_u,u,v,w,x;
- double v1,v2,v12;
-
- // Check for invalid input values
-
- if (a <= 0.0) throw new IllegalArgumentException();
- if (lambda <= 0.0) new IllegalArgumentException();
-
- if (a < 1.0) { // CASE A: Acceptance rejection algorithm gs
- b = 1.0 + 0.36788794412 * a; // Step 1
- for(;;) {
- p = b * randomGenerator.raw();
- if (p <= 1.0) { // Step 2. Case gds <= 1
- gds = Math.exp(Math.log(p) / a);
- if (Math.log(randomGenerator.raw()) <= -gds) return(gds/lambda);
- }
- else { // Step 3. Case gds > 1
- gds = - Math.log ((b - p) / a);
- if (Math.log(randomGenerator.raw()) <= ((a - 1.0) * Math.log(gds))) return(gds/lambda);
- }
- }
- }
-
- else { // CASE B: Acceptance complement algorithm gd (gaussian distribution, box muller transformation)
- if (a != aa) { // Step 1. Preparations
- aa = a;
- ss = a - 0.5;
- s = Math.sqrt(ss);
- d = 5.656854249 - 12.0 * s;
- }
- // Step 2. Normal deviate
- do {
- v1 = 2.0 * randomGenerator.raw() - 1.0;
- v2 = 2.0 * randomGenerator.raw() - 1.0;
- v12 = v1*v1 + v2*v2;
- } while ( v12 > 1.0 );
- t = v1*Math.sqrt(-2.0*Math.log(v12)/v12);
- x = s + 0.5 * t;
- gds = x * x;
- if (t >= 0.0) return(gds/lambda); // Immediate acceptance
-
- u = randomGenerator.raw(); // Step 3. Uniform random number
- if (d * u <= t * t * t) return(gds/lambda); // Squeeze acceptance
-
- if (a != aaa) { // Step 4. Set-up for hat case
- aaa = a;
- r = 1.0 / a;
- q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) *
- r + q3) * r + q2) * r + q1) * r;
- if (a > 3.686) {
- if (a > 13.022) {
- b = 1.77;
- si = 0.75;
- c = 0.1515 / s;
- }
- else {
- b = 1.654 + 0.0076 * ss;
- si = 1.68 / s + 0.275;
- c = 0.062 / s + 0.024;
- }
- }
- else {
- b = 0.463 + s - 0.178 * ss;
- si = 1.235;
- c = 0.195 / s - 0.079 + 0.016 * s;
- }
- }
- if (x > 0.0) { // Step 5. Calculation of q
- v = t / (s + s); // Step 6.
- if (Math.abs(v) > 0.25) {
- q = q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log(1.0 + v);
- }
- else {
- q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) *
- v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
- } // Step 7. Quotient acceptance
- if (Math.log(1.0 - u) <= q) return(gds/lambda);
- }
-
- for(;;) { // Step 8. Double exponential deviate t
- do {
- e = -Math.log(randomGenerator.raw());
- u = randomGenerator.raw();
- u = u + u - 1.0;
- sign_u = (u > 0)? 1.0 : -1.0;
- t = b + (e * si) * sign_u;
- } while (t <= -0.71874483771719); // Step 9. Rejection of t
- v = t / (s + s); // Step 10. New q(t)
- if (Math.abs(v) > 0.25) {
- q = q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log(1.0 + v);
- }
- else {
- q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) *
- v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
- }
- if (q <= 0.0) continue; // Step 11.
- if (q > 0.5) {
- w = Math.exp(q) - 1.0;
- }
- else {
- w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) *
- q + e1) * q;
- } // Step 12. Hat acceptance
- if ( c * u * sign_u <= w * Math.exp(e - 0.5 * t * t)) {
- x = s + 0.5 * t;
- return(x*x/lambda);
- }
- }
- }
+ double a = alpha;
+ double aa = -1.0, aaa = -1.0,
+ b=0.0, c=0.0, d=0.0, e, r, s=0.0, si=0.0, ss=0.0, q0=0.0,
+ q1 = 0.0416666664, q2 = 0.0208333723, q3 = 0.0079849875,
+ q4 = 0.0015746717, q5 = -0.0003349403, q6 = 0.0003340332,
+ q7 = 0.0006053049, q8 = -0.0004701849, q9 = 0.0001710320,
+ a1 = 0.333333333, a2 = -0.249999949, a3 = 0.199999867,
+ a4 =-0.166677482, a5 = 0.142873973, a6 =-0.124385581,
+ a7 = 0.110368310, a8 = -0.112750886, a9 = 0.104089866,
+ e1 = 1.000000000, e2 = 0.499999994, e3 = 0.166666848,
+ e4 = 0.041664508, e5 = 0.008345522, e6 = 0.001353826,
+ e7 = 0.000247453;
+
+ double gds,p,q,t,sign_u,u,v,w,x;
+ double v1,v2,v12;
+
+ // Check for invalid input values
+
+ if (a <= 0.0) throw new IllegalArgumentException();
+ if (lambda <= 0.0) new IllegalArgumentException();
+
+ if (a < 1.0) { // CASE A: Acceptance rejection algorithm gs
+ b = 1.0 + 0.36788794412 * a; // Step 1
+ for(;;) {
+ p = b * randomGenerator.raw();
+ if (p <= 1.0) { // Step 2. Case gds <= 1
+ gds = Math.exp(Math.log(p) / a);
+ if (Math.log(randomGenerator.raw()) <= -gds) return(gds/lambda);
+ }
+ else { // Step 3. Case gds > 1
+ gds = - Math.log ((b - p) / a);
+ if (Math.log(randomGenerator.raw()) <= ((a - 1.0) * Math.log(gds))) return(gds/lambda);
+ }
+ }
+ }
+
+ else { // CASE B: Acceptance complement algorithm gd (gaussian distribution, box muller transformation)
+ if (a != aa) { // Step 1. Preparations
+ aa = a;
+ ss = a - 0.5;
+ s = Math.sqrt(ss);
+ d = 5.656854249 - 12.0 * s;
+ }
+ // Step 2. Normal deviate
+ do {
+ v1 = 2.0 * randomGenerator.raw() - 1.0;
+ v2 = 2.0 * randomGenerator.raw() - 1.0;
+ v12 = v1*v1 + v2*v2;
+ } while ( v12 > 1.0 );
+ t = v1*Math.sqrt(-2.0*Math.log(v12)/v12);
+ x = s + 0.5 * t;
+ gds = x * x;
+ if (t >= 0.0) return(gds/lambda); // Immediate acceptance
+
+ u = randomGenerator.raw(); // Step 3. Uniform random number
+ if (d * u <= t * t * t) return(gds/lambda); // Squeeze acceptance
+
+ if (a != aaa) { // Step 4. Set-up for hat case
+ aaa = a;
+ r = 1.0 / a;
+ q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) *
+ r + q3) * r + q2) * r + q1) * r;
+ if (a > 3.686) {
+ if (a > 13.022) {
+ b = 1.77;
+ si = 0.75;
+ c = 0.1515 / s;
+ }
+ else {
+ b = 1.654 + 0.0076 * ss;
+ si = 1.68 / s + 0.275;
+ c = 0.062 / s + 0.024;
+ }
+ }
+ else {
+ b = 0.463 + s - 0.178 * ss;
+ si = 1.235;
+ c = 0.195 / s - 0.079 + 0.016 * s;
+ }
+ }
+ if (x > 0.0) { // Step 5. Calculation of q
+ v = t / (s + s); // Step 6.
+ if (Math.abs(v) > 0.25) {
+ q = q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log(1.0 + v);
+ }
+ else {
+ q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) *
+ v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
+ } // Step 7. Quotient acceptance
+ if (Math.log(1.0 - u) <= q) return(gds/lambda);
+ }
+
+ for(;;) { // Step 8. Double exponential deviate t
+ do {
+ e = -Math.log(randomGenerator.raw());
+ u = randomGenerator.raw();
+ u = u + u - 1.0;
+ sign_u = (u > 0)? 1.0 : -1.0;
+ t = b + (e * si) * sign_u;
+ } while (t <= -0.71874483771719); // Step 9. Rejection of t
+ v = t / (s + s); // Step 10. New q(t)
+ if (Math.abs(v) > 0.25) {
+ q = q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log(1.0 + v);
+ }
+ else {
+ q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) *
+ v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
+ }
+ if (q <= 0.0) continue; // Step 11.
+ if (q > 0.5) {
+ w = Math.exp(q) - 1.0;
+ }
+ else {
+ w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) *
+ q + e1) * q;
+ } // Step 12. Hat acceptance
+ if ( c * u * sign_u <= w * Math.exp(e - 0.5 * t * t)) {
+ x = s + 0.5 * t;
+ return(x*x/lambda);
+ }
+ }
+ }
}
/**
* Returns the probability distribution function.
*/
public double pdf(double x) {
- if (x < 0) throw new IllegalArgumentException();
- if (x == 0) {
- if (alpha == 1.0) return 1.0/lambda;
- else return 0.0;
- }
- if (alpha == 1.0) return Math.exp(-x/lambda)/lambda;
+ if (x < 0) throw new IllegalArgumentException();
+ if (x == 0) {
+ if (alpha == 1.0) return 1.0/lambda;
+ else return 0.0;
+ }
+ if (alpha == 1.0) return Math.exp(-x/lambda)/lambda;
- return Math.exp((alpha-1.0) * Math.log(x/lambda) - x/lambda - Fun.logGamma(alpha)) / lambda;
+ return Math.exp((alpha-1.0) * Math.log(x/lambda) - x/lambda - Fun.logGamma(alpha)) / lambda;
}
/**
* Sets the mean and variance.
* @throws IllegalArgumentException if <tt>alpha <= 0.0 || lambda <= 0.0</tt>.
*/
public void setState(double alpha, double lambda) {
- if (alpha <= 0.0) throw new IllegalArgumentException();
- if (lambda <= 0.0) throw new IllegalArgumentException();
- this.alpha = alpha;
- this.lambda = lambda;
+ if (alpha <= 0.0) throw new IllegalArgumentException();
+ if (lambda <= 0.0) throw new IllegalArgumentException();
+ this.alpha = alpha;
+ this.lambda = lambda;
}
/**
* Returns a random number from the distribution.
* @throws IllegalArgumentException if <tt>alpha <= 0.0 || lambda <= 0.0</tt>.
*/
public static double staticNextDouble(double alpha, double lambda) {
- synchronized (shared) {
- return shared.nextDouble(alpha,lambda);
- }
+ synchronized (shared) {
+ return shared.nextDouble(alpha,lambda);
+ }
}
/**
* Returns a String representation of the receiver.
*/
public String toString() {
- return this.getClass().getName()+"("+alpha+","+lambda+")";
+ return this.getClass().getName()+"("+alpha+","+lambda+")";
}
/**
* Sets the uniform random number generated shared by all <b>static</b> methods.
* @param randomGenerator the new uniform random number generator to be shared.
*/
private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
- synchronized (shared) {
- shared.setRandomGenerator(randomGenerator);
- }
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
}
}
Modified: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/HyperGeometric.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/HyperGeometric.java?rev=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/HyperGeometric.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/HyperGeometric.java Wed Nov 25 03:31:47 2009
@@ -32,43 +32,41 @@
* H. Zechner (1994): Efficient sampling from continuous and discrete unimodal distributions,
* Doctoral Dissertation, 156 pp., Technical University Graz, Austria.
*
- * @author wolfgang.hoschek@cern.ch
- * @version 1.0, 09/24/99
*/
/**
* @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
*/
@Deprecated
public class HyperGeometric extends AbstractDiscreteDistribution {
- protected int my_N;
- protected int my_s;
- protected int my_n;
-
- // cached vars shared by hmdu(...) and hprs(...)
- private int N_last = -1, M_last = -1, n_last = -1;
- private int N_Mn, m;
-
- // cached vars for hmdu(...)
- private int mp, b;
- private double Mp, np, fm;
-
- // cached vars for hprs(...)
- private int k2, k4, k1, k5;
- private double dl, dr, r1, r2, r4, r5, ll, lr, c_pm,
- f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
-
-
- // The uniform random number generated shared by all <b>static</b> methods.
- protected static HyperGeometric shared = new HyperGeometric(1,1,1,makeDefaultGenerator());
+ protected int my_N;
+ protected int my_s;
+ protected int my_n;
+
+ // cached vars shared by hmdu(...) and hprs(...)
+ private int N_last = -1, M_last = -1, n_last = -1;
+ private int N_Mn, m;
+
+ // cached vars for hmdu(...)
+ private int mp, b;
+ private double Mp, np, fm;
+
+ // cached vars for hprs(...)
+ private int k2, k4, k1, k5;
+ private double dl, dr, r1, r2, r4, r5, ll, lr, c_pm,
+ f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
+
+
+ // The uniform random number generated shared by all <b>static</b> methods.
+ protected static HyperGeometric shared = new HyperGeometric(1,1,1,makeDefaultGenerator());
/**
* Constructs a HyperGeometric distribution.
*/
public HyperGeometric(int N, int s, int n, RandomEngine randomGenerator) {
- setRandomGenerator(randomGenerator);
- setState(N,s,n);
+ setRandomGenerator(randomGenerator);
+ setState(N,s,n);
}
private static double fc_lnpk(int k, int N_Mn, int M, int n) {
- return(Arithmetic.logFactorial(k) + Arithmetic.logFactorial(M - k) + Arithmetic.logFactorial(n - k) + Arithmetic.logFactorial(N_Mn + k));
+ return(Arithmetic.logFactorial(k) + Arithmetic.logFactorial(M - k) + Arithmetic.logFactorial(n - k) + Arithmetic.logFactorial(N_Mn + k));
}
/**
* Returns a random number from the distribution.
@@ -79,209 +77,209 @@
double p, nu, c, d, U;
if (N != N_last || M != M_last || n != n_last) { // set-up */
- N_last = N;
- M_last = M;
- n_last = n;
-
- Mp = (double) (M + 1);
- np = (double) (n + 1); N_Mn = N - M - n;
-
- p = Mp / (N + 2.0);
- nu = np * p; /* mode, real */
- if ((m = (int) nu) == nu && p == 0.5) { /* mode, integer */
- mp = m--;
- }
- else {
- mp = m + 1; /* mp = m + 1 */
- }
+ N_last = N;
+ M_last = M;
+ n_last = n;
+
+ Mp = (double) (M + 1);
+ np = (double) (n + 1); N_Mn = N - M - n;
+
+ p = Mp / (N + 2.0);
+ nu = np * p; /* mode, real */
+ if ((m = (int) nu) == nu && p == 0.5) { /* mode, integer */
+ mp = m--;
+ }
+ else {
+ mp = m + 1; /* mp = m + 1 */
+ }
/* mode probability, using the external function flogfak(k) = ln(k!) */
- fm = Math.exp(Arithmetic.logFactorial(N - M) - Arithmetic.logFactorial(N_Mn + m) - Arithmetic.logFactorial(n - m)
- + Arithmetic.logFactorial(M) - Arithmetic.logFactorial(M - m) - Arithmetic.logFactorial(m)
- - Arithmetic.logFactorial(N) + Arithmetic.logFactorial(N - n) + Arithmetic.logFactorial(n) );
+ fm = Math.exp(Arithmetic.logFactorial(N - M) - Arithmetic.logFactorial(N_Mn + m) - Arithmetic.logFactorial(n - m)
+ + Arithmetic.logFactorial(M) - Arithmetic.logFactorial(M - m) - Arithmetic.logFactorial(m)
+ - Arithmetic.logFactorial(N) + Arithmetic.logFactorial(N - n) + Arithmetic.logFactorial(n) );
/* safety bound - guarantees at least 17 significant decimal digits */
/* b = min(n, (long int)(nu + k*c')) */
- b = (int) (nu + 11.0 * Math.sqrt(nu * (1.0 - p) * (1.0 - n/(double)N) + 1.0));
- if (b > n) b = n;
- }
-
- for (;;) {
- if ((U = randomGenerator.raw() - fm) <= 0.0) return(m);
- c = d = fm;
+ b = (int) (nu + 11.0 * Math.sqrt(nu * (1.0 - p) * (1.0 - n/(double)N) + 1.0));
+ if (b > n) b = n;
+ }
+
+ for (;;) {
+ if ((U = randomGenerator.raw() - fm) <= 0.0) return(m);
+ c = d = fm;
/* down- and upward search from the mode */
- for (I = 1; I <= m; I++) {
- K = mp - I; /* downward search */
- c *= (double)K/(np - K) * ((double)(N_Mn + K)/(Mp - K));
- if ((U -= c) <= 0.0) return(K - 1);
-
- K = m + I; /* upward search */
- d *= (np - K)/(double)K * ((Mp - K)/(double)(N_Mn + K));
- if ((U -= d) <= 0.0) return(K);
- }
+ for (I = 1; I <= m; I++) {
+ K = mp - I; /* downward search */
+ c *= (double)K/(np - K) * ((double)(N_Mn + K)/(Mp - K));
+ if ((U -= c) <= 0.0) return(K - 1);
+
+ K = m + I; /* upward search */
+ d *= (np - K)/(double)K * ((Mp - K)/(double)(N_Mn + K));
+ if ((U -= d) <= 0.0) return(K);
+ }
/* upward search from K = 2m + 1 to K = b */
- for (K = mp + m; K <= b; K++) {
- d *= (np - K)/(double)K * ((Mp - K)/(double)(N_Mn + K));
- if ((U -= d) <= 0.0) return(K);
- }
- }
+ for (K = mp + m; K <= b; K++) {
+ d *= (np - K)/(double)K * ((Mp - K)/(double)(N_Mn + K));
+ if ((U -= d) <= 0.0) return(K);
+ }
+ }
}
/**
* Returns a random number from the distribution.
*/
protected int hprs(int N, int M, int n, RandomEngine randomGenerator) {
- int Dk, X, V;
- double Mp, np, p, nu, U, Y, W; /* (X, Y) <-> (V, W) */
+ int Dk, X, V;
+ double Mp, np, p, nu, U, Y, W; /* (X, Y) <-> (V, W) */
- if (N != N_last || M != M_last || n != n_last) { /* set-up */
- N_last = N;
- M_last = M;
- n_last = n;
-
- Mp = (double) (M + 1);
- np = (double) (n + 1); N_Mn = N - M - n;
-
- p = Mp / (N + 2.0); nu = np * p; // main parameters
-
- // approximate deviation of reflection points k2, k4 from nu - 1/2
- U = Math.sqrt(nu * (1.0 - p) * (1.0 - (n + 2.0)/(N + 3.0)) + 0.25);
-
- // mode m, reflection points k2 and k4, and points k1 and k5, which
- // delimit the centre region of h(x)
- // k2 = ceil (nu - 1/2 - U), k1 = 2*k2 - (m - 1 + delta_ml)
- // k4 = floor(nu - 1/2 + U), k5 = 2*k4 - (m + 1 - delta_mr)
+ if (N != N_last || M != M_last || n != n_last) { /* set-up */
+ N_last = N;
+ M_last = M;
+ n_last = n;
+
+ Mp = (double) (M + 1);
+ np = (double) (n + 1); N_Mn = N - M - n;
+
+ p = Mp / (N + 2.0); nu = np * p; // main parameters
+
+ // approximate deviation of reflection points k2, k4 from nu - 1/2
+ U = Math.sqrt(nu * (1.0 - p) * (1.0 - (n + 2.0)/(N + 3.0)) + 0.25);
+
+ // mode m, reflection points k2 and k4, and points k1 and k5, which
+ // delimit the centre region of h(x)
+ // k2 = ceil (nu - 1/2 - U), k1 = 2*k2 - (m - 1 + delta_ml)
+ // k4 = floor(nu - 1/2 + U), k5 = 2*k4 - (m + 1 - delta_mr)
- m = (int) nu;
- k2 = (int) Math.ceil(nu - 0.5 - U); if (k2 >= m) k2 = m - 1;
- k4 = (int) (nu - 0.5 + U);
- k1 = k2 + k2 - m + 1; // delta_ml = 0
- k5 = k4 + k4 - m; // delta_mr = 1
-
- // range width of the critical left and right centre region
- dl = (double) (k2 - k1);
- dr = (double) (k5 - k4);
-
- // recurrence constants r(k) = p(k)/p(k-1) at k = k1, k2, k4+1, k5+1
- r1 = (np/(double) k1 - 1.0) * (Mp - k1)/(double)(N_Mn + k1);
- r2 = (np/(double) k2 - 1.0) * (Mp - k2)/(double)(N_Mn + k2);
- r4 = (np/(double)(k4 + 1) - 1.0) * (M - k4)/(double)(N_Mn + k4 + 1);
- r5 = (np/(double)(k5 + 1) - 1.0) * (M - k5)/(double)(N_Mn + k5 + 1);
-
- // reciprocal values of the scale parameters of expon. tail envelopes
- ll = Math.log(r1); // expon. tail left //
- lr = -Math.log(r5); // expon. tail right //
-
- // hypergeom. constant, necessary for computing function values f(k)
- c_pm = fc_lnpk(m, N_Mn, M, n);
-
- // function values f(k) = p(k)/p(m) at k = k2, k4, k1, k5
- f2 = Math.exp(c_pm - fc_lnpk(k2, N_Mn, M, n));
- f4 = Math.exp(c_pm - fc_lnpk(k4, N_Mn, M, n));
- f1 = Math.exp(c_pm - fc_lnpk(k1, N_Mn, M, n));
- f5 = Math.exp(c_pm - fc_lnpk(k5, N_Mn, M, n));
-
- // area of the two centre and the two exponential tail regions
- // area of the two immediate acceptance regions between k2, k4
- p1 = f2 * (dl + 1.0); // immed. left
- p2 = f2 * dl + p1; // centre left
- p3 = f4 * (dr + 1.0) + p2; // immed. right
- p4 = f4 * dr + p3; // centre right
- p5 = f1 / ll + p4; // expon. tail left
- p6 = f5 / lr + p5; // expon. tail right
- }
-
- for (;;) {
- // generate uniform number U -- U(0, p6)
- // case distinction corresponding to U
- if ((U = randomGenerator.raw() * p6) < p2) { // centre left
-
- // immediate acceptance region R2 = [k2, m) *[0, f2), X = k2, ... m -1
- if ((W = U - p1) < 0.0) return(k2 + (int)(U/f2));
- // immediate acceptance region R1 = [k1, k2)*[0, f1), X = k1, ... k2-1
- if ((Y = W / dl) < f1 ) return(k1 + (int)(W/f1));
-
- // computation of candidate X < k2, and its counterpart V > k2
- // either squeeze-acceptance of X or acceptance-rejection of V
- Dk = (int)(dl * randomGenerator.raw()) + 1;
- if (Y <= f2 - Dk * (f2 - f2/r2)) { // quick accept of
- return(k2 - Dk); // X = k2 - Dk
- }
- if ((W = f2 + f2 - Y) < 1.0) { // quick reject of V
- V = k2 + Dk;
- if (W <= f2 + Dk * (1.0 - f2)/(dl + 1.0)) { // quick accept of
- return(V); // V = k2 + Dk
- }
- if (Math.log(W) <= c_pm - fc_lnpk(V, N_Mn, M, n)) {
- return(V); // final accept of V
- }
- }
- X = k2 - Dk;
- }
- else if (U < p4) { // centre right
-
- // immediate acceptance region R3 = [m, k4+1)*[0, f4), X = m, ... k4
- if ((W = U - p3) < 0.0) return(k4 - (int)((U - p2)/f4));
- // immediate acceptance region R4 = [k4+1, k5+1)*[0, f5)
- if ((Y = W / dr) < f5 ) return(k5 - (int)(W/f5));
-
- // computation of candidate X > k4, and its counterpart V < k4
- // either squeeze-acceptance of X or acceptance-rejection of V
- Dk = (int)(dr * randomGenerator.raw()) + 1;
- if (Y <= f4 - Dk * (f4 - f4*r4)) { // quick accept of
- return(k4 + Dk); // X = k4 + Dk
- }
- if ((W = f4 + f4 - Y) < 1.0) { // quick reject of V
- V = k4 - Dk;
- if (W <= f4 + Dk * (1.0 - f4)/dr) { // quick accept of
- return(V); // V = k4 - Dk
- }
- if (Math.log(W) <= c_pm - fc_lnpk(V, N_Mn, M, n)) {
- return(V); // final accept of V
- }
- }
- X = k4 + Dk;
- }
- else {
- Y = randomGenerator.raw();
- if (U < p5) { // expon. tail left
- Dk = (int)(1.0 - Math.log(Y)/ll);
- if ((X = k1 - Dk) < 0) continue; // 0 <= X <= k1 - 1
- Y *= (U - p4) * ll; // Y -- U(0, h(x))
- if (Y <= f1 - Dk * (f1 - f1/r1)) {
- return(X); // quick accept of X
- }
- }
- else { // expon. tail right
- Dk = (int)(1.0 - Math.log(Y)/lr);
- if ((X = k5 + Dk) > n ) continue; // k5 + 1 <= X <= n
- Y *= (U - p5) * lr; // Y -- U(0, h(x)) /
- if (Y <= f5 - Dk * (f5 - f5*r5)) {
- return(X); // quick accept of X
- }
- }
- }
-
- // acceptance-rejection test of candidate X from the original area
- // test, whether Y <= f(X), with Y = U*h(x) and U -- U(0, 1)
- // log f(X) = log( m! (M - m)! (n - m)! (N - M - n + m)! )
- // - log( X! (M - X)! (n - X)! (N - M - n + X)! )
- // by using an external function for log k!
- if (Math.log(Y) <= c_pm - fc_lnpk(X, N_Mn, M, n)) return(X);
- }
+ m = (int) nu;
+ k2 = (int) Math.ceil(nu - 0.5 - U); if (k2 >= m) k2 = m - 1;
+ k4 = (int) (nu - 0.5 + U);
+ k1 = k2 + k2 - m + 1; // delta_ml = 0
+ k5 = k4 + k4 - m; // delta_mr = 1
+
+ // range width of the critical left and right centre region
+ dl = (double) (k2 - k1);
+ dr = (double) (k5 - k4);
+
+ // recurrence constants r(k) = p(k)/p(k-1) at k = k1, k2, k4+1, k5+1
+ r1 = (np/(double) k1 - 1.0) * (Mp - k1)/(double)(N_Mn + k1);
+ r2 = (np/(double) k2 - 1.0) * (Mp - k2)/(double)(N_Mn + k2);
+ r4 = (np/(double)(k4 + 1) - 1.0) * (M - k4)/(double)(N_Mn + k4 + 1);
+ r5 = (np/(double)(k5 + 1) - 1.0) * (M - k5)/(double)(N_Mn + k5 + 1);
+
+ // reciprocal values of the scale parameters of expon. tail envelopes
+ ll = Math.log(r1); // expon. tail left //
+ lr = -Math.log(r5); // expon. tail right //
+
+ // hypergeom. constant, necessary for computing function values f(k)
+ c_pm = fc_lnpk(m, N_Mn, M, n);
+
+ // function values f(k) = p(k)/p(m) at k = k2, k4, k1, k5
+ f2 = Math.exp(c_pm - fc_lnpk(k2, N_Mn, M, n));
+ f4 = Math.exp(c_pm - fc_lnpk(k4, N_Mn, M, n));
+ f1 = Math.exp(c_pm - fc_lnpk(k1, N_Mn, M, n));
+ f5 = Math.exp(c_pm - fc_lnpk(k5, N_Mn, M, n));
+
+ // area of the two centre and the two exponential tail regions
+ // area of the two immediate acceptance regions between k2, k4
+ p1 = f2 * (dl + 1.0); // immed. left
+ p2 = f2 * dl + p1; // centre left
+ p3 = f4 * (dr + 1.0) + p2; // immed. right
+ p4 = f4 * dr + p3; // centre right
+ p5 = f1 / ll + p4; // expon. tail left
+ p6 = f5 / lr + p5; // expon. tail right
+ }
+
+ for (;;) {
+ // generate uniform number U -- U(0, p6)
+ // case distinction corresponding to U
+ if ((U = randomGenerator.raw() * p6) < p2) { // centre left
+
+ // immediate acceptance region R2 = [k2, m) *[0, f2), X = k2, ... m -1
+ if ((W = U - p1) < 0.0) return(k2 + (int)(U/f2));
+ // immediate acceptance region R1 = [k1, k2)*[0, f1), X = k1, ... k2-1
+ if ((Y = W / dl) < f1 ) return(k1 + (int)(W/f1));
+
+ // computation of candidate X < k2, and its counterpart V > k2
+ // either squeeze-acceptance of X or acceptance-rejection of V
+ Dk = (int)(dl * randomGenerator.raw()) + 1;
+ if (Y <= f2 - Dk * (f2 - f2/r2)) { // quick accept of
+ return(k2 - Dk); // X = k2 - Dk
+ }
+ if ((W = f2 + f2 - Y) < 1.0) { // quick reject of V
+ V = k2 + Dk;
+ if (W <= f2 + Dk * (1.0 - f2)/(dl + 1.0)) { // quick accept of
+ return(V); // V = k2 + Dk
+ }
+ if (Math.log(W) <= c_pm - fc_lnpk(V, N_Mn, M, n)) {
+ return(V); // final accept of V
+ }
+ }
+ X = k2 - Dk;
+ }
+ else if (U < p4) { // centre right
+
+ // immediate acceptance region R3 = [m, k4+1)*[0, f4), X = m, ... k4
+ if ((W = U - p3) < 0.0) return(k4 - (int)((U - p2)/f4));
+ // immediate acceptance region R4 = [k4+1, k5+1)*[0, f5)
+ if ((Y = W / dr) < f5 ) return(k5 - (int)(W/f5));
+
+ // computation of candidate X > k4, and its counterpart V < k4
+ // either squeeze-acceptance of X or acceptance-rejection of V
+ Dk = (int)(dr * randomGenerator.raw()) + 1;
+ if (Y <= f4 - Dk * (f4 - f4*r4)) { // quick accept of
+ return(k4 + Dk); // X = k4 + Dk
+ }
+ if ((W = f4 + f4 - Y) < 1.0) { // quick reject of V
+ V = k4 - Dk;
+ if (W <= f4 + Dk * (1.0 - f4)/dr) { // quick accept of
+ return(V); // V = k4 - Dk
+ }
+ if (Math.log(W) <= c_pm - fc_lnpk(V, N_Mn, M, n)) {
+ return(V); // final accept of V
+ }
+ }
+ X = k4 + Dk;
+ }
+ else {
+ Y = randomGenerator.raw();
+ if (U < p5) { // expon. tail left
+ Dk = (int)(1.0 - Math.log(Y)/ll);
+ if ((X = k1 - Dk) < 0) continue; // 0 <= X <= k1 - 1
+ Y *= (U - p4) * ll; // Y -- U(0, h(x))
+ if (Y <= f1 - Dk * (f1 - f1/r1)) {
+ return(X); // quick accept of X
+ }
+ }
+ else { // expon. tail right
+ Dk = (int)(1.0 - Math.log(Y)/lr);
+ if ((X = k5 + Dk) > n ) continue; // k5 + 1 <= X <= n
+ Y *= (U - p5) * lr; // Y -- U(0, h(x)) /
+ if (Y <= f5 - Dk * (f5 - f5*r5)) {
+ return(X); // quick accept of X
+ }
+ }
+ }
+
+ // acceptance-rejection test of candidate X from the original area
+ // test, whether Y <= f(X), with Y = U*h(x) and U -- U(0, 1)
+ // log f(X) = log( m! (M - m)! (n - m)! (N - M - n + m)! )
+ // - log( X! (M - X)! (n - X)! (N - M - n + X)! )
+ // by using an external function for log k!
+ if (Math.log(Y) <= c_pm - fc_lnpk(X, N_Mn, M, n)) return(X);
+ }
}
/**
* Returns a random number from the distribution.
*/
public int nextInt() {
- return nextInt(this.my_N, this.my_s, this.my_n, this.randomGenerator);
+ return nextInt(this.my_N, this.my_s, this.my_n, this.randomGenerator);
}
/**
* Returns a random number from the distribution; bypasses the internal state.
*/
public int nextInt(int N, int s, int n) {
- return nextInt(N,s,n,this.randomGenerator);
+ return nextInt(N,s,n,this.randomGenerator);
}
/**
* Returns a random number from the distribution; bypasses the internal state.
@@ -329,66 +327,66 @@
* long integer N , M , n. *
* *
******************************************************************/
- int Nhalf, n_le_Nhalf, M_le_Nhalf, K;
+ int Nhalf, n_le_Nhalf, M_le_Nhalf, K;
- Nhalf = N / 2;
- n_le_Nhalf = (n <= Nhalf) ? n : N - n;
- M_le_Nhalf = (M <= Nhalf) ? M : N - M;
-
- if ((n*M/N) < 10) {
- K = (n_le_Nhalf <= M_le_Nhalf)
- ? hmdu(N, M_le_Nhalf, n_le_Nhalf, randomGenerator)
- : hmdu(N, n_le_Nhalf, M_le_Nhalf, randomGenerator);
- }
- else {
- K = (n_le_Nhalf <= M_le_Nhalf)
- ? hprs(N, M_le_Nhalf, n_le_Nhalf, randomGenerator)
- : hprs(N, n_le_Nhalf, M_le_Nhalf, randomGenerator);
- }
-
- if (n <= Nhalf) {
- return (M <= Nhalf) ? K : n - K;
- }
- else {
- return (M <= Nhalf) ? M - K : n - N + M + K;
- }
+ Nhalf = N / 2;
+ n_le_Nhalf = (n <= Nhalf) ? n : N - n;
+ M_le_Nhalf = (M <= Nhalf) ? M : N - M;
+
+ if ((n*M/N) < 10) {
+ K = (n_le_Nhalf <= M_le_Nhalf)
+ ? hmdu(N, M_le_Nhalf, n_le_Nhalf, randomGenerator)
+ : hmdu(N, n_le_Nhalf, M_le_Nhalf, randomGenerator);
+ }
+ else {
+ K = (n_le_Nhalf <= M_le_Nhalf)
+ ? hprs(N, M_le_Nhalf, n_le_Nhalf, randomGenerator)
+ : hprs(N, n_le_Nhalf, M_le_Nhalf, randomGenerator);
+ }
+
+ if (n <= Nhalf) {
+ return (M <= Nhalf) ? K : n - K;
+ }
+ else {
+ return (M <= Nhalf) ? M - K : n - N + M + K;
+ }
}
/**
* Returns the probability distribution function.
*/
public double pdf(int k) {
- return Arithmetic.binomial(my_s, k) * Arithmetic.binomial(my_N - my_s, my_n - k)
- / Arithmetic.binomial(my_N, my_n);
+ return Arithmetic.binomial(my_s, k) * Arithmetic.binomial(my_N - my_s, my_n - k)
+ / Arithmetic.binomial(my_N, my_n);
}
/**
* Sets the parameters.
*/
public void setState(int N, int s, int n) {
- this.my_N = N;
- this.my_s = s;
- this.my_n = n;
+ this.my_N = N;
+ this.my_s = s;
+ this.my_n = n;
}
/**
* Returns a random number from the distribution.
*/
public static double staticNextInt(int N, int M, int n) {
- synchronized (shared) {
- return shared.nextInt(N,M,n);
- }
+ synchronized (shared) {
+ return shared.nextInt(N,M,n);
+ }
}
/**
* Returns a String representation of the receiver.
*/
public String toString() {
- return this.getClass().getName()+"("+my_N+","+my_s+","+my_n+")";
+ return this.getClass().getName()+"("+my_N+","+my_s+","+my_n+")";
}
/**
* Sets the uniform random number generated shared by all <b>static</b> methods.
* @param randomGenerator the new uniform random number generator to be shared.
*/
private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
- synchronized (shared) {
- shared.setRandomGenerator(randomGenerator);
- }
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
}
}
Modified: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Hyperbolic.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Hyperbolic.java?rev=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Hyperbolic.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Hyperbolic.java Wed Nov 25 03:31:47 2009
@@ -26,37 +26,35 @@
* <p>
* L. Devroye (1986): Non-Uniform Random Variate Generation, Springer Verlag, New York.
*
- * @author wolfgang.hoschek@cern.ch
- * @version 1.0, 09/24/99
*/
/**
* @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
*/
@Deprecated
public class Hyperbolic extends AbstractContinousDistribution {
- protected double alpha;
- protected double beta;
+ protected double alpha;
+ protected double beta;
- // cached values shared for generateHyperbolic(...)
- protected double a_setup = 0.0, b_setup = -1.0;
- protected double x, u, v, e;
- protected double hr, hl, s, pm, pr, samb, pmr, mpa_1, mmb_1;
+ // cached values shared for generateHyperbolic(...)
+ protected double a_setup = 0.0, b_setup = -1.0;
+ protected double x, u, v, e;
+ protected double hr, hl, s, pm, pr, samb, pmr, mpa_1, mmb_1;
- // The uniform random number generated shared by all <b>static</b> methods.
- protected static Hyperbolic shared = new Hyperbolic(10.0,10.0,makeDefaultGenerator());
+ // The uniform random number generated shared by all <b>static</b> methods.
+ protected static Hyperbolic shared = new Hyperbolic(10.0,10.0,makeDefaultGenerator());
/**
* Constructs a Beta distribution.
*/
public Hyperbolic(double alpha, double beta, RandomEngine randomGenerator) {
- setRandomGenerator(randomGenerator);
- setState(alpha,beta);
+ setRandomGenerator(randomGenerator);
+ setState(alpha,beta);
}
/**
* Returns a random number from the distribution.
*/
public double nextDouble() {
- return nextDouble(alpha, beta);
+ return nextDouble(alpha, beta);
}
/**
* Returns a hyperbolic distributed random number; bypasses the internal state.
@@ -79,94 +77,94 @@
* unsigned long integer *seed. *
* *
******************************************************************/
- double a = alpha;
- double b = beta;
+ double a = alpha;
+ double b = beta;
- if ((a_setup != a) || (b_setup != b)) { // SET-UP
- double mpa, mmb, mode;
- double amb;
- double a_, b_, a_1, b_1, pl;
- double help_1, help_2;
- amb = a*a - b*b; // a^2 - b^2
- samb = Math.sqrt(amb); // -log(f(mode))
- mode = b/samb; // mode
- help_1 = a*Math.sqrt(2.0*samb + 1.0);
- help_2 = b*(samb + 1.0);
- mpa = (help_2 + help_1)/amb; // fr^-1(exp(-sqrt(a^2 - b^2) - 1.0))
- mmb = (help_2 - help_1)/amb; // fl^-1(exp(-sqrt(a^2 - b^2) - 1.0))
- a_ = mpa - mode;
- b_ = -mmb + mode;
- hr = -1.0/(-a*mpa/Math.sqrt(1.0 + mpa*mpa) + b);
- hl = 1.0/(-a*mmb/Math.sqrt(1.0 + mmb*mmb) + b);
- a_1 = a_ - hr;
- b_1 = b_ - hl;
- mmb_1 = mode - b_1; // lower border
- mpa_1 = mode + a_1; // upper border
-
- s = (a_ + b_);
- pm = (a_1 + b_1)/s;
- pr = hr/s;
- pmr = pm + pr;
-
- a_setup = a;
- b_setup = b;
- }
-
- // GENERATOR
- for(;;) {
- u = randomGenerator.raw();
- v = randomGenerator.raw();
- if (u <= pm) { // Rejection with a uniform majorizing function
- // over the body of the distribution
- x = mmb_1 + u*s;
- if (Math.log(v) <= (-a*Math.sqrt(1.0 + x*x) + b*x + samb)) break;
- }
- else {
- if (u <= pmr) { // Rejection with an exponential envelope on the
- // right side of the mode
- e = -Math.log((u-pm)/pr);
- x = mpa_1 + hr*e;
- if ((Math.log(v) - e) <= (-a*Math.sqrt(1.0 + x*x) + b*x + samb)) break;
- }
- else { // Rejection with an exponential envelope on the
- // left side of the mode
- e = Math.log((u-pmr)/(1.0 - pmr));
- x = mmb_1 + hl*e;
- if ((Math.log(v) + e) <= (-a*Math.sqrt(1.0 + x*x) + b*x + samb)) break;
- }
- }
- }
+ if ((a_setup != a) || (b_setup != b)) { // SET-UP
+ double mpa, mmb, mode;
+ double amb;
+ double a_, b_, a_1, b_1, pl;
+ double help_1, help_2;
+ amb = a*a - b*b; // a^2 - b^2
+ samb = Math.sqrt(amb); // -log(f(mode))
+ mode = b/samb; // mode
+ help_1 = a*Math.sqrt(2.0*samb + 1.0);
+ help_2 = b*(samb + 1.0);
+ mpa = (help_2 + help_1)/amb; // fr^-1(exp(-sqrt(a^2 - b^2) - 1.0))
+ mmb = (help_2 - help_1)/amb; // fl^-1(exp(-sqrt(a^2 - b^2) - 1.0))
+ a_ = mpa - mode;
+ b_ = -mmb + mode;
+ hr = -1.0/(-a*mpa/Math.sqrt(1.0 + mpa*mpa) + b);
+ hl = 1.0/(-a*mmb/Math.sqrt(1.0 + mmb*mmb) + b);
+ a_1 = a_ - hr;
+ b_1 = b_ - hl;
+ mmb_1 = mode - b_1; // lower border
+ mpa_1 = mode + a_1; // upper border
+
+ s = (a_ + b_);
+ pm = (a_1 + b_1)/s;
+ pr = hr/s;
+ pmr = pm + pr;
+
+ a_setup = a;
+ b_setup = b;
+ }
+
+ // GENERATOR
+ for(;;) {
+ u = randomGenerator.raw();
+ v = randomGenerator.raw();
+ if (u <= pm) { // Rejection with a uniform majorizing function
+ // over the body of the distribution
+ x = mmb_1 + u*s;
+ if (Math.log(v) <= (-a*Math.sqrt(1.0 + x*x) + b*x + samb)) break;
+ }
+ else {
+ if (u <= pmr) { // Rejection with an exponential envelope on the
+ // right side of the mode
+ e = -Math.log((u-pm)/pr);
+ x = mpa_1 + hr*e;
+ if ((Math.log(v) - e) <= (-a*Math.sqrt(1.0 + x*x) + b*x + samb)) break;
+ }
+ else { // Rejection with an exponential envelope on the
+ // left side of the mode
+ e = Math.log((u-pmr)/(1.0 - pmr));
+ x = mmb_1 + hl*e;
+ if ((Math.log(v) + e) <= (-a*Math.sqrt(1.0 + x*x) + b*x + samb)) break;
+ }
+ }
+ }
- return(x);
+ return(x);
}
/**
* Sets the parameters.
*/
public void setState(double alpha, double beta) {
- this.alpha = alpha;
- this.beta = beta;
+ this.alpha = alpha;
+ this.beta = beta;
}
/**
* Returns a random number from the distribution.
*/
public static double staticNextDouble(double alpha, double beta) {
- synchronized (shared) {
- return shared.nextDouble(alpha,beta);
- }
+ synchronized (shared) {
+ return shared.nextDouble(alpha,beta);
+ }
}
/**
* Returns a String representation of the receiver.
*/
public String toString() {
- return this.getClass().getName()+"("+alpha+","+beta+")";
+ return this.getClass().getName()+"("+alpha+","+beta+")";
}
/**
* Sets the uniform random number generated shared by all <b>static</b> methods.
* @param randomGenerator the new uniform random number generator to be shared.
*/
private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
- synchronized (shared) {
- shared.setRandomGenerator(randomGenerator);
- }
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
}
}
Modified: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Logarithmic.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Logarithmic.java?rev=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Logarithmic.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Logarithmic.java Wed Nov 25 03:31:47 2009
@@ -27,33 +27,31 @@
* <p>
* A.W. Kemp (1981): Efficient generation of logarithmically distributed pseudo-random variables, Appl. Statist. 30, 249-253.
*
- * @author wolfgang.hoschek@cern.ch
- * @version 1.0, 09/24/99
*/
/**
* @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
*/
@Deprecated
public class Logarithmic extends AbstractContinousDistribution {
- protected double my_p;
+ protected double my_p;
- // cached vars for method nextDouble(a) (for performance only)
- private double t,h,a_prev = -1.0;
+ // cached vars for method nextDouble(a) (for performance only)
+ private double t,h,a_prev = -1.0;
- // The uniform random number generated shared by all <b>static</b> methods.
- protected static Logarithmic shared = new Logarithmic(0.5,makeDefaultGenerator());
+ // The uniform random number generated shared by all <b>static</b> methods.
+ protected static Logarithmic shared = new Logarithmic(0.5,makeDefaultGenerator());
/**
* Constructs a Logarithmic distribution.
*/
public Logarithmic(double p, RandomEngine randomGenerator) {
- setRandomGenerator(randomGenerator);
- setState(p);
+ setRandomGenerator(randomGenerator);
+ setState(p);
}
/**
* Returns a random number from the distribution.
*/
public double nextDouble() {
- return nextDouble(this.my_p);
+ return nextDouble(this.my_p);
}
/**
* Returns a random number from the distribution; bypasses the internal state.
@@ -92,66 +90,66 @@
* unsigned long integer *seed. *
* *
******************************************************************/
- double u,v,p,q;
- int k;
+ double u,v,p,q;
+ int k;
- if (a != a_prev) { // Set-up
- a_prev = a;
- if (a<0.97) t = -a / Math.log(1.0 - a);
- else h=Math.log(1.0 - a);
- }
+ if (a != a_prev) { // Set-up
+ a_prev = a;
+ if (a<0.97) t = -a / Math.log(1.0 - a);
+ else h=Math.log(1.0 - a);
+ }
- u=randomGenerator.raw();
- if (a<0.97) { // Inversion/Chop-down
- k = 1;
- p = t;
- while (u > p) {
- //System.out.println("u="+u+", p="+p);
- u -= p;
- k++;
- p *= a * (k-1.0)/(double)k;
- }
- return k;
- }
+ u=randomGenerator.raw();
+ if (a<0.97) { // Inversion/Chop-down
+ k = 1;
+ p = t;
+ while (u > p) {
+ //System.out.println("u="+u+", p="+p);
+ u -= p;
+ k++;
+ p *= a * (k-1.0)/(double)k;
+ }
+ return k;
+ }
- if (u > a) return 1; // Transformation
- u=randomGenerator.raw();
- v = u;
- q = 1.0 - Math.exp(v * h);
- if ( u <= q * q) {
- k = (int) (1 + Math.log(u) / Math.log(q));
- return k;
- }
- if (u > q) return 1;
- return 2;
+ if (u > a) return 1; // Transformation
+ u=randomGenerator.raw();
+ v = u;
+ q = 1.0 - Math.exp(v * h);
+ if ( u <= q * q) {
+ k = (int) (1 + Math.log(u) / Math.log(q));
+ return k;
+ }
+ if (u > q) return 1;
+ return 2;
}
/**
* Sets the distribution parameter.
*/
public void setState(double p) {
- this.my_p = p;
+ this.my_p = p;
}
/**
* Returns a random number from the distribution.
*/
public static double staticNextDouble(double p) {
- synchronized (shared) {
- return shared.nextDouble(p);
- }
+ synchronized (shared) {
+ return shared.nextDouble(p);
+ }
}
/**
* Returns a String representation of the receiver.
*/
public String toString() {
- return this.getClass().getName()+"("+my_p+")";
+ return this.getClass().getName()+"("+my_p+")";
}
/**
* Sets the uniform random number generated shared by all <b>static</b> methods.
* @param randomGenerator the new uniform random number generator to be shared.
*/
private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
- synchronized (shared) {
- shared.setRandomGenerator(randomGenerator);
- }
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
}
}