Return-Path: 
 <commits-return-23424-apmail-commons-commits-archive=commons.apache.org@commons.apache.org>
X-Original-To: apmail-commons-commits-archive@minotaur.apache.org
Delivered-To: apmail-commons-commits-archive@minotaur.apache.org
Received: from mail.apache.org (hermes.apache.org [140.211.11.3])
	by minotaur.apache.org (Postfix) with SMTP id F229C9113
	for <apmail-commons-commits-archive@minotaur.apache.org>;
 Sun,  6 Nov 2011 07:18:23 +0000 (UTC)
Received: (qmail 10384 invoked by uid 500); 6 Nov 2011 07:18:23 -0000
Delivered-To: apmail-commons-commits-archive@commons.apache.org
Received: (qmail 10331 invoked by uid 500); 6 Nov 2011 07:18:19 -0000
Mailing-List: contact commits-help@commons.apache.org; run by ezmlm
Precedence: bulk
List-Help: <mailto:commits-help@commons.apache.org>
List-Unsubscribe: <mailto:commits-unsubscribe@commons.apache.org>
List-Post: <mailto:commits@commons.apache.org>
List-Id: <commits.commons.apache.org>
Reply-To: dev@commons.apache.org
Delivered-To: mailing list commits@commons.apache.org
Received: (qmail 10323 invoked by uid 99); 6 Nov 2011 07:18:17 -0000
Received: from nike.apache.org (HELO nike.apache.org) (192.87.106.230)
    by apache.org (qpsmtpd/0.29) with ESMTP; Sun, 06 Nov 2011 07:18:17 +0000
X-ASF-Spam-Status: No, hits=-2000.0 required=5.0
	tests=ALL_TRUSTED
X-Spam-Check-By: apache.org
Received: from [140.211.11.4] (HELO eris.apache.org) (140.211.11.4)
    by apache.org (qpsmtpd/0.29) with ESMTP; Sun, 06 Nov 2011 07:18:13 +0000
Received: from eris.apache.org (localhost [127.0.0.1])
	by eris.apache.org (Postfix) with ESMTP id 7DF0A2388A6E
	for <commits@commons.apache.org>; Sun,  6 Nov 2011 07:17:51 +0000 (UTC)
Content-Type: text/plain; charset="utf-8"
MIME-Version: 1.0
Content-Transfer-Encoding: 8bit
Subject: svn commit: r1198165 - in /commons/proper/math/trunk: ./
 src/main/java/org/apache/commons/math/exception/util/
 src/main/java/org/apache/commons/math/random/
 src/main/resources/META-INF/localization/ src/site/xdoc/
 src/test/java/org/apache/commons/math/...
Date: Sun, 06 Nov 2011 07:17:51 -0000
To: commits@commons.apache.org
From: psteitz@apache.org
X-Mailer: svnmailer-1.0.8-patched
Message-Id: <20111106071751.7DF0A2388A6E@eris.apache.org>
X-Virus-Checked: Checked by ClamAV on apache.org

Author: psteitz
Date: Sun Nov  6 07:17:50 2011
New Revision: 1198165

URL: http://svn.apache.org/viewvc?rev=1198165&view=rev
Log:
Added stable random generator based on Chambers-Mallows-Stuck method.  Contributed by Pavel Ryhzov.  JIRA: MATH-462.

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/StableRandomGenerator.java   (with props)
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/random/StableRandomGeneratorTest.java   (with props)
Modified:
    commons/proper/math/trunk/pom.xml
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java
    commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties
    commons/proper/math/trunk/src/site/xdoc/changes.xml

Modified: commons/proper/math/trunk/pom.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/pom.xml?rev=1198165&r1=1198164&r2=1198165&view=diff
==============================================================================
--- commons/proper/math/trunk/pom.xml (original)
+++ commons/proper/math/trunk/pom.xml Sun Nov  6 07:17:50 2011
@@ -225,6 +225,9 @@
       <name>Matthew Rowles</name>
     </contributor>
     <contributor>
+      <name>Pavel Ryzhov</name>
+    </contributor>
+    <contributor>
       <name>Joni Salonen</name>
     </contributor>
     <contributor>

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java?rev=1198165&r1=1198164&r2=1198165&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java Sun Nov  6 07:17:50 2011
@@ -271,6 +271,8 @@ public enum LocalizedFormats implements 
     OUT_OF_RANGE_ROOT_OF_UNITY_INDEX("out of range root of unity index {0} (must be in [{1};{2}])"),
     OUT_OF_RANGE("out of range"), /* keep */
     OUT_OF_RANGE_SIMPLE("{0} out of [{1}, {2}] range"), /* keep */
+    OUT_OF_RANGE_LEFT("{0} out of ({1}, {2}] range"),
+    OUT_OF_RANGE_RIGHT("{0} out of [{1}, {2}) range"),
     OUTLINE_BOUNDARY_LOOP_OPEN("an outline boundary loop is open"),
     OVERFLOW("overflow"), /* keep */
     OVERFLOW_IN_FRACTION("overflow in fraction {0}/{1}, cannot negate"),

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/StableRandomGenerator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/StableRandomGenerator.java?rev=1198165&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/StableRandomGenerator.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/StableRandomGenerator.java Sun Nov  6 07:17:50 2011
@@ -0,0 +1,131 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.random;
+
+import org.apache.commons.math.exception.NullArgumentException;
+import org.apache.commons.math.exception.OutOfRangeException;
+import org.apache.commons.math.exception.util.LocalizedFormats;
+import org.apache.commons.math.util.FastMath;
+
+/**
+ * <p>This class provides a stable normalized random generator. It samples from a stable
+ * distribution with location parameter 0 and scale 1.</p>
+ *
+ * <p>The implementation uses the Chambers-Mallows-Stuck method as described in
+ * <i>Handbook of computational statistics: concepts and methods</i> by
+ * James E. Gentle, Wolfgang HŠrdle, Yuichi Mori.</p>
+ *
+ * @version $Id$
+ * @since 3.0
+ */
+public class StableRandomGenerator implements NormalizedRandomGenerator {
+    /** Underlying generator. */
+    private final RandomGenerator generator;
+
+    /** stability parameter */
+    private final double alpha;
+
+    /** skewness parameter */
+    private final double beta;
+
+    /** cache of expression value used in generation */
+    private final double zeta;
+
+    /**
+     * Create a new generator.
+     *
+     * @param generator underlying random generator to use
+     * @param alpha Stability parameter. Must be in range (0, 2]
+     * @param beta Skewness parameter. Must be in range [-1, 1]
+     */
+    public StableRandomGenerator(final RandomGenerator generator, double alpha,
+            double beta) {
+        if (generator == null) {
+            throw new NullArgumentException();
+        }
+
+        if (!(alpha > 0d && alpha <= 2d)) {
+            throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_LEFT,
+                    alpha, 0, 2);
+        }
+
+        if (!(beta >= -1d && beta <= 1d)) {
+            throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_SIMPLE,
+                    beta, -1, 1);
+        }
+
+        this.generator = generator;
+        this.alpha = alpha;
+        this.beta = beta;
+        if (alpha < 2d && beta != 0d) {
+            zeta = beta * FastMath.tan(FastMath.PI * alpha / 2);
+        } else {
+            zeta = 0d;
+        }
+    }
+
+    /**
+     * Generate a random scalar with zero location and unit scale.
+     *
+     * @return a random scalar with zero location and unit scale
+     */
+    public double nextNormalizedDouble() {
+        // we need 2 uniform random numbers to calculate omega and phi
+        double omega = -FastMath.log(generator.nextDouble());
+        double phi = FastMath.PI * (generator.nextDouble() - 0.5);
+
+        // Normal distribution case (Box-Muller algorithm)
+        if (alpha == 2d) {
+            return FastMath.sqrt(2d * omega) * FastMath.sin(phi);
+        }
+
+        double x;
+        // when beta = 0, zeta is zero as well
+        // Thus we can exclude it from the formula
+        if (beta == 0d) {
+            // Cauchy distribution case
+            if (alpha == 1d) {
+                x = FastMath.tan(phi);
+            } else {
+                x = FastMath.pow(omega * FastMath.cos((1 - alpha) * phi),
+                    1d / alpha - 1d) *
+                    FastMath.sin(alpha * phi) /
+                    FastMath.pow(FastMath.cos(phi), 1d / alpha);
+            }
+        } else {
+            // Generic stable distribution
+            double cosPhi = FastMath.cos(phi);
+            // to avoid rounding errors around alpha = 1
+            if (FastMath.abs(alpha - 1d) > 1e-8) {
+                double alphaPhi = alpha * phi;
+                double invAlphaPhi = phi - alphaPhi;
+                x = (FastMath.sin(alphaPhi) + zeta * FastMath.cos(alphaPhi)) / cosPhi *
+                    (FastMath.cos(invAlphaPhi) + zeta * FastMath.sin(invAlphaPhi)) /
+                     FastMath.pow(omega * cosPhi, (1 - alpha) / alpha);
+            } else {
+                double betaPhi = FastMath.PI / 2 + beta * phi;
+                x = 2d / FastMath.PI * (betaPhi * FastMath.tan(phi) - beta *
+                    FastMath.log(FastMath.PI / 2d * omega * cosPhi / betaPhi));
+
+                if (alpha != 1d) {
+                    x = x + beta * FastMath.tan(FastMath.PI * alpha / 2);
+                }
+            }
+        }
+        return x;
+    }
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/StableRandomGenerator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/StableRandomGenerator.java
------------------------------------------------------------------------------
    svn:keywords = Date Revision Id

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/StableRandomGenerator.java
------------------------------------------------------------------------------
    svn:mime-type = text/plain

Modified: commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties?rev=1198165&r1=1198164&r2=1198165&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties (original)
+++ commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties Sun Nov  6 07:17:50 2011
@@ -237,6 +237,8 @@ SIGNIFICANCE_LEVEL = niveau de significa
 OUT_OF_ORDER_ABSCISSA_ARRAY = les abscisses doivent \u00eatre en ordre strictement croissant, mais l''\u00e9l\u00e9ment {0} vaut {1} alors que l''\u00e9l\u00e9ment {2} vaut {3}
 OUT_OF_RANGE_ROOT_OF_UNITY_INDEX = l''indice de racine de l''unit\u00e9 {0} est hors du domaine autoris\u00e9 [{1};{2}]
 OUT_OF_RANGE_SIMPLE = {0} hors du domaine [{1}, {2}]
+OUT_OF_RANGE_LEFT = {0} hors du domaine ({1}, {2}]
+OUT_OF_RANGE_RIGHT = {0} hors du domaine [{1}, {2})
 OUT_OF_RANGE = hors domaine
 OVERFLOW = d\u00e9passement de capacit\u00e9
 OVERFLOW_IN_FRACTION = d\u00e9passement de capacit\u00e9 pour la fraction {0}/{1}, son signe ne peut \u00eatre chang\u00e9

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1198165&r1=1198164&r2=1198165&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Sun Nov  6 07:17:50 2011
@@ -52,6 +52,9 @@ The <action> type attribute can be add,u
     If the output is not quite correct, check for invisible trailing spaces!
      -->
     <release version="3.0" date="TBD" description="TBD">
+      <action dev="psteitz" type="add" issue="MATH-462" due-to="Pavel Ryzhof">
+        Added stable random generator based on Chambers-Mallows-Stuck method.
+      </action>
       <action dev="psteitz" type="update" issue="MATH-701">
         Changed the default seed used for RandomDataImpl, AbstractWell and MersenneTwister
         PRNGs to add the system identity hash code of the instance to the current system

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math/random/StableRandomGeneratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/random/StableRandomGeneratorTest.java?rev=1198165&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/random/StableRandomGeneratorTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/random/StableRandomGeneratorTest.java Sun Nov  6 07:17:50 2011
@@ -0,0 +1,128 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.random;
+
+import junit.framework.TestCase;
+
+import org.apache.commons.math.exception.OutOfRangeException;
+import org.apache.commons.math.stat.StatUtils;
+import org.apache.commons.math.stat.descriptive.DescriptiveStatistics;
+
+/**
+ * The class <code>StableRandomGeneratorTest</code> contains tests for the class
+ * {@link <code>StableRandomGenerator</code>}
+ * 
+ * @version $Revision$
+ */
+public class StableRandomGeneratorTest extends TestCase {
+
+    private RandomGenerator rg = new Well19937c(100);
+    private final static int sampleSize = 10000;
+
+    /**
+     * Construct new test instance
+     * 
+     * @param name the test name
+     */
+    public StableRandomGeneratorTest(String name) {
+        super(name);
+    }
+
+    /**
+     * Run the double nextDouble() method test Due to leptokurtic property the
+     * acceptance range is widened.
+     */
+    public void testNextDouble() {
+        StableRandomGenerator generator = new StableRandomGenerator(rg, 1.3,
+                0.1);
+        double[] sample = new double[2 * sampleSize];
+        for (int i = 0; i < sample.length; ++i) {
+            sample[i] = generator.nextNormalizedDouble();
+        }
+        assertEquals(0.0, StatUtils.mean(sample), 0.2);
+    }
+
+    /**
+     * If alpha = 2, than it must be Gaussian distribution
+     */
+    public void testGaussianCase() {
+        StableRandomGenerator generator = new StableRandomGenerator(rg, 2d, 0.0);
+
+        double[] sample = new double[sampleSize];
+        for (int i = 0; i < sample.length; ++i) {
+            sample[i] = generator.nextNormalizedDouble();
+        }
+        assertEquals(0.0, StatUtils.mean(sample), 0.02);
+        assertEquals(1.0, StatUtils.variance(sample), 0.02);
+    }
+
+    /**
+     * If alpha = 1, than it must be Cauchy distribution
+     */
+    public void testCauchyCase() {
+        StableRandomGenerator generator = new StableRandomGenerator(rg, 1d, 0.0);
+        DescriptiveStatistics summary = new DescriptiveStatistics();
+
+        for (int i = 0; i < sampleSize; ++i) {
+            double sample = generator.nextNormalizedDouble();
+            summary.addValue(sample);
+        }
+
+        // Standard Cauchy distribution should have zero median and mode
+        double median = summary.getPercentile(50);
+        assertEquals(0.0, median, 0.2);
+    }
+
+    /**
+     * Input parameter range tests
+     */
+    public void testAlphaRangeBelowZero() {
+        try {
+            StableRandomGenerator generator = new StableRandomGenerator(rg,
+                    -1.0, 0.0);
+        } catch (OutOfRangeException e) {
+            assertEquals(-1.0, e.getArgument());
+        }
+    }
+
+    public void testAlphaRangeAboveTwo() {
+        try {
+            StableRandomGenerator generator = new StableRandomGenerator(rg,
+                    3.0, 0.0);
+        } catch (OutOfRangeException e) {
+            assertEquals(3.0, e.getArgument());
+        }
+    }
+
+    public void testBetaRangeBelowMinusOne() {
+        try {
+            StableRandomGenerator generator = new StableRandomGenerator(rg,
+                    1.0, -2.0);
+        } catch (OutOfRangeException e) {
+            assertEquals(-2.0, e.getArgument());
+        }
+    }
+
+    public void testBetaRangeAboveOne() {
+        try {
+            StableRandomGenerator generator = new StableRandomGenerator(rg,
+                    1.0, 2.0);
+        } catch (OutOfRangeException e) {
+            assertEquals(2.0, e.getArgument());
+        }
+    }
+}
\ No newline at end of file

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math/random/StableRandomGeneratorTest.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math/random/StableRandomGeneratorTest.java
------------------------------------------------------------------------------
    svn:keywords = Date Revision Id

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math/random/StableRandomGeneratorTest.java
------------------------------------------------------------------------------
    svn:mime-type = text/plain