Mercurial > hg > octave-nkf

comparison scripts/plot/hist.m @ 4880:b9662e2ceb6b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2004-04-23 16:13:49 by jwe]
author jwe
date Fri, 23 Apr 2004 16:13:49 +0000
parents 0358ed4394f5
children c08cb1098afc
comparison
equal deleted inserted replaced
4879:013350fee837 4880:b9662e2ceb6b
47 47
48 if (nargin < 1 || nargin > 3) 48 if (nargin < 1 || nargin > 3)
49 usage ("[nn, xx] = hist (y, x, norm)"); 49 usage ("[nn, xx] = hist (y, x, norm)");
50 endif 50 endif
51 51
52 if (isvector (y)) 52 transpose = rows (y) == 1;
53 if (transpose)
54 y = y(:);
55 endif
56
57 if (isreal (y))
53 max_val = max (y); 58 max_val = max (y);
54 min_val = min (y); 59 min_val = min (y);
55 else 60 else
56 error ("hist: first argument must be a vector"); 61 error ("hist: first argument must be a vector");
57 endif 62 endif
58 63
59 if (nargin == 1) 64 if (nargin == 1)
60 n = 10; 65 n = 10;
61 delta = (max_val - min_val) / (n-1) / 2; 66 x = [0.5:n]'/n;
62 cutoff = linspace (min_val+delta, max_val-delta, n-1); 67 x = x * (max_val - min_val) + ones(size(x)) * min_val;
63 else 68 else
64 ## nargin is either 2 or 3 69 ## nargin is either 2 or 3
65 if (isscalar (x)) 70 if (isscalar (x))
66 n = x; 71 n = x;
67 if (n <= 0) 72 if (n <= 0)
68 error ("hist: number of bins must be positive"); 73 error ("hist: number of bins must be positive");
69 endif 74 endif
70 delta = (max_val - min_val) / (n-1) / 2; 75 x = [0.5:n]'/n;
71 cutoff = linspace (min_val+delta, max_val-delta, n-1); 76 x = x * (max_val - min_val) + ones(size(x)) * min_val;
72 elseif (isvector (x)) 77 elseif (isreal (x))
78 if (isvector (x))
79 x = x(:);
80 endif
73 tmp = sort (x); 81 tmp = sort (x);
74 if (any (tmp != x)) 82 if (any (tmp != x))
75 warning ("hist: bin values not sorted on input"); 83 warning ("hist: bin values not sorted on input");
76 x = tmp; 84 x = tmp;
77 endif 85 endif
78 cutoff = (x(1:end-1) + x(2:end)) / 2;
79 n = length (x);
80 else 86 else
81 error ("hist: second argument must be a scalar or a vector"); 87 error ("hist: second argument must be a scalar or a vector");
82 endif 88 endif
83 endif 89 endif
84 90
85 if (n < 30) 91 cutoff = (x(1:end-1,:) + x(2:end,:)) / 2;
92 n = rows (x);
93 if (n < 30 && columns (x) == 1)
86 ## The following algorithm works fastest for n less than about 30. 94 ## The following algorithm works fastest for n less than about 30.
87 chist = [zeros(n,1); length(y)]; 95 chist = zeros (n+1, columns (y));
88 for i = 1:n-1 96 for i = 1:n-1
89 chist(i+1) = sum (y < cutoff(i)); 97 chist(i+1,:) = sum (y <= cutoff(i));
90 endfor 98 endfor
99 chist(n+1,:) = rows (y);
91 else 100 else
92 ## The following algorithm works fastest for n greater than about 30. 101 ## The following algorithm works fastest for n greater than about 30.
93 ## Put cutoff elements between boundaries, integrate over all 102 ## Put cutoff elements between boundaries, integrate over all
94 ## elements, keep totals at boundaries. 103 ## elements, keep totals at boundaries.
95 [s, idx] = sort ([cutoff(:); y(:)]); 104 [s, idx] = sort ([y; cutoff]);
96 chist = cumsum(idx>=n); 105 len = rows (y);
97 chist = [0; chist(idx<n); chist(end)]; 106 chist = cumsum (idx <= len);
107 t1 = zeros (1, columns (y));
108 t2 = reshape (chist(idx > len), size (cutoff));
109 t3 = chist(end,:);
110 chist = [t1; t2; t3];
98 endif 111 endif
99 112
100 freq= diff(chist)'; 113 freq = diff (chist);
101 114
102 if (nargin == 3) 115 if (nargin == 3)
103 ## Normalise the histogram. 116 ## Normalise the histogram.
104 freq = freq / length (y) * norm; 117 freq = freq / rows (y) * norm;
105 endif 118 endif
106 119
107 if (nargout > 0) 120 if (nargout > 0)
108 nn = freq; 121 if (transpose)
109 xx = x; 122 nn = freq';
123 xx = x';
124 else
125 nn = freq;
126 xx = x;
127 endif
110 else 128 else
111 bar (x, freq); 129 bar (x, freq);
112 endif 130 endif
113 131
114 endfunction 132 endfunction
115 133
116 %!test 134 %!test
117 %! for n = [1, 10, 30, 100, 1000] 135 %! [nn,xx]=hist([1:4],3);
118 %! assert( sum(hist([1:n], [1:n])) == n ); 136 %! assert(xx, [1.5,2.5,3.5]);
137 %! assert(nn, [2,1,1]);
138 %!test
139 %! [nn,xx]=hist([1:4]',3);
140 %! assert(xx, [1.5,2.5,3.5]');
141 %! assert(nn, [2,1,1]');
142 %!test
143 %! [nn,xx]=hist([[1:4]',[1:4]'],3);
144 %! assert(xx, [[1.5,2.5,3.5]',[1.5,2.5,3.5]']);
145 %! assert(nn, [[2,1,1]',[2,1,1]']);
146 %!assert(hist(1,1),1);
147 %!test
148 %! for n = [10, 30, 100, 1000]
149 %! assert( sum(hist([1:n], n)), n );
150 %! assert( sum(hist([1:n], [2:n-1])), n);
151 %! assert( sum(hist([1:n], [1:n])), n );
152 %! assert( sum(hist([1:n], 29)), n);
153 %! assert( sum(hist([1:n], 30)), n);
119 %! endfor 154 %! endfor