Mercurial > hg > octave-nkf

comparison scripts/specfun/erfinv.m @ 3266:3add30d33bd0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1999-09-03 17:15:40 by jwe]
author jwe
date Fri, 03 Sep 1999 17:15:41 +0000
parents e330cb788508
children 6923abb04e16
comparison
equal deleted inserted replaced
3265:a2b3a1413d28 3266:3add30d33bd0
1 ## Copyright (C) 1995, 1996, 1997 Kurt Hornik 1 ## Copyright (C) 1995, 1996 Kurt Hornik
2 ## 2 ##
3 ## This program is free software; you can redistribute it and/or modify 3 ## This program is free software; you can redistribute it and/or modify
4 ## it under the terms of the GNU General Public License as published by 4 ## it under the terms of the GNU General Public License as published by
5 ## the Free Software Foundation; either version 2, or (at your option) 5 ## the Free Software Foundation; either version 2, or (at your option)
6 ## any later version. 6 ## any later version.
20 20
21 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> 21 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
22 ## Created: 27 September 1994 22 ## Created: 27 September 1994
23 ## Adapted-By: jwe 23 ## Adapted-By: jwe
24 24
25 function y = erfinv (x) 25 function [y, iterations] = erfinv (x)
26 26
27 if (nargin != 1) 27 if (nargin != 1)
28 usage ("erfinv (x)"); 28 usage ("erfinv (x)");
29 endif 29 endif
30 30
31 maxit = 100;
32 tol = eps;
33
34 iterations = 0;
35
31 [m, n] = size (x); 36 [m, n] = size (x);
32 x = reshape (x, m * n, 1); 37 x = reshape (x, m * n, 1);
33 y = zeros (m * n, 1); 38 y = zeros (m * n, 1);
34 39
35 i = find ((x < -1) | (x > 1)); 40 i = find ((x < -1) | (x > 1));
36 if any (i) 41 if any (i)
37 y(i) = NaN * ones (length (i), 1); 42 y(i) = NaN * ones (length (i), 1);
38 endif 43 endif
39 44
40 t = find (x == -1); 45 t = find (x == -1);
41 y (t) = (-Inf) * ones (size (t)); 46 y (t) = (-Inf) * ones (size (t));
42 47
43 t = find (x == 1); 48 t = find (x == 1);
44 y (t) = Inf * ones (size (t)); 49 y (t) = Inf * ones (size (t));
45 50
46 i = find ((x > -1) & (x < 1)); 51 i = find ((x > -1) & (x < 1));
47 if any (i) 52 if any (i)
48 s = sqrt (pi) / 2; 53 s = sqrt (pi) / 2;
49 z_old = ones (length (i), 1); 54 z_old = ones (length (i), 1);
50 z_new = zeros (length (i), 1); 55 z_new = sqrt (-log (1 - abs (x(i)))) .* sign (x(i));
51 while (any (any (abs (z_new - z_old) > 2 * eps))) 56 while (any (abs (erf (z_old) - x(i)) > tol * abs (x(i))))
52 z_old = z_new; 57 z_old = z_new;
53 z_new = z_old - (erf (z_old) - x(i)) .* exp (z_old.^2) * s; 58 z_new = z_old - (erf (z_old) - x(i)) .* exp (z_old.^2) * s;
59 if (++iterations > maxit)
60 warning ("erfinv: iteration limit exceeded");
61 break;
62 endif
54 endwhile 63 endwhile
55 y(i) = z_new; 64 y(i) = z_new;
56 endif 65 endif
57 66
58 y = reshape (y, m, n); 67 y = reshape (y, m, n);