View Javadoc
1   /*
2    * Licensed to the Hipparchus project under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.hipparchus.distribution.discrete;
18  
19  import org.junit.Assert;
20  import org.junit.Test;
21  
22  public class SaddlePointExpansionTest {
23  
24      @Test
25      public void testSmallInteger() {
26          for (int n = 4; n < 16; ++n) {
27               Assert.assertEquals(alternateStirlingErrorImplementation(n, 8),
28                                  SaddlePointExpansion.getStirlingError(n),
29                                  1.0e-10);
30          }
31      }
32  
33      @Test
34      public void testSmallNonInteger() {
35          for (double z = 3.75; z < 14.8; z += 1.0) {
36              Assert.assertEquals(alternateStirlingErrorImplementation(z, 12),
37                                  SaddlePointExpansion.getStirlingError(z),
38                                  1.0e-10);
39          }
40      }
41  
42      @Test
43      public void testLargeValues() {
44          for (double z = 15.25; z < 25.5; z += 1.0) {
45              Assert.assertEquals(alternateStirlingErrorImplementation(z, 21),
46                                  SaddlePointExpansion.getStirlingError(z),
47                                  1.0e-15);
48          }
49      }
50  
51      @Test
52      public void testSpecialValues() {
53          Assert.assertEquals(0.0, SaddlePointExpansion.logBinomialProbability(0, 0, 0.6, 0.4), 1.0e-15);
54          Assert.assertEquals(Double.NEGATIVE_INFINITY, SaddlePointExpansion.logBinomialProbability(1, 0, 0.6, 0.4), 1.0e-15);
55      }
56  
57      private double alternateStirlingErrorImplementation(final double z, final int kMax) {
58          // numerators of coefficients in Stirling's expansion for log(Γ(z)) from https://oeis.org/A046968
59          final double[] A046968 = {
60              1.0, -1.0, 1.0, -1.0, 1.0, -691.0, 1.0,
61              -3617.0, 43867.0, -174611.0, 77683.0, -236364091.0, 657931.0, -3392780147.0,
62              1723168255201.0, -7709321041217.0, 151628697551.0, -26315271553053477373.0,
63              154210205991661.0, -261082718496449122051.0, 1520097643918070802691.0
64          };
65          // denominators of coefficients in Stirling's expansion for log(Γ(z)) from https://oeis.org/A046969
66          final double[] A046969 = {
67              12.0, 360.0, 1260.0, 1680.0, 1188.0, 360360.0, 156.0,
68              122400.0, 244188.0, 125400.0, 5796.0, 1506960.0, 300.0, 93960.0,
69              2492028.0, 505920.0, 396.0, 2418179400.0, 444.0, 21106800.0, 3109932.0
70          };
71          double zk   = z;
72          double expansion = 0.0;
73          for (int i = 0; i < kMax; ++i) {
74              expansion += A046968[i] / (zk * A046969[i]);
75              zk *= z * z;
76          }
77          return expansion;
78      }
79  }