Mercurial > hg > octave-lyh

comparison scripts/statistics/distributions/beta_inv.m @ 4854:4b0f3b055331

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2004-04-07 02:37:05 by jwe]
author jwe
date Wed, 07 Apr 2004 02:38:06 +0000
parents 38c61cbf086c
children 4c8a2e4e0717
comparison
equal deleted inserted replaced
4853:66b3cce2bf37 4854:4b0f3b055331
31 31
32 if (nargin != 3) 32 if (nargin != 3)
33 usage ("beta_inv (x, a, b)"); 33 usage ("beta_inv (x, a, b)");
34 endif 34 endif
35 35
36 [retval, x, a, b] = common_size (x, a, b); 36 if (!isscalar (a) || !isscalar(b))
37 if (retval > 0) 37 [retval, x, a, b] = common_size (x, a, b);
38 error ("beta_inv: x, a and b must be of common size or scalars"); 38 if (retval > 0)
39 error ("beta_inv: x, a and b must be of common size or scalars");
40 endif
39 endif 41 endif
40 42
41 [r, c] = size (x); 43 sz = size (x);
42 s = r * c; 44 inv = zeros (sz);
43 x = reshape (x, s, 1);
44 a = reshape (a, s, 1);
45 b = reshape (b, s, 1);
46 inv = zeros (s, 1);
47 45
48 k = find ((x < 0) | (x > 1) | !(a > 0) | !(b > 0) | isnan (x)); 46 k = find ((x < 0) | (x > 1) | !(a > 0) | !(b > 0) | isnan (x));
49 if (any (k)) 47 if (any (k))
50 inv (k) = NaN * ones (length (k), 1); 48 inv (k) = NaN;
51 endif 49 endif
52 50
53 k = find ((x == 1) & (a > 0) & (b > 0)); 51 k = find ((x == 1) & (a > 0) & (b > 0));
54 if (any (k)) 52 if (any (k))
55 inv (k) = ones (length (k), 1); 53 inv (k) = 1;
56 endif 54 endif
57 55
58 k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0)); 56 k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0));
59 if (any (k)) 57 if (any (k))
60 a = a (k); 58 if (!isscalar(a) || !isscalar(b))
61 b = b (k); 59 a = a (k);
60 b = b (k);
61 y = a ./ (a + b);
62 else
63 y = a / (a + b) * ones (size (k));
64 endif
62 x = x (k); 65 x = x (k);
63 y = a ./ (a + b);
64 l = find (y < eps); 66 l = find (y < eps);
65 if (any (l)) 67 if (any (l))
66 y(l) = sqrt (eps) * ones (length (l), 1); 68 y(l) = sqrt (eps) * ones (length (l), 1);
67 endif 69 endif
68 l = find (y > 1 - eps); 70 l = find (y > 1 - eps);
90 endfor 92 endfor
91 93
92 inv (k) = y_new; 94 inv (k) = y_new;
93 endif 95 endif
94 96
95 inv = reshape (inv, r, c);
96
97 endfunction 97 endfunction