nalgebra/base/statistics.rs
1use crate::allocator::Allocator;
2use crate::storage::RawStorage;
3use crate::{Const, DefaultAllocator, Dim, Matrix, OVector, RowOVector, Scalar, VectorView, U1};
4use num::{One, Zero};
5use simba::scalar::{ClosedAddAssign, ClosedMulAssign, Field, SupersetOf};
6use std::mem::MaybeUninit;
7
8/// # Folding on columns and rows
9impl<T: Scalar, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> {
10 /// Returns a row vector where each element is the result of the application of `f` on the
11 /// corresponding column of the original matrix.
12 #[inline]
13 #[must_use]
14 pub fn compress_rows(
15 &self,
16 f: impl Fn(VectorView<'_, T, R, S::RStride, S::CStride>) -> T,
17 ) -> RowOVector<T, C>
18 where
19 DefaultAllocator: Allocator<U1, C>,
20 {
21 let ncols = self.shape_generic().1;
22 let mut res = Matrix::uninit(Const::<1>, ncols);
23
24 for i in 0..ncols.value() {
25 // TODO: avoid bound checking of column.
26 // Safety: all indices are in range.
27 unsafe {
28 *res.get_unchecked_mut((0, i)) = MaybeUninit::new(f(self.column(i)));
29 }
30 }
31
32 // Safety: res is now fully initialized.
33 unsafe { res.assume_init() }
34 }
35
36 /// Returns a column vector where each element is the result of the application of `f` on the
37 /// corresponding column of the original matrix.
38 ///
39 /// This is the same as `self.compress_rows(f).transpose()`.
40 #[inline]
41 #[must_use]
42 pub fn compress_rows_tr(
43 &self,
44 f: impl Fn(VectorView<'_, T, R, S::RStride, S::CStride>) -> T,
45 ) -> OVector<T, C>
46 where
47 DefaultAllocator: Allocator<C>,
48 {
49 let ncols = self.shape_generic().1;
50 let mut res = Matrix::uninit(ncols, Const::<1>);
51
52 for i in 0..ncols.value() {
53 // TODO: avoid bound checking of column.
54 // Safety: all indices are in range.
55 unsafe {
56 *res.vget_unchecked_mut(i) = MaybeUninit::new(f(self.column(i)));
57 }
58 }
59
60 // Safety: res is now fully initialized.
61 unsafe { res.assume_init() }
62 }
63
64 /// Returns a column vector resulting from the folding of `f` on each column of this matrix.
65 #[inline]
66 #[must_use]
67 pub fn compress_columns(
68 &self,
69 init: OVector<T, R>,
70 f: impl Fn(&mut OVector<T, R>, VectorView<'_, T, R, S::RStride, S::CStride>),
71 ) -> OVector<T, R>
72 where
73 DefaultAllocator: Allocator<R>,
74 {
75 let mut res = init;
76
77 for i in 0..self.ncols() {
78 f(&mut res, self.column(i))
79 }
80
81 res
82 }
83}
84
85/// # Common statistics operations
86impl<T: Scalar, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> {
87 /*
88 *
89 * Sum computation.
90 *
91 */
92 /// The sum of all the elements of this matrix.
93 ///
94 /// # Example
95 ///
96 /// ```
97 /// # use nalgebra::Matrix2x3;
98 ///
99 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
100 /// 4.0, 5.0, 6.0);
101 /// assert_eq!(m.sum(), 21.0);
102 /// ```
103 #[inline]
104 #[must_use]
105 pub fn sum(&self) -> T
106 where
107 T: ClosedAddAssign + Zero,
108 {
109 self.iter().cloned().fold(T::zero(), |a, b| a + b)
110 }
111
112 /// The sum of all the rows of this matrix.
113 ///
114 /// Use `.row_sum_tr` if you need the result in a column vector instead.
115 ///
116 /// # Example
117 ///
118 /// ```
119 /// # use nalgebra::{Matrix2x3, Matrix3x2};
120 /// # use nalgebra::{RowVector2, RowVector3};
121 ///
122 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
123 /// 4.0, 5.0, 6.0);
124 /// assert_eq!(m.row_sum(), RowVector3::new(5.0, 7.0, 9.0));
125 ///
126 /// let mint = Matrix3x2::new(1, 2,
127 /// 3, 4,
128 /// 5, 6);
129 /// assert_eq!(mint.row_sum(), RowVector2::new(9,12));
130 /// ```
131 #[inline]
132 #[must_use]
133 pub fn row_sum(&self) -> RowOVector<T, C>
134 where
135 T: ClosedAddAssign + Zero,
136 DefaultAllocator: Allocator<U1, C>,
137 {
138 self.compress_rows(|col| col.sum())
139 }
140
141 /// The sum of all the rows of this matrix. The result is transposed and returned as a column vector.
142 ///
143 /// # Example
144 ///
145 /// ```
146 /// # use nalgebra::{Matrix2x3, Matrix3x2};
147 /// # use nalgebra::{Vector2, Vector3};
148 ///
149 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
150 /// 4.0, 5.0, 6.0);
151 /// assert_eq!(m.row_sum_tr(), Vector3::new(5.0, 7.0, 9.0));
152 ///
153 /// let mint = Matrix3x2::new(1, 2,
154 /// 3, 4,
155 /// 5, 6);
156 /// assert_eq!(mint.row_sum_tr(), Vector2::new(9, 12));
157 /// ```
158 #[inline]
159 #[must_use]
160 pub fn row_sum_tr(&self) -> OVector<T, C>
161 where
162 T: ClosedAddAssign + Zero,
163 DefaultAllocator: Allocator<C>,
164 {
165 self.compress_rows_tr(|col| col.sum())
166 }
167
168 /// The sum of all the columns of this matrix.
169 ///
170 /// # Example
171 ///
172 /// ```
173 /// # use nalgebra::{Matrix2x3, Matrix3x2};
174 /// # use nalgebra::{Vector2, Vector3};
175 ///
176 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
177 /// 4.0, 5.0, 6.0);
178 /// assert_eq!(m.column_sum(), Vector2::new(6.0, 15.0));
179 ///
180 /// let mint = Matrix3x2::new(1, 2,
181 /// 3, 4,
182 /// 5, 6);
183 /// assert_eq!(mint.column_sum(), Vector3::new(3, 7, 11));
184 /// ```
185 #[inline]
186 #[must_use]
187 pub fn column_sum(&self) -> OVector<T, R>
188 where
189 T: ClosedAddAssign + Zero,
190 DefaultAllocator: Allocator<R>,
191 {
192 let nrows = self.shape_generic().0;
193 self.compress_columns(OVector::zeros_generic(nrows, Const::<1>), |out, col| {
194 *out += col;
195 })
196 }
197
198 /*
199 *
200 * Product computation.
201 *
202 */
203 /// The product of all the elements of this matrix.
204 ///
205 /// # Example
206 ///
207 /// ```
208 /// # use nalgebra::Matrix2x3;
209 ///
210 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
211 /// 4.0, 5.0, 6.0);
212 /// assert_eq!(m.product(), 720.0);
213 /// ```
214 #[inline]
215 #[must_use]
216 pub fn product(&self) -> T
217 where
218 T: ClosedMulAssign + One,
219 {
220 self.iter().cloned().fold(T::one(), |a, b| a * b)
221 }
222
223 /// The product of all the rows of this matrix.
224 ///
225 /// Use `.row_sum_tr` if you need the result in a column vector instead.
226 ///
227 /// # Example
228 ///
229 /// ```
230 /// # use nalgebra::{Matrix2x3, Matrix3x2};
231 /// # use nalgebra::{RowVector2, RowVector3};
232 ///
233 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
234 /// 4.0, 5.0, 6.0);
235 /// assert_eq!(m.row_product(), RowVector3::new(4.0, 10.0, 18.0));
236 ///
237 /// let mint = Matrix3x2::new(1, 2,
238 /// 3, 4,
239 /// 5, 6);
240 /// assert_eq!(mint.row_product(), RowVector2::new(15, 48));
241 /// ```
242 #[inline]
243 #[must_use]
244 pub fn row_product(&self) -> RowOVector<T, C>
245 where
246 T: ClosedMulAssign + One,
247 DefaultAllocator: Allocator<U1, C>,
248 {
249 self.compress_rows(|col| col.product())
250 }
251
252 /// The product of all the rows of this matrix. The result is transposed and returned as a column vector.
253 ///
254 /// # Example
255 ///
256 /// ```
257 /// # use nalgebra::{Matrix2x3, Matrix3x2};
258 /// # use nalgebra::{Vector2, Vector3};
259 ///
260 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
261 /// 4.0, 5.0, 6.0);
262 /// assert_eq!(m.row_product_tr(), Vector3::new(4.0, 10.0, 18.0));
263 ///
264 /// let mint = Matrix3x2::new(1, 2,
265 /// 3, 4,
266 /// 5, 6);
267 /// assert_eq!(mint.row_product_tr(), Vector2::new(15, 48));
268 /// ```
269 #[inline]
270 #[must_use]
271 pub fn row_product_tr(&self) -> OVector<T, C>
272 where
273 T: ClosedMulAssign + One,
274 DefaultAllocator: Allocator<C>,
275 {
276 self.compress_rows_tr(|col| col.product())
277 }
278
279 /// The product of all the columns of this matrix.
280 ///
281 /// # Example
282 ///
283 /// ```
284 /// # use nalgebra::{Matrix2x3, Matrix3x2};
285 /// # use nalgebra::{Vector2, Vector3};
286 ///
287 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
288 /// 4.0, 5.0, 6.0);
289 /// assert_eq!(m.column_product(), Vector2::new(6.0, 120.0));
290 ///
291 /// let mint = Matrix3x2::new(1, 2,
292 /// 3, 4,
293 /// 5, 6);
294 /// assert_eq!(mint.column_product(), Vector3::new(2, 12, 30));
295 /// ```
296 #[inline]
297 #[must_use]
298 pub fn column_product(&self) -> OVector<T, R>
299 where
300 T: ClosedMulAssign + One,
301 DefaultAllocator: Allocator<R>,
302 {
303 let nrows = self.shape_generic().0;
304 self.compress_columns(
305 OVector::repeat_generic(nrows, Const::<1>, T::one()),
306 |out, col| {
307 out.component_mul_assign(&col);
308 },
309 )
310 }
311
312 /*
313 *
314 * Variance computation.
315 *
316 */
317 /// The variance of all the elements of this matrix.
318 ///
319 /// # Example
320 ///
321 /// ```
322 /// # #[macro_use] extern crate approx;
323 /// # use nalgebra::Matrix2x3;
324 ///
325 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
326 /// 4.0, 5.0, 6.0);
327 /// assert_relative_eq!(m.variance(), 35.0 / 12.0, epsilon = 1.0e-8);
328 /// ```
329 #[inline]
330 #[must_use]
331 pub fn variance(&self) -> T
332 where
333 T: Field + SupersetOf<f64>,
334 {
335 if self.is_empty() {
336 T::zero()
337 } else {
338 let n_elements: T = crate::convert(self.len() as f64);
339 let mean = self.mean();
340
341 self.iter().cloned().fold(T::zero(), |acc, x| {
342 acc + (x.clone() - mean.clone()) * (x - mean.clone())
343 }) / n_elements
344 }
345 }
346
347 /// The variance of all the rows of this matrix.
348 ///
349 /// Use `.row_variance_tr` if you need the result in a column vector instead.
350 /// # Example
351 ///
352 /// ```
353 /// # use nalgebra::{Matrix2x3, RowVector3};
354 ///
355 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
356 /// 4.0, 5.0, 6.0);
357 /// assert_eq!(m.row_variance(), RowVector3::new(2.25, 2.25, 2.25));
358 /// ```
359 #[inline]
360 #[must_use]
361 pub fn row_variance(&self) -> RowOVector<T, C>
362 where
363 T: Field + SupersetOf<f64>,
364 DefaultAllocator: Allocator<U1, C>,
365 {
366 self.compress_rows(|col| col.variance())
367 }
368
369 /// The variance of all the rows of this matrix. The result is transposed and returned as a column vector.
370 ///
371 /// # Example
372 ///
373 /// ```
374 /// # use nalgebra::{Matrix2x3, Vector3};
375 ///
376 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
377 /// 4.0, 5.0, 6.0);
378 /// assert_eq!(m.row_variance_tr(), Vector3::new(2.25, 2.25, 2.25));
379 /// ```
380 #[inline]
381 #[must_use]
382 pub fn row_variance_tr(&self) -> OVector<T, C>
383 where
384 T: Field + SupersetOf<f64>,
385 DefaultAllocator: Allocator<C>,
386 {
387 self.compress_rows_tr(|col| col.variance())
388 }
389
390 /// The variance of all the columns of this matrix.
391 ///
392 /// # Example
393 ///
394 /// ```
395 /// # #[macro_use] extern crate approx;
396 /// # use nalgebra::{Matrix2x3, Vector2};
397 ///
398 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
399 /// 4.0, 5.0, 6.0);
400 /// assert_relative_eq!(m.column_variance(), Vector2::new(2.0 / 3.0, 2.0 / 3.0), epsilon = 1.0e-8);
401 /// ```
402 #[inline]
403 #[must_use]
404 pub fn column_variance(&self) -> OVector<T, R>
405 where
406 T: Field + SupersetOf<f64>,
407 DefaultAllocator: Allocator<R>,
408 {
409 let (nrows, ncols) = self.shape_generic();
410
411 let mut mean = self.column_mean();
412 mean.apply(|e| *e = -(e.clone() * e.clone()));
413
414 let denom = T::one() / crate::convert::<_, T>(ncols.value() as f64);
415 self.compress_columns(mean, |out, col| {
416 for i in 0..nrows.value() {
417 unsafe {
418 let val = col.vget_unchecked(i);
419 *out.vget_unchecked_mut(i) += denom.clone() * val.clone() * val.clone()
420 }
421 }
422 })
423 }
424
425 /*
426 *
427 * Mean computation.
428 *
429 */
430 /// The mean of all the elements of this matrix.
431 ///
432 /// # Example
433 ///
434 /// ```
435 /// # use nalgebra::Matrix2x3;
436 ///
437 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
438 /// 4.0, 5.0, 6.0);
439 /// assert_eq!(m.mean(), 3.5);
440 /// ```
441 #[inline]
442 #[must_use]
443 pub fn mean(&self) -> T
444 where
445 T: Field + SupersetOf<f64>,
446 {
447 if self.is_empty() {
448 T::zero()
449 } else {
450 self.sum() / crate::convert(self.len() as f64)
451 }
452 }
453
454 /// The mean of all the rows of this matrix.
455 ///
456 /// Use `.row_mean_tr` if you need the result in a column vector instead.
457 ///
458 /// # Example
459 ///
460 /// ```
461 /// # use nalgebra::{Matrix2x3, RowVector3};
462 ///
463 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
464 /// 4.0, 5.0, 6.0);
465 /// assert_eq!(m.row_mean(), RowVector3::new(2.5, 3.5, 4.5));
466 /// ```
467 #[inline]
468 #[must_use]
469 pub fn row_mean(&self) -> RowOVector<T, C>
470 where
471 T: Field + SupersetOf<f64>,
472 DefaultAllocator: Allocator<U1, C>,
473 {
474 self.compress_rows(|col| col.mean())
475 }
476
477 /// The mean of all the rows of this matrix. The result is transposed and returned as a column vector.
478 ///
479 /// # Example
480 ///
481 /// ```
482 /// # use nalgebra::{Matrix2x3, Vector3};
483 ///
484 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
485 /// 4.0, 5.0, 6.0);
486 /// assert_eq!(m.row_mean_tr(), Vector3::new(2.5, 3.5, 4.5));
487 /// ```
488 #[inline]
489 #[must_use]
490 pub fn row_mean_tr(&self) -> OVector<T, C>
491 where
492 T: Field + SupersetOf<f64>,
493 DefaultAllocator: Allocator<C>,
494 {
495 self.compress_rows_tr(|col| col.mean())
496 }
497
498 /// The mean of all the columns of this matrix.
499 ///
500 /// # Example
501 ///
502 /// ```
503 /// # use nalgebra::{Matrix2x3, Vector2};
504 ///
505 /// let m = Matrix2x3::new(1.0, 2.0, 3.0,
506 /// 4.0, 5.0, 6.0);
507 /// assert_eq!(m.column_mean(), Vector2::new(2.0, 5.0));
508 /// ```
509 #[inline]
510 #[must_use]
511 pub fn column_mean(&self) -> OVector<T, R>
512 where
513 T: Field + SupersetOf<f64>,
514 DefaultAllocator: Allocator<R>,
515 {
516 let (nrows, ncols) = self.shape_generic();
517 let denom = T::one() / crate::convert::<_, T>(ncols.value() as f64);
518 self.compress_columns(OVector::zeros_generic(nrows, Const::<1>), |out, col| {
519 out.axpy(denom.clone(), &col, T::one())
520 })
521 }
522}