Results 1 to 3 of 3

Thread: Bogus results from FNL7.0 Stat Library Function BETIN()

  1. #1
    Senior Member mecej4's Avatar
    Join Date
    Dec 2009
    Posts
    139

    Thumbs up Bogus results from FNL7.0 Stat Library Function BETIN()

    The function BETIN(Q, p, q) returns the value of the inverse of the regularized beta function, I_x(p,q), using the notation of https://en.wikipedia.org/wiki/Beta_distribution . For large values of the parameter q, the routine from IMSL FNL 7.0 returns inaccurate or even completely bogus values.

    The FNL documentation for FIN, the corresponding inverse for the F-distribution, states that FIN(Q, p1, q1) is obtained using the well-known transformation (see https://en.wikipedia.org/wiki/F-distribution ) from the results from BETIN. Curiously, if I use the transformation in reverse, by calling FIN first, I obtain reasonable, if not always accurate, results, for the beta-inverse values.

    The following Fortran code will produce illustrative results.

    Code:
          program xibeta
          USE BETIN_INT
          USE FIN_INT
          IMPLICIT   NONE
          REAL       QIN(7), X1, X2, Y, Quant,d1,d2
          integer i,j
          data QIN/1E3,2E3,5E3,1E4,2E4,5E4,1E5/
    !
          Quant=0.9999                        ! quantile value
          d1=2.0                              ! PIN of BETIN
          DO I=1,7
             d2 = QIN(i)                      ! QIN of BETIN
             X1 = BETIN(Quant,d1,d2)          ! Result from BETIN
             Y  = FIN(Quant,2.0*d1,2.0*d2)
             X2 = 1.0/(1.0+d2/(Y*d1))         ! Result from FIN
             WRITE (*,10)d2,X1,X2
          END DO
       10 format(F7.0,2x,1p2E12.5)
          END
    Here are the results from FNL7, for p ("PIN") = 2.0 and various large values of q ("QIN") for the quantile value Q = 0.9999, with an additional column containing the results from Maple-17 (command: with(Statistics): Quantile(('Beta')(p, q), Q); )

    Code:
       q      BETIN (FNL)  From FIN (FNL)   Maple-17
       
      1000.   1.16814E-2    1.16815E-2    1.16817E-2
      2000.   5.85949E-3    5.85941E-3    5.85948E-3
      5000.   2.34829E-3    2.34814E-3    2.34828E-3
     10000.   1.17485E-3    1.17498E-3    1.17489E-3
     20000.   4.33734E-3    5.87595E-4    5.87631E-4
     50000.   3.76282E-1    2.35025E-4    2.35097E-4
    100000.   3.79713E-1    1.17527E-4    1.17556E-4
    The last two values in the second column are bogus. Earlier values in the same column are reasonable, but not accurate.
    Last edited by mecej4; 10-06-2016 at 03:39 AM.

  2. #2
    Junior Member
    Join Date
    Feb 2009
    Location
    Houston, TX
    Posts
    29
    Please report this problem to customer support. We will check if this exists in IMSL Fortran 7.1. If the issue occurs in 7.1 we will create a bug report to fix in the upcoming release. By reporting directly to customer support we will contact you directly to let you know the issue is fixed.
    Last edited by jlo; 10-20-2016 at 04:20 PM.
    Jennifer Locke
    Product Manager, IMSL, Views, JViews, Elixir and Stingray
    P: 713.954.6454

  3. #3
    Administrator
    Join Date
    Feb 2008
    Posts
    24
    When the included source code is converted to double precision BETIN returns
    results identical to FIN and Maple-17. It is our recommendation that the
    double precision version of BETIN be used.


    Code:
        program xibeta
        USE BETIN_INT
        USE FIN_INT
        IMPLICIT NONE
    
    
        REAL(kind(1d0)) QIN(7), X1, X2, Y, Quant, d1, d2
        integer i,j
        data QIN/1.d3, 2.d3, 5.d3, 1.d4, 2.d4, 5.d4, 1.d5/
    !
        Quant=0.9999d0 ! quantile value
        d1=2.0d0 ! PIN of BETIN
        DO I=1,7
            d2 = QIN ! QIN of BETIN
            X1 = BETIN(Quant, d1, d2) ! Result from BETIN
            Y = FIN(Quant, 2.0d0*d1, 2.0d0*d2)
            X2 = 1.0d0 / (1.0d0 + d2 / (Y * d1)) ! Result from FIN
            WRITE (*,10) d2, X1, X2
        END DO
      10 format(F7.0, 2x, 1p2E20.8)
        END
    Here are the results:

    Code:
    Double Precision Results
    =============================================
      1000.   1.16817E-02 1.16817E-02
      2000.   5.85948E-03 5.85948E-03
      5000.   2.34828E-03 2.34828E-03
     10000.   1.17489E-03 1.17489E-03
     20000.   5.87631E-04 5.87631E-04
     50000.   2.35097E-04 2.35097E-04
    100000.   1.17556E-04 1.17556E-04
    =============================================
    Jeremy Dean
    Rogue Wave Software
    Technical Support
    support@roguewave.com

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •