IdentifiantMot de passe
Loading...
Mot de passe oublié ?Je m'inscris ! (gratuit)
Navigation

Inscrivez-vous gratuitement
pour pouvoir participer, suivre les réponses en temps réel, voter pour les messages, poser vos propres questions et recevoir la newsletter

Fortran Discussion :

Précision et utilisation fonction gamma


Sujet :

Fortran

  1. #1
    Membre régulier
    Homme Profil pro
    Inscrit en
    Juillet 2006
    Messages
    79
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : Belgique

    Informations forums :
    Inscription : Juillet 2006
    Messages : 79
    Points : 89
    Points
    89
    Par défaut Précision et utilisation fonction gamma
    Bonjour,

    j'ai un programme dans lequel je cherche à calculer une p-value associée à une loi gamma.
    J'utilise des routines classiques pour ce faire. cf l'ECM.

    Mon problème est que pour certaines valeurs de paramètres (e.g. celles des lignes 6-8) j'obtiens une valeur de 1 à l'appel de gammad.J'ai longtemps pensé que le problème venait de l'impossibilité de stocker des valeurs avec une précision supérieure. Après une exécution avec gdb, je me suis aperçu que j'avais des calculs avec une précision très élevée, jusqu'à la ligne 130, où l'opération :
    one -EXP(ARG)
    retourne 1 ! alors que la représentation de exp(ARG) est différente de 0
    (j'ai ajouté un write (ligne 132), pour connaître la valeur de exp(ARG) j'ai une valeur d'environ 4.32E-32).

    Je n'arrive pas bien à comprendre ce qui se passe. Si quelqu'un avait une idée et un moyen de conserver des valeurs de gamma avec de grandes précisions, je lui en serais très reconnaissant.

    Par avance merci

    Code : Sélectionner tout - Visualiser dans une fenêtre à part
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    131
    132
    133
    134
    135
    136
    137
    138
    139
    140
    141
    142
    143
    144
    145
    146
    147
    148
    149
    150
    151
    152
    153
    154
    155
    156
    157
    158
    159
    160
    161
    162
    163
    164
    165
    166
    167
    168
    169
    170
    171
    172
    173
    174
    175
    176
    177
    178
    179
    180
    181
    182
    183
    184
    185
    186
    187
    188
    189
    190
    191
    192
    193
    194
    195
    196
    197
    198
    199
    200
    201
    202
    203
    204
    205
    206
    207
    208
    209
    210
    211
    212
    213
    214
    215
    216
    217
    218
    219
    220
    221
    222
    223
    224
    225
    226
    227
    228
    229
    230
    231
    232
    233
    234
    235
    236
    237
    238
    239
    240
    241
    242
    243
    244
    245
    246
    247
    248
    249
    250
    251
    252
    253
    254
    255
    256
    257
    258
    259
    260
    261
    262
    263
    264
    265
    266
    267
    268
    269
    270
    271
    272
    273
    274
    275
    276
    277
    278
    279
    280
    281
    282
    283
    284
    285
    286
    287
    288
    289
    290
    291
    292
    293
    294
    295
    296
    297
    298
    299
    300
    301
    302
    303
    304
    305
    306
    307
    308
    309
    310
    311
    312
    313
    314
    315
    316
    317
    318
    319
    320
    321
    322
    323
    324
    325
    326
    327
    328
    329
    330
    331
    332
    333
    334
    335
    336
    337
    338
     
    program testGamma
    implicit None
    integer,parameter::dp=selected_real_kind(10,30)
    real ( kind = dp ) bound,scale,shape,p,q,score
     
    shape=6
    scale=20
    score=1800
     
    write(*,*)(score/scale),shape,"  ---  ",1.0_dp-gammad(score/scale,shape)
     
    contains 
     
      !======================================================================
      ! Gamma density function
      !======================================================================
     
      FUNCTION gammad(x, p) RESULT(fn_val)
     
        !      ALGORITHM AS239  APPL. STATIST. (1988) VOL. 37, NO. 3
     
        !      Computation of the Incomplete Gamma Integral
     
        !      Auxiliary functions required: ALOGAM = logarithm of the gamma
        !      function, and ALNORM = algorithm AS66
     
        ! ELF90-compatible version by Alan Miller
        ! Latest revision - 27 October 2000
     
        ! N.B. Argument IFAULT has been removed
     
        IMPLICIT NONE
        INTEGER, PARAMETER    :: dp = SELECTED_REAL_KIND(12, 60)
        REAL (dp), INTENT(IN) :: x, p
        REAL (dp)             :: fn_val
     
        ! Local variables
        REAL (dp)             :: pn1, pn2, pn3, pn4, pn5, pn6, arg, c, rn, a, b, an
        REAL (dp), PARAMETER  :: zero = 0.d0, one = 1.d0, two = 2.d0, &
             oflo = 1.d+37, three = 3.d0, nine = 9.d0, &
             tol = 1.d-20, xbig = 1.d+8, plimit = 1000.d0, &
             elimit = -88.d0
        INTEGER               :: ifault
        ! EXTERNAL alogam, alnorm
     
        fn_val = zero
     
        !      Check that we have valid values for X and P
     
        IF (p <= zero .OR. x < zero) THEN
           WRITE(*, *) 'AS239: Either p <= 0 or x < 0'
           RETURN
        END IF
        IF (x == zero) RETURN
     
        !      Use a normal approximation if P > PLIMIT
     
        IF (p > plimit) THEN
           pn1 = three * SQRT(p) * ((x / p) ** (one / three) + one /(nine * p) - one)
           fn_val = alnorm(pn1, .false.)
           RETURN
        END IF
     
        !      If X is extremely large compared to P then set fn_val = 1
     
        IF (x > xbig) THEN
           fn_val = one
           RETURN
        END IF
     
        IF (x <= one .OR. x < p) THEN
     
           !      Use Pearson's series expansion.
           !      (Note that P is not large enough to force overflow in ALOGAM).
           !      No need to test IFAULT on exit since P > 0.
     
           arg = p * LOG(x) - x - alogam(p + one, ifault)
           c = one
           fn_val = one
           a = p
    40     a = a + one
           c = c * x / a
           fn_val = fn_val + c
           IF (c > tol) GO TO 40
           arg = arg + LOG(fn_val)
           fn_val = zero
           IF (arg >= elimit) fn_val = EXP(arg)
     
        ELSE
     
           !      Use a continued fraction expansion
     
           arg = p * LOG(x) - x - alogam(p, ifault)
           a = one - p
           b = a + x + one
           c = zero
           pn1 = one
           pn2 = x
           pn3 = x + one
           pn4 = x * b
           fn_val = pn3 / pn4
    60     a = a + one
           b = b + two
           c = c + one
           an = a * c
           pn5 = b * pn3 - an * pn1
           pn6 = b * pn4 - an * pn2
           IF (ABS(pn6) > zero) THEN
              rn = pn5 / pn6
              IF (ABS(fn_val - rn) <= MIN(tol, tol * rn)) GO TO 80
              fn_val = rn
           END IF
     
           pn1 = pn3
           pn2 = pn4
           pn3 = pn5
           pn4 = pn6
           IF (ABS(pn5) >= oflo) THEN
     
              !      Re-scale terms in continued fraction if terms are large
     
              pn1 = pn1 / oflo
              pn2 = pn2 / oflo
              pn3 = pn3 / oflo
              pn4 = pn4 / oflo
           END IF
           GO TO 60
    80     arg = arg + LOG(fn_val)
           fn_val = one
           IF (arg >= elimit) fn_val = one - EXP(arg)
        END IF
        write(*,*)EXP(arg),fn_val
        RETURN
      END FUNCTION gammad
     
     
      function alogam ( x, ifault )
     
        !*****************************************************************************80
        !
        !! ALOGAM computes the logarithm of the Gamma function.
        !
        !  Modified:
        !
        !    28 March 1999
        !
        !  Author:
        !
        !    Original FORTRAN77 version by Malcolm Pike, David Hill.
        !    FORTRAN90 version by John Burkardt.
        !
        !  Reference:
        !
        !    Malcolm Pike, David Hill,
        !    Algorithm 291:
        !    Logarithm of Gamma Function,
        !    Communications of the ACM,
        !    Volume 9, Number 9, September 1966, page 684.
        !
        !  Parameters:
        !
        !    Input, real ( kind = 8 ) X, the argument of the Gamma function.
        !    X should be greater than 0.
        !
        !    Output, integer ( kind = 4 ) IFAULT, error flag.
        !    0, no error.
        !    1, X <= 0.
        !
        !    Output, real ( kind = 8 ) ALOGAM, the logarithm of the Gamma
        !    function of X.
        !
        implicit none
        INTEGER, PARAMETER    :: dp = SELECTED_REAL_KIND(12, 60)
        real    ( kind = dp ) alogam
        real    ( kind = dp ) f
        integer :: ifault
        real    ( kind = dp ) x
        real    ( kind = dp ) y
        real    ( kind = dp ) z
     
        if ( x <= 0.0D+00 ) then
           ifault = 1
           alogam = 0.0D+00
           return
        end if
     
        ifault = 0
        y = x
     
        if ( x < 7.0D+00 ) then
     
           f = 1.0D+00
           z = y
     
           do while ( z < 7.0D+00 )
              f = f * z
              z = z + 1.0D+00
           end do
     
           y = z
           f = - log ( f )
     
        else
     
           f = 0.0D+00
     
        end if
     
        z = 1.0D+00 / y / y
     
        alogam = f + ( y - 0.5D+00 ) * log ( y ) - y &
             + 0.918938533204673D+00 + &
             ((( &
             - 0.000595238095238D+00   * z &
             + 0.000793650793651D+00 ) * z &
             - 0.002777777777778D+00 ) * z &
             + 0.083333333333333D+00 ) / y
     
        return
      end function alogam
     
      !  Algorithm AS66 Applied Statistics (1973) vol.22, no.3
     
      !  Evaluates the tail area of the standardised normal curve
      !  from x to infinity if upper is .true. or
      !  from minus infinity to x if upper is .false.
     
      ! ELF90-compatible version by Alan Miller
      ! Latest revision - 29 November 2001
     
      function alnorm ( x, upper )
     
        !*****************************************************************************80
        !
        !! ALNORM computes the cumulative density of the standard normal distribution.
        !
        !  Licensing:
        !
        !    This code is distributed under the GNU LGPL license.
        !
        !  Modified:
        !
        !    17 January 2008
        !
        !  Author:
        !
        !    Original FORTRAN77 version by David Hill
        !    MATLAB version by John Burkardt
        !
        !  Reference:
        !
        !    David Hill,
        !    Algorithm AS 66:
        !    The Normal Integral,
        !    Applied Statistics,
        !    Volume 22, Number 3, 1973, pages 424-427.
        !
        !  Parameters:
        !
        !    Input, real X, is one endpoint of the semi-infinite interval
        !    over which the integration takes place.
        !
        !    Input, logical UPPER, determines whether the upper or lower
        !    interval is to be integrated:
        !    1 => integrate from X to + Infinity;
        !    0 => integrate from - Infinity to X.
        !
        !    Output, real VALUE, the integral of the standard normal
        !    distribution over the desired interval.
        !
        implicit none
        INTEGER, PARAMETER    :: dp = SELECTED_REAL_KIND(12, 60)
        real(kind=dp), parameter ::  a1 = 5.75885480458
        real(kind=dp), parameter ::  a2 = 2.62433121679
        real(kind=dp), parameter ::  a3 = 5.92885724438
        real(kind=dp), parameter ::  b1 = -29.8213557807
        real(kind=dp), parameter ::  b2 = 48.6959930692
        real(kind=dp), parameter ::  c1 = -0.000000038052
        real(kind=dp), parameter ::  c2 = 0.000398064794
        real(kind=dp), parameter ::  c3 = -0.151679116635
        real(kind=dp), parameter ::  c4 = 4.8385912808
        real(kind=dp), parameter ::  c5 = 0.742380924027
        real(kind=dp), parameter ::  c6 = 3.99019417011
        real(kind=dp), parameter ::  con = 1.28
        real(kind=dp), parameter ::  d1 = 1.00000615302
        real(kind=dp), parameter ::  d2 = 1.98615381364
        real(kind=dp), parameter ::  d3 = 5.29330324926
        real(kind=dp), parameter ::  d4 = -15.1508972451
        real(kind=dp), parameter ::  d5 = 30.789933034
        real(kind=dp), parameter ::  ltone = 7.0
        real(kind=dp), parameter ::  p = 0.39894228044
        real(kind=dp), parameter ::  q = 0.39990348504
        real(kind=dp), parameter ::  r = 0.398942280385
        real(kind=dp), parameter ::  utzero = 18.66
        real(kind=dp) ::x,y,z,alnorm
        logical ::up,upper
     
        up = upper
        z = x
     
        if (z .lt. 0.0D+00 )then
           up = .not. up
           z = - z
        end if
     
        if ( z .gt. ltone .and. ( ( .not. up ) .or. utzero .lt. z ) ) then
           if ( up ) then
              alnorm = 0.0D+00
           else
              alnorm = 1.0D+00
           end if
     
           return
     
        end if
     
        y = 0.5D+00 * z * z
     
        if ( z .le. con ) then
     
           alnorm = 0.5D+00 - z * (p-q*y/(y+a1+b1/(y+a2+b2/( y + a3 ))))
     
        else
     
           alnorm = r * exp ( - y )/(z+c1+d1/(z+c2+d2/(z+c3+d3/(z+c4+d4/(z+c5+d5/(z+c6))))))
     
        end if
     
        if ( .not. up ) then
           alnorm = 1.0D+00 - alnorm
        end if
     
        return
      end function alnorm
     
     
    end program testGamma
    Notes supplémentaires:

    Je ne travaille que sous Linux 3.0.0 x86_64), je n'arrive à compiler l'ECM qu'avec gfortran 4.6 (ifort12.1 ne semble pas aimer le write dans dans la fonction, et sans cette instruction je ne sais pas comment est représenté exp(ARG) à la ligne 130)

    Pour des valeurs plus faibles des paramètres de gammad,
    shape=6
    scale=20
    score=1000
    la fonction gammad marche encore bien et donne des précisions de l'ordre de 1E-16.

  2. #2
    Modérateur

    Profil pro
    Inscrit en
    Août 2006
    Messages
    974
    Détails du profil
    Informations personnelles :
    Localisation : Canada

    Informations forums :
    Inscription : Août 2006
    Messages : 974
    Points : 1 346
    Points
    1 346
    Par défaut
    En double précision, tu as environ 15 chiffres de précision. Alors, il est tout à fait normal d'avoir 1 - exp(ARG) égal à 1 si exp(Arg) vaut 4.32E-32. La résultat 0,999 999 999 999 999 999 999 999 999 999 n'est simplement (ou doublement !) pas représentable.

    Tu peux essayer en précision quadruple, mais 1 - 10^-32 est proche de la limite. Voir http://en.wikipedia.org/wiki/Quadrup...g-point_format

  3. #3
    Membre régulier
    Homme Profil pro
    Inscrit en
    Juillet 2006
    Messages
    79
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : Belgique

    Informations forums :
    Inscription : Juillet 2006
    Messages : 79
    Points : 89
    Points
    89
    Par défaut
    Merci bien Sylvain,

    ceci confirme ma première impression. En fait ce qui m'a induit en erreur c'est qu'avec la ligne 132:

    Code : Sélectionner tout - Visualiser dans une fenêtre à part
     write(*,*)EXP(arg),fn_val
    J'arrivais a afficher des valeurs nettement plus faibles (de l'ordre de 10E-250) ! Le biais de raisonnement c'est que cette valeur est affichable mais pas forcément représentable en machine.
    J'avais un vague espoir qu'avec une invocation des modules IEEE, j'aurais une solution standard et plus propre. Mais a priori non, c'est cela ?

    Me reste deux solutions, ou je stocke une valeur intermédiaire représentable et ne fait la transformation qu'au moment de l'écriture....ou j'indique que les solutions sont bornées à la précision permise par un réel quadruple précision !

    En tout une fois de plus merci beaucoup Sylvain.

  4. #4
    Modérateur

    Profil pro
    Inscrit en
    Août 2006
    Messages
    974
    Détails du profil
    Informations personnelles :
    Localisation : Canada

    Informations forums :
    Inscription : Août 2006
    Messages : 974
    Points : 1 346
    Points
    1 346
    Par défaut
    Citation Envoyé par yogitetradim Voir le message
    ...En fait ce qui m'a induit en erreur c'est qu'avec la ligne 132:

    Code : Sélectionner tout - Visualiser dans une fenêtre à part
     write(*,*)EXP(arg),fn_val
    J'arrivais a afficher des valeurs nettement plus faibles (de l'ordre de 10E-250) ! Le biais de raisonnement c'est que cette valeur est affichable mais pas forcément représentable en machine.
    Le problème n'est pas la représentation interne d'un (très) petit nombre (si il est affichable, c'est qu'il est représentable). C'est le résultat de la différence entre un nombre très grand (1) et très petit (10^-32) qui n'est pas représentable sans perte de précision. La perte étant en fait la différence par rapport au nombre très grand. Tu peux le tester sur ta calculette : 1 - 10^-32 retourne 1, pas 0,9999999999999999...

    Si tu tiens à conserver un tel résultat, tu devras soit augmenter la précision, soit retourner le résultat intermédiaire (exp(Arg)).

Discussions similaires

  1. Pb de blancs dans utilisation fonction FtpFindFirstFile
    Par AlvinTheMaker dans le forum MFC
    Réponses: 2
    Dernier message: 06/04/2005, 13h33
  2. Réponses: 6
    Dernier message: 24/02/2005, 10h44
  3. [GIMP] [Script-FU] Utilisation fonction gimp-curves-spline
    Par narmataru dans le forum Autres langages
    Réponses: 1
    Dernier message: 09/02/2005, 18h25
  4. [Débutant] Aide utilisation fonctions :(
    Par trakiss dans le forum Débuter
    Réponses: 10
    Dernier message: 27/08/2004, 16h59
  5. Utilisation fonction définie dans un .Dll
    Par jeab. dans le forum Windows
    Réponses: 5
    Dernier message: 23/03/2004, 17h23

Partager

Partager
  • Envoyer la discussion sur Viadeo
  • Envoyer la discussion sur Twitter
  • Envoyer la discussion sur Google
  • Envoyer la discussion sur Facebook
  • Envoyer la discussion sur Digg
  • Envoyer la discussion sur Delicious
  • Envoyer la discussion sur MySpace
  • Envoyer la discussion sur Yahoo