Mercurial > octave

comparison scripts/specfun/gammainc.m @ 28945:6e460773bdda

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use newlines after "function" and before "endfunction" for clarity. * quad2d.m, rng.m, __ok_cancel_dlg__.m, __unimplemented__.m, get_first_help_sentence.m, imformats.m, vectorize.m, ordeig.m, inputParser.m, tar_is_bsd.m, publish.m, axis.m, legend.m, polar.m, __plt__.m, print.m, __add_default_menu__.m, __gnuplot_draw_axes__.m, struct2hdl.m, movfun.m, gmres.m, pcg.m, __alltohandles__.m, tfqmr.m, betaincinv.m, gammainc.m, gammaincinv.m, mean.m, weboptions.m: Use newlines after "function" and before "endfunction" for clarity.
author Rik <rik@octave.org>
date Fri, 16 Oct 2020 09:27:56 -0700
parents d39b09b4c5db
children 7854d5752dd2
comparison
equal deleted inserted replaced
28944:d39b09b4c5db 28945:6e460773bdda
297 endif 297 endif
298 endfunction 298 endfunction
299 299
300 ## a == 1. 300 ## a == 1.
301 function y = gammainc_a1 (x, tail) 301 function y = gammainc_a1 (x, tail)
302
302 if (strcmp (tail, "lower")) 303 if (strcmp (tail, "lower"))
303 if (abs (x) < 1/2) 304 if (abs (x) < 1/2)
304 y = - expm1 (-x); 305 y = - expm1 (-x);
305 else 306 else
306 y = 1 - exp (-x); 307 y = 1 - exp (-x);
314 y = (exp (x) - 1) ./ x; 315 y = (exp (x) - 1) ./ x;
315 endif 316 endif
316 else 317 else
317 y = 1 ./ x; 318 y = 1 ./ x;
318 endif 319 endif
320
319 endfunction 321 endfunction
320 322
321 ## positive integer a; exp (x) and a! both under 1/eps 323 ## positive integer a; exp (x) and a! both under 1/eps
322 ## uses closed-form expressions for nonnegative integer a 324 ## uses closed-form expressions for nonnegative integer a
323 ## -- http://mathworld.wolfram.com/IncompleteGammaFunction.html. 325 ## -- http://mathworld.wolfram.com/IncompleteGammaFunction.html.
324 function y = gammainc_an (x, a, tail) 326 function y = gammainc_an (x, a, tail)
327
325 y = t = ones (size (x), class (x)); 328 y = t = ones (size (x), class (x));
326 i = 1; 329 i = 1;
327 while (any (a(:) > i)) 330 while (any (a(:) > i))
328 jj = (a > i); 331 jj = (a > i);
329 t(jj) .*= (x(jj) / i); 332 t(jj) .*= (x(jj) / i);
337 elseif (strcmp (tail, "scaledlower")) 340 elseif (strcmp (tail, "scaledlower"))
338 y = (1 - exp (-x) .* y) ./ D(x, a); 341 y = (1 - exp (-x) .* y) ./ D(x, a);
339 elseif (strcmp (tail, "scaledupper")) 342 elseif (strcmp (tail, "scaledupper"))
340 y .*= exp (-x) ./ D(x, a); 343 y .*= exp (-x) ./ D(x, a);
341 endif 344 endif
345
342 endfunction 346 endfunction
343 347
344 ## x + 0.25 < a | x < 0 | abs(x) < 1. 348 ## x + 0.25 < a | x < 0 | abs(x) < 1.
345 ## Numerical Recipes in Fortran 77 (6.2.5) 349 ## Numerical Recipes in Fortran 77 (6.2.5)
346 ## series 350 ## series
347 function y = gammainc_s (x, a, tail) 351 function y = gammainc_s (x, a, tail)
352
348 if (strcmp (tail, "scaledlower") || strcmp (tail, "scaledupper")) 353 if (strcmp (tail, "scaledlower") || strcmp (tail, "scaledupper"))
349 y = ones (size (x), class (x)); 354 y = ones (size (x), class (x));
350 term = x ./ (a + 1); 355 term = x ./ (a + 1);
351 else 356 else
352 ## Of course it is possible to scale at the end, but some tests fail. 357 ## Of course it is possible to scale at the end, but some tests fail.
365 if (strcmp (tail, "upper")) 370 if (strcmp (tail, "upper"))
366 y = 1 - y; 371 y = 1 - y;
367 elseif (strcmp (tail, "scaledupper")) 372 elseif (strcmp (tail, "scaledupper"))
368 y = 1 ./ D (x,a) - y; 373 y = 1 ./ D (x,a) - y;
369 endif 374 endif
375
370 endfunction 376 endfunction
371 377
372 ## x positive and large relative to a 378 ## x positive and large relative to a
373 ## NRF77 (6.2.7) 379 ## NRF77 (6.2.7)
374 ## Gamma (a,x)/Gamma (a) 380 ## Gamma (a,x)/Gamma (a)
375 ## Lentz's algorithm 381 ## Lentz's algorithm
376 ## __gammainc__ in libinterp/corefcn/__gammainc__.cc 382 ## __gammainc__ in libinterp/corefcn/__gammainc__.cc
377 function y = gammainc_l (x, a, tail) 383 function y = gammainc_l (x, a, tail)
384
378 y = __gammainc__ (x, a); 385 y = __gammainc__ (x, a);
379 if (strcmp (tail, "lower")) 386 if (strcmp (tail, "lower"))
380 y = 1 - y .* D (x, a); 387 y = 1 - y .* D (x, a);
381 elseif (strcmp (tail, "upper")) 388 elseif (strcmp (tail, "upper"))
382 y .*= D (x, a); 389 y .*= D (x, a);
383 elseif (strcmp (tail, "scaledlower")) 390 elseif (strcmp (tail, "scaledlower"))
384 y = 1 ./ D (x, a) - y; 391 y = 1 ./ D (x, a) - y;
385 endif 392 endif
393
386 endfunction 394 endfunction
387 395
388 ## Compute exp(-x)*x^a/Gamma(a+1) in a stable way for x and a large. 396 ## Compute exp(-x)*x^a/Gamma(a+1) in a stable way for x and a large.
389 ## 397 ##
390 ## L. Knusel, Computation of the Chi-square and Poisson distribution, 398 ## L. Knusel, Computation of the Chi-square and Poisson distribution,
391 ## SIAM J. Sci. Stat. Comput., 7(3), 1986 399 ## SIAM J. Sci. Stat. Comput., 7(3), 1986
392 ## which quotes Section 5, Abramowitz&Stegun 6.1.40, 6.1.41. 400 ## which quotes Section 5, Abramowitz&Stegun 6.1.40, 6.1.41.
393 function y = D (x, a) 401 function y = D (x, a)
402
394 athresh = 10; # FIXME: can this be better tuned? 403 athresh = 10; # FIXME: can this be better tuned?
395 y = zeros (size (x), class (x)); 404 y = zeros (size (x), class (x));
396 405
397 todo = true (size (x)); 406 todo = true (size (x));
398 todo(x == 0) = false; 407 todo(x == 0) = false;