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:

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.

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

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:

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.

Citation:

Envoyé par yogitetradim

...En fait ce qui m'a induit en erreur c'est qu'avec la ligne 132:

Code:

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)).