Merge branch 'main' into dataframe

2025-10-04 07:49:26 +00:00 · 2025-08-24 21:29:31 +01:00
parent 39a95e63d9 c53693fa7b
commit 811c153eaf
59 changed files with 2838 additions and 205 deletions
--- a/.github/htmldocs/index.html
+++ b/.github/htmldocs/index.html
@@ -58,6 +58,14 @@
        <h2>A lightweight dataframe & math toolkit for Rust</h2>
        <hr style="border: 1px solid #d4d4d4; margin: 20px 0;">
        <p>
+
+            🐙 <a href="https://github.com/Magnus167/rustframe">GitHub</a>
+            <br><br>
+            
+            📖 <a href="https://magnus167.github.io/rustframe/user-guide">User Guide</a>
+            <br><br>
+            
+
            📚 <a href="https://magnus167.github.io/rustframe/docs">Docs</a> |
            📊 <a href="https://magnus167.github.io/rustframe/benchmark-report/">Benchmarks</a>

@@ -65,8 +73,7 @@
            🦀 <a href="https://crates.io/crates/rustframe">Crates.io</a> |
            🔖 <a href="https://docs.rs/rustframe/latest/rustframe/">docs.rs</a>
            <br><br>
-            🐙 <a href="https://github.com/Magnus167/rustframe">GitHub</a> |
-            🌐 <a href="https://gitea.nulltech.uk/Magnus167/rustframe">Gitea mirror</a>
+            <!-- 🌐 <a href="https://gitea.nulltech.uk/Magnus167/rustframe">Gitea mirror</a> -->
        </p>
    </main>
 </body>
--- a/.github/runners/runner-x64/Dockerfile
+++ b/.github/runners/runner-x64/Dockerfile
@@ -7,7 +7,7 @@ ARG DEBIAN_FRONTEND=noninteractive

 RUN apt update -y && apt upgrade -y && useradd -m docker
 RUN apt install -y --no-install-recommends \
-    curl jq git unzip \
+    curl jq git zip unzip \
    # dev dependencies
    build-essential libssl-dev libffi-dev python3 python3-venv python3-dev python3-pip \
    # dot net core dependencies
--- a/.github/scripts/run_examples.sh
+++ b/.github/scripts/run_examples.sh
@@ -0,0 +1,16 @@
+cargo build --release --examples 
+
+for ex in examples/*.rs; do
+  name=$(basename "$ex" .rs)
+  echo
+  echo "🟡 Running example: $name"
+
+  if ! cargo run --release --example "$name" -- --debug; then
+    echo
+    echo "❌ Example '$name' failed. Aborting."
+    exit 1
+  fi
+done
+
+echo
+echo "✅ All examples ran successfully."
--- a/.github/workflows/docs-and-testcov.yml
+++ b/.github/workflows/docs-and-testcov.yml
@@ -151,7 +151,8 @@ jobs:
          mkdir -p target/doc/docs
          mv target/doc/rustframe/* target/doc/docs/

-          mkdir output
+          echo "<meta http-equiv=\"refresh\" content=\"0; url=../docs/index.html\">" > target/doc/rustframe/index.html
+
          cp tarpaulin-report.html target/doc/docs/
          cp tarpaulin-report.json target/doc/docs/
          cp tarpaulin-badge.json target/doc/docs/
@@ -164,16 +165,30 @@ jobs:
          # copy the benchmark report to the output directory
          cp -r benchmark-report target/doc/

+          mkdir output
+          cp -r target/doc/* output/
+
+      - name: Build user guide
+        run: |
+          cargo binstall mdbook
+          bash ./docs/build.sh
+
+      - name: Copy user guide to output directory
+        run: |
+          mkdir output/user-guide
+          cp -r docs/book/* output/user-guide/
+
      - name: Add index.html to output directory
        run: |
-          cp .github/htmldocs/index.html target/doc/index.html
-          cp .github/rustframe_logo.png target/doc/rustframe_logo.png
+          cp .github/htmldocs/index.html output/index.html
+          cp .github/rustframe_logo.png output/rustframe_logo.png

      - name: Upload Pages artifact
        # if: github.event_name == 'push' || github.event_name == 'workflow_dispatch'
        uses: actions/upload-pages-artifact@v3
        with:
-          path: target/doc/
+          # path: target/doc/
+          path: output/

      - name: Deploy to GitHub Pages
        # if: github.event_name == 'push' || github.event_name == 'workflow_dispatch'
--- a/.github/workflows/run-unit-tests.yml
+++ b/.github/workflows/run-unit-tests.yml
@@ -12,14 +12,12 @@ concurrency:

 jobs:
  pick-runner:
-
    if: github.event.pull_request.draft == false

    runs-on: ubuntu-latest
    outputs:
      runner: ${{ steps.choose.outputs.use-runner }}
    steps:
-
      - uses: actions/checkout@v4
      - id: choose
        uses: ./.github/actions/runner-fallback
@@ -27,7 +25,6 @@ jobs:
          primary-runner: "self-hosted"
          fallback-runner: "ubuntu-latest"
          github-token: ${{ secrets.CUSTOM_GH_TOKEN }}
-          

  run-unit-tests:
    needs: pick-runner
@@ -56,6 +53,20 @@ jobs:
      - name: Test docs generation
        run: cargo doc --no-deps --release

+      - name: Test examples
+        run: cargo test --examples --release
+
+      - name: Run all examples
+        run: |
+          for example in examples/*.rs; do
+            name=$(basename "$example" .rs)
+            echo "Running example: $name"
+            cargo run --release --example "$name" -- --debug || exit 1
+          done
+
+      - name: Cargo test all targets
+        run: cargo test --all-targets --release
+
      - name: Upload coverage to Codecov
        uses: codecov/codecov-action@v3
        with:
@@ -67,3 +78,8 @@ jobs:
        uses: codecov/test-results-action@v1
        with:
          token: ${{ secrets.CODECOV_TOKEN }}
+
+      - name: Test build user guide
+        run: |
+          cargo binstall mdbook
+          bash ./docs/build.sh
--- a/.gitignore
+++ b/.gitignore
@@ -16,4 +16,6 @@ data/

 tarpaulin-report.*

-.github/htmldocs/rustframe_logo.png
+.github/htmldocs/rustframe_logo.png
+
+docs/book/
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,10 +1,12 @@
 [package]
 name = "rustframe"
-version = "0.0.1-a.0"
+authors = ["Palash Tyagi (https://github.com/Magnus167)"]
+version = "0.0.1-a.20250805"
 edition = "2021"
 license = "GPL-3.0-or-later"
 readme = "README.md"
-description = "A simple dataframe library"
+description = "A simple dataframe and math toolkit"
+documentation = "https://magnus167.github.io/rustframe/"

 [lib]
 name = "rustframe"
@@ -14,14 +16,10 @@ crate-type = ["cdylib", "lib"]
 [dependencies]
 chrono = "^0.4.10"
 criterion = { version = "0.5", features = ["html_reports"], optional = true }
-rand = "^0.9.1"

 [features]
 bench = ["dep:criterion"]

-# [dev-dependencies]
-# criterion = { version = "0.5", features = ["html_reports"], optional = true }
-
 [[bench]]
 name = "benchmarks"
 harness = false
--- a/README.md
+++ b/README.md
@@ -1,40 +1,63 @@
 # rustframe

-<!-- # <img align="center" alt="Rustframe" src=".github/rustframe_logo.png" height="50px" /> rustframe -->
-
-<!-- though the centre tag doesn't work as it would noramlly, it achieves the desired effect -->
-
-📚 [Docs](https://magnus167.github.io/rustframe/) | 🐙 [GitHub](https://github.com/Magnus167/rustframe) | 🌐 [Gitea mirror](https://gitea.nulltech.uk/Magnus167/rustframe) | 🦀 [Crates.io](https://crates.io/crates/rustframe) | 🔖 [docs.rs](https://docs.rs/rustframe/latest/rustframe/)
+🐙 [GitHub](https://github.com/Magnus167/rustframe) | 📚 [Docs](https://magnus167.github.io/rustframe/) | 📖 [User Guide](https://magnus167.github.io/rustframe/user-guide/) | 🦀 [Crates.io](https://crates.io/crates/rustframe) | 🔖 [docs.rs](https://docs.rs/rustframe/latest/rustframe/)

 <!-- [![Last commit](https://img.shields.io/endpoint?url=https://magnus167.github.io/rustframe/rustframe/last-commit-date.json)](https://github.com/Magnus167/rustframe) -->

 [![codecov](https://codecov.io/gh/Magnus167/rustframe/graph/badge.svg?token=J7ULJEFTVI)](https://codecov.io/gh/Magnus167/rustframe)
 [![Coverage](https://img.shields.io/endpoint?url=https://magnus167.github.io/rustframe/docs/tarpaulin-badge.json)](https://magnus167.github.io/rustframe/docs/tarpaulin-report.html)
+[![gitea-mirror](https://img.shields.io/badge/git_mirror-blue)](https://gitea.nulltech.uk/Magnus167/rustframe)

 ---

 ## Rustframe: _A lightweight dataframe & math toolkit for Rust_

-Rustframe provides intuitive dataframe, matrix, and series operations small-to-mid scale data analysis and manipulation.
+Rustframe provides intuitive dataframe, matrix, and series operations for data analysis and manipulation.

-Rustframe keeps things simple, safe, and readable. It is handy for quick numeric experiments and small analytical tasks, but it is **not** meant to compete with powerhouse crates like `polars` or `ndarray`.
+Rustframe keeps things simple, safe, and readable. It is handy for quick numeric experiments and small analytical tasks as well as for educational purposes. It is designed to be easy to use and understand, with a clean API implemented in 100% safe Rust.
+
+Rustframe is an educational project, and is not intended for production use. It is **not** meant to compete with powerhouse crates like `polars` or `ndarray`. It is a work in progress, and the API is subject to change. There are no guarantees of stability or performance, and it is not optimized for large datasets or high-performance computing.

 ### What it offers

- **Math that reads like math** - element‑wise `+`, `−`, `×`, `÷` on entire frames or scalars.
- **Broadcast & reduce** - sum, product, any/all across rows or columns without boilerplate.
- **Boolean masks made simple** - chain comparisons, combine with `&`/`|`, get a tidy `BoolMatrix` back.
- **Date‑centric row index** - business‑day ranges and calendar slicing built in.
- **Pure safe Rust** - 100 % safe, zero `unsafe`.
+- **Matrix operations** - Element-wise arithmetic, boolean logic, transpose, and more.
+- **Math that reads like math** - element-wise `+`, `−`, `×`, `÷` on entire frames or scalars.
+- **Frames** - Column major data structure for single-type data, with labeled columns and typed row indices.
+- **Compute module** - Implements various statistical computations and machine learning models.
+- **Random number utils** - Built-in pseudo and cryptographically secure generators for simulations.
+- **[Coming Soon]** _DataFrame_ - Multi-type data structure for heterogeneous data, with labeled columns and typed row indices.
+
+#### Matrix and Frame functionality
+
+- **Matrix operations** - Element-wise arithmetic, boolean logic, transpose, and more.
+- **Frame operations** - Column manipulation, sorting, and more.
+
+#### Compute Module
+
+The `compute` module provides implementations for various statistical computations and machine learning models.
+
+**Statistics, Data Analysis, and Machine Learning:**
+
+- Correlation analysis
+- Descriptive statistics
+- Distributions
+- Inferential statistics
+
+- Dense Neural Networks
+- Gaussian Naive Bayes
+- K-Means Clustering
+- Linear Regression
+- Logistic Regression
+- Principal Component Analysis

 ### Heads up

 - **Not memory‑efficient (yet)** - footprint needs work.
- **Feature set still small** - expect missing pieces.
+- **The feature set is still limited** - expect missing pieces.

-### On the horizon
+### Somewhere down the line

- Optional GPU help (Vulkan or similar) for heavier workloads.
+- Optional GPU acceleration (Vulkan or similar) for heavier workloads.
 - Straightforward Python bindings using `pyo3`.

 ---
@@ -51,7 +74,7 @@ use rustframe::{

 let n_periods = 4;

-// Four business days starting 2024‑01‑02
+// Four business days starting 2024-01-02
 let dates: Vec<NaiveDate> =
    BDatesList::from_n_periods("2024-01-02".to_string(), DateFreq::Daily, n_periods)
        .unwrap()
@@ -86,13 +109,13 @@ let result: Matrix<f64> = result / 2.0; // divide by scalar
 let check: bool = result.eq_elem(ma.clone()).all();
 assert!(check);

-// The above math can also be written as:
+// Alternatively:
 let check: bool = (&(&(&(&ma + 1.0) - 1.0) * 2.0) / 2.0)
    .eq_elem(ma.clone())
    .all();
 assert!(check);

-// The above math can also be written as:
+// or even as:
 let check: bool = ((((ma.clone() + 1.0) - 1.0) * 2.0) / 2.0)
    .eq_elem(ma.clone())
    .all();
@@ -103,10 +126,6 @@ let mc: Matrix<f64> = Matrix::from_cols(vec![vec![1.0, 2.0], vec![3.0, 4.0]]);
 let md: Matrix<f64> = Matrix::from_cols(vec![vec![5.0, 6.0], vec![7.0, 8.0]]);
 let mul_result: Matrix<f64> = mc.matrix_mul(&md);
 // Expected:
-// 1*5 + 3*6 = 5 + 18 = 23
-// 2*5 + 4*6 = 10 + 24 = 34
-// 1*7 + 3*8 = 7 + 24 = 31
-// 2*7 + 4*8 = 14 + 32 = 46
 assert_eq!(mul_result.data(), &[23.0, 34.0, 31.0, 46.0]);

 // Dot product (alias for matrix_mul for FloatMatrix)
@@ -115,14 +134,7 @@ assert_eq!(dot_result, mul_result);

 // Transpose
 let original_matrix: Matrix<f64> = Matrix::from_cols(vec![vec![1.0, 2.0, 3.0], vec![4.0, 5.0, 6.0]]);
-// Original:
-// 1 4
-// 2 5
-// 3 6
 let transposed_matrix: Matrix<f64> = original_matrix.transpose();
-// Transposed:
-// 1 2 3
-// 4 5 6
 assert_eq!(transposed_matrix.rows(), 2);
 assert_eq!(transposed_matrix.cols(), 3);
 assert_eq!(transposed_matrix.data(), &[1.0, 4.0, 2.0, 5.0, 3.0, 6.0]);
@@ -131,10 +143,6 @@ assert_eq!(transposed_matrix.data(), &[1.0, 4.0, 2.0, 5.0, 3.0, 6.0]);
 let matrix = Matrix::from_cols(vec![vec![1.0, 2.0, 3.0], vec![4.0, 5.0, 6.0]]);
 // Map function to double each value
 let mapped_matrix = matrix.map(|x| x * 2.0);
-// Expected data after mapping
-// 2 8
-// 4 10
-// 6 12
 assert_eq!(mapped_matrix.data(), &[2.0, 4.0, 6.0, 8.0, 10.0, 12.0]);

 // Zip
@@ -142,9 +150,6 @@ let a = Matrix::from_cols(vec![vec![1.0, 2.0], vec![3.0, 4.0]]); // 2x2 matrix
 let b = Matrix::from_cols(vec![vec![5.0, 6.0], vec![7.0, 8.0]]); // 2x2 matrix
                                                                   // Zip function to add corresponding elements
 let zipped_matrix = a.zip(&b, |x, y| x + y);
-// Expected data after zipping
-// 6 10
-// 8 12
 assert_eq!(zipped_matrix.data(), &[6.0, 8.0, 10.0, 12.0]);
 ```

@@ -275,7 +280,7 @@ println!(
 // df.drop_column("non_existent_col"); // Uncomment to see panic
 ```

-### More examples
+## More examples

 See the [examples](./examples/) directory for some demonstrations of Rustframe's syntax and functionality.

@@ -290,3 +295,45 @@ E.g. to run the `game_of_life` example:
 ```bash
 cargo run --example game_of_life
 ```
+
+More demos:
+
+```bash
+cargo run --example linear_regression
+cargo run --example logistic_regression
+cargo run --example k_means
+cargo run --example pca
+cargo run --example stats_overview
+cargo run --example descriptive_stats
+cargo run --example correlation
+cargo run --example inferential_stats
+cargo run --example distributions
+```
+
+To simply list all available examples, you can run:
+
+```bash
+# this technically raises an error, but it will list all examples
+cargo run --example
+```
+
+Each demo runs a couple of mini-scenarios showcasing the APIs.
+
+## Running benchmarks
+
+To run the benchmarks, use:
+
+```bash
+cargo bench --features "bench"
+```
+
+## Building the user-guide
+
+To build the user guide, use:
+
+```bash
+cargo binstall mdbook
+bash docs/build.sh
+```
+
+This will generate the user guide in the `docs/book` directory.
--- a/docs/book.toml
+++ b/docs/book.toml
@@ -0,0 +1,7 @@
+[book]
+title = "Rustframe User Guide"
+authors = ["Palash Tyagi (https://github.com/Magnus167)"]
+description = "Guided journey through Rustframe capabilities."
+
+[build]
+build-dir = "book"
--- a/docs/build.sh
+++ b/docs/build.sh
@@ -0,0 +1,7 @@
+#!/usr/bin/env sh
+# Build and test the Rustframe user guide using mdBook.
+set -e
+
+cd docs
+bash gen.sh "$@"
+cd ..
--- a/docs/gen.sh
+++ b/docs/gen.sh
@@ -0,0 +1,14 @@
+#!/usr/bin/env sh
+
+set -e
+
+cargo clean
+
+cargo build --manifest-path ../Cargo.toml
+
+mdbook test -L ../target/debug/deps "$@"
+
+mdbook build "$@"
+
+cargo build
+# cargo build --release
--- a/docs/src/SUMMARY.md
+++ b/docs/src/SUMMARY.md
@@ -0,0 +1,7 @@
+# Summary
+
+- [Introduction](./introduction.md)
+- [Data Manipulation](./data-manipulation.md)
+- [Compute Features](./compute.md)
+- [Machine Learning](./machine-learning.md)
+- [Utilities](./utilities.md)
--- a/docs/src/compute.md
+++ b/docs/src/compute.md
@@ -0,0 +1,222 @@
+# Compute Features
+
+The `compute` module hosts numerical routines for exploratory data analysis.
+It covers descriptive statistics, correlations, probability distributions and
+some basic inferential tests.
+
+## Basic Statistics
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::stats::{mean, mean_horizontal, mean_vertical, stddev, median, population_variance, percentile};
+use rustframe::matrix::Matrix;
+
+let m = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+assert_eq!(mean(&m), 2.5);
+assert_eq!(stddev(&m), 1.118033988749895);
+assert_eq!(median(&m), 2.5);
+assert_eq!(population_variance(&m), 1.25);
+assert_eq!(percentile(&m, 50.0), 3.0);
+// column averages returned as 1 x n matrix
+let row_means = mean_horizontal(&m);
+assert_eq!(row_means.data(), &[2.0, 3.0]);
+let col_means = mean_vertical(&m);
+assert_eq!(col_means.data(), & [1.5, 3.5]);
+```
+
+### Axis-specific Operations
+
+Operations can be applied along specific axes (rows or columns):
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::stats::{mean_vertical, mean_horizontal, stddev_vertical, stddev_horizontal};
+use rustframe::matrix::Matrix;
+
+// 3x2 matrix
+let m = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0], 3, 2);
+
+// Mean along columns (vertical) - returns 1 x cols matrix
+let col_means = mean_vertical(&m);
+assert_eq!(col_means.shape(), (1, 2));
+assert_eq!(col_means.data(), &[3.0, 4.0]); // [(1+3+5)/3, (2+4+6)/3]
+
+// Mean along rows (horizontal) - returns rows x 1 matrix
+let row_means = mean_horizontal(&m);
+assert_eq!(row_means.shape(), (3, 1));
+assert_eq!(row_means.data(), &[1.5, 3.5, 5.5]); // [(1+2)/2, (3+4)/2, (5+6)/2]
+
+// Standard deviation along columns
+let col_stddev = stddev_vertical(&m);
+assert_eq!(col_stddev.shape(), (1, 2));
+
+// Standard deviation along rows
+let row_stddev = stddev_horizontal(&m);
+assert_eq!(row_stddev.shape(), (3, 1));
+```
+
+## Correlation
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::stats::{pearson, covariance};
+use rustframe::matrix::Matrix;
+
+let x = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+let y = Matrix::from_vec(vec![2.0, 4.0, 6.0, 8.0], 2, 2);
+let corr = pearson(&x, &y);
+let cov = covariance(&x, &y);
+assert!((corr - 1.0).abs() < 1e-8);
+assert!((cov - 2.5).abs() < 1e-8);
+```
+
+## Covariance
+
+### `covariance`
+
+Computes the population covariance between two equally sized matrices by flattening
+their values.
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::stats::covariance;
+use rustframe::matrix::Matrix;
+
+let x = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+let y = Matrix::from_vec(vec![2.0, 4.0, 6.0, 8.0], 2, 2);
+let cov = covariance(&x, &y);
+assert!((cov - 2.5).abs() < 1e-8);
+```
+
+### `covariance_vertical`
+
+Evaluates covariance between columns (i.e. across rows) and returns a matrix of
+column pair covariances.
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::stats::covariance_vertical;
+use rustframe::matrix::Matrix;
+
+let m = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+let cov = covariance_vertical(&m);
+assert_eq!(cov.shape(), (2, 2));
+assert!(cov.data().iter().all(|&v| (v - 1.0).abs() < 1e-8));
+```
+
+### `covariance_horizontal`
+
+Computes covariance between rows (i.e. across columns) returning a matrix that
+describes how each pair of rows varies together.
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::stats::covariance_horizontal;
+use rustframe::matrix::Matrix;
+
+let m = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+let cov = covariance_horizontal(&m);
+assert_eq!(cov.shape(), (2, 2));
+assert!(cov.data().iter().all(|&v| (v - 0.25).abs() < 1e-8));
+```
+
+### `covariance_matrix`
+
+Builds a covariance matrix either between columns (`Axis::Col`) or rows
+(`Axis::Row`). Each entry represents how two series co-vary.
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::stats::covariance_matrix;
+use rustframe::matrix::{Axis, Matrix};
+
+let data = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+
+// Covariance between columns
+let cov_cols = covariance_matrix(&data, Axis::Col);
+assert!((cov_cols.get(0, 0) - 2.0).abs() < 1e-8);
+
+// Covariance between rows
+let cov_rows = covariance_matrix(&data, Axis::Row);
+assert!((cov_rows.get(0, 1) + 0.5).abs() < 1e-8);
+```
+
+## Distributions
+
+Probability distribution helpers are available for common PDFs and CDFs.
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::stats::distributions::normal_pdf;
+use rustframe::matrix::Matrix;
+
+let x = Matrix::from_vec(vec![0.0, 1.0], 1, 2);
+let pdf = normal_pdf(x, 0.0, 1.0);
+assert_eq!(pdf.data().len(), 2);
+```
+
+### Additional Distributions
+
+Rustframe provides several other probability distributions:
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::stats::distributions::{normal_cdf, binomial_pmf, binomial_cdf, poisson_pmf};
+use rustframe::matrix::Matrix;
+
+// Normal distribution CDF
+let x = Matrix::from_vec(vec![0.0, 1.0], 1, 2);
+let cdf = normal_cdf(x, 0.0, 1.0);
+assert_eq!(cdf.data().len(), 2);
+
+// Binomial distribution PMF
+// Probability of k successes in n trials with probability p
+let k = Matrix::from_vec(vec![0_u64, 1, 2, 3], 1, 4);
+let pmf = binomial_pmf(3, k.clone(), 0.5);
+assert_eq!(pmf.data().len(), 4);
+
+// Binomial distribution CDF
+let cdf = binomial_cdf(3, k, 0.5);
+assert_eq!(cdf.data().len(), 4);
+
+// Poisson distribution PMF
+// Probability of k events with rate parameter lambda
+let k = Matrix::from_vec(vec![0_u64, 1, 2], 1, 3);
+let pmf = poisson_pmf(2.0, k);
+assert_eq!(pmf.data().len(), 3);
+```
+
+### Inferential Statistics
+
+Rustframe provides several inferential statistical tests:
+
+```rust
+# extern crate rustframe;
+use rustframe::matrix::Matrix;
+use rustframe::compute::stats::inferential::{t_test, chi2_test, anova};
+
+// Two-sample t-test
+let sample1 = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0], 1, 5);
+let sample2 = Matrix::from_vec(vec![6.0, 7.0, 8.0, 9.0, 10.0], 1, 5);
+let (t_statistic, p_value) = t_test(&sample1, &sample2);
+assert!((t_statistic + 5.0).abs() < 1e-5);
+assert!(p_value > 0.0 && p_value < 1.0);
+
+// Chi-square test of independence
+let observed = Matrix::from_vec(vec![12.0, 5.0, 8.0, 10.0], 2, 2);
+let (chi2_statistic, p_value) = chi2_test(&observed);
+assert!(chi2_statistic > 0.0);
+assert!(p_value > 0.0 && p_value < 1.0);
+
+// One-way ANOVA
+let group1 = Matrix::from_vec(vec![1.0, 2.0, 3.0], 1, 3);
+let group2 = Matrix::from_vec(vec![2.0, 3.0, 4.0], 1, 3);
+let group3 = Matrix::from_vec(vec![3.0, 4.0, 5.0], 1, 3);
+let groups = vec![&group1, &group2, &group3];
+let (f_statistic, p_value) = anova(groups);
+assert!(f_statistic > 0.0);
+assert!(p_value > 0.0 && p_value < 1.0);
+```
+
+With the basics covered, explore predictive models in the
+[machine learning](./machine-learning.md) chapter.
--- a/docs/src/data-manipulation.md
+++ b/docs/src/data-manipulation.md
@@ -0,0 +1,157 @@
+# Data Manipulation
+
+Rustframe's `Frame` type couples tabular data with
+column labels and a typed row index. Frames expose a familiar API for loading
+data, selecting rows or columns and performing aggregations.
+
+## Creating a Frame
+
+```rust
+# extern crate rustframe;
+use rustframe::frame::{Frame, RowIndex};
+use rustframe::matrix::Matrix;
+
+let data = Matrix::from_cols(vec![vec![1.0, 2.0], vec![3.0, 4.0]]);
+let frame = Frame::new(data, vec!["A", "B"], None);
+assert_eq!(frame["A"], vec![1.0, 2.0]);
+```
+
+## Indexing Rows
+
+Row labels can be integers, dates or a default range. Retrieving a row returns a
+view that lets you inspect values by column name or position.
+
+```rust
+# extern crate rustframe;
+# extern crate chrono;
+use chrono::NaiveDate;
+use rustframe::frame::{Frame, RowIndex};
+use rustframe::matrix::Matrix;
+
+let d = |y, m, d| NaiveDate::from_ymd_opt(y, m, d).unwrap();
+let data = Matrix::from_cols(vec![vec![1.0, 2.0], vec![3.0, 4.0]]);
+let index = RowIndex::Date(vec![d(2024, 1, 1), d(2024, 1, 2)]);
+let mut frame = Frame::new(data, vec!["A", "B"], Some(index));
+assert_eq!(frame.get_row_date(d(2024, 1, 2))["B"], 4.0);
+
+// mutate by row key
+frame.get_row_date_mut(d(2024, 1, 1)).set_by_index(0, 9.0);
+assert_eq!(frame.get_row_date(d(2024, 1, 1))["A"], 9.0);
+```
+
+## Column operations
+
+Columns can be inserted, renamed, removed or reordered in place.
+
+```rust
+# extern crate rustframe;
+use rustframe::frame::{Frame, RowIndex};
+use rustframe::matrix::Matrix;
+
+let data = Matrix::from_cols(vec![vec![1, 2], vec![3, 4]]);
+let mut frame = Frame::new(data, vec!["X", "Y"], Some(RowIndex::Range(0..2)));
+
+frame.add_column("Z", vec![5, 6]);
+frame.rename("Y", "W");
+let removed = frame.delete_column("X");
+assert_eq!(removed, vec![1, 2]);
+frame.sort_columns();
+assert_eq!(frame.columns(), &["W", "Z"]);
+```
+
+## Aggregations
+
+Any numeric aggregation available on `Matrix` is forwarded to `Frame`.
+
+```rust
+# extern crate rustframe;
+use rustframe::frame::Frame;
+use rustframe::matrix::{Matrix, SeriesOps};
+
+let frame = Frame::new(Matrix::from_cols(vec![vec![1.0, 2.0], vec![3.0, 4.0]]), vec!["A", "B"], None);
+assert_eq!(frame.sum_vertical(), vec![3.0, 7.0]);
+assert_eq!(frame.sum_horizontal(), vec![4.0, 6.0]);
+```
+
+## Matrix Operations
+
+```rust
+# extern crate rustframe;
+use rustframe::matrix::Matrix;
+
+let data1 = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+let data2 = Matrix::from_vec(vec![5.0, 6.0, 7.0, 8.0], 2, 2);
+
+let sum = data1.clone() + data2.clone();
+assert_eq!(sum.data(), vec![6.0, 8.0, 10.0, 12.0]);
+
+let product = data1.clone() * data2.clone();
+assert_eq!(product.data(), vec![5.0, 12.0, 21.0, 32.0]);
+
+let scalar_product = data1.clone() * 2.0;
+assert_eq!(scalar_product.data(), vec![2.0, 4.0, 6.0, 8.0]);
+
+let equals = data1 == data1.clone();
+assert_eq!(equals, true);
+```
+
+### Advanced Matrix Operations
+
+Matrices support a variety of advanced operations:
+
+```rust
+# extern crate rustframe;
+use rustframe::matrix::{Matrix, SeriesOps};
+
+// Matrix multiplication (dot product)
+let a = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+let b = Matrix::from_vec(vec![5.0, 6.0, 7.0, 8.0], 2, 2);
+let product = a.matrix_mul(&b);
+assert_eq!(product.data(), vec![23.0, 34.0, 31.0, 46.0]);
+
+// Transpose
+let m = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+let transposed = m.transpose();
+assert_eq!(transposed.data(), vec![1.0, 3.0, 2.0, 4.0]);
+
+// Map function over all elements
+let m = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+let squared = m.map(|x| x * x);
+assert_eq!(squared.data(), vec![1.0, 4.0, 9.0, 16.0]);
+
+// Zip two matrices with a function
+let a = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+let b = Matrix::from_vec(vec![5.0, 6.0, 7.0, 8.0], 2, 2);
+let zipped = a.zip(&b, |x, y| x + y);
+assert_eq!(zipped.data(), vec![6.0, 8.0, 10.0, 12.0]);
+```
+
+### Matrix Reductions
+
+Matrices support various reduction operations:
+
+```rust
+# extern crate rustframe;
+use rustframe::matrix::{Matrix, SeriesOps};
+
+let m = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0], 3, 2);
+
+// Sum along columns (vertical)
+let col_sums = m.sum_vertical();
+assert_eq!(col_sums, vec![9.0, 12.0]); // [1+3+5, 2+4+6]
+
+// Sum along rows (horizontal)
+let row_sums = m.sum_horizontal();
+assert_eq!(row_sums, vec![3.0, 7.0, 11.0]); // [1+2, 3+4, 5+6]
+
+// Cumulative sum along columns
+let col_cumsum = m.cumsum_vertical();
+assert_eq!(col_cumsum.data(), vec![1.0, 4.0, 9.0, 2.0, 6.0, 12.0]);
+
+// Cumulative sum along rows
+let row_cumsum = m.cumsum_horizontal();
+assert_eq!(row_cumsum.data(), vec![1.0, 3.0, 5.0, 3.0, 7.0, 11.0]);
+```
+
+With the basics covered, continue to the [compute features](./compute.md)
+chapter for statistics and analytics.
--- a/docs/src/introduction.md
+++ b/docs/src/introduction.md
@@ -0,0 +1,40 @@
+# Introduction
+
+🐙 [GitHub](https://github.com/Magnus167/rustframe) | 📚 [Docs](https://magnus167.github.io/rustframe/) | 📖 [User Guide](https://magnus167.github.io/rustframe/user-guide/) | 🦀 [Crates.io](https://crates.io/crates/rustframe) | 🔖 [docs.rs](https://docs.rs/rustframe/latest/rustframe/)
+
+Welcome to the **Rustframe User Guide**. Rustframe is a lightweight dataframe
+and math toolkit for Rust written in 100% safe Rust. It focuses on keeping the
+API approachable while offering handy features for small analytical or
+educational projects.
+
+Rustframe bundles:
+
+- column‑labelled frames built on a fast column‑major matrix
+- familiar element‑wise math and aggregation routines
+- a growing `compute` module for statistics and machine learning
+- utilities for dates and random numbers
+
+```rust
+# extern crate rustframe;
+use rustframe::{frame::Frame, matrix::{Matrix, SeriesOps}};
+
+let data = Matrix::from_cols(vec![vec![1.0, 2.0], vec![3.0, 4.0]]);
+let frame = Frame::new(data, vec!["A", "B"], None);
+
+// Perform column wise aggregation
+assert_eq!(frame.sum_vertical(), vec![3.0, 7.0]);
+```
+
+## Resources
+
+- [GitHub repository](https://github.com/Magnus167/rustframe)
+- [Crates.io](https://crates.io/crates/rustframe) & [API docs](https://docs.rs/rustframe)
+- [Code coverage](https://codecov.io/gh/Magnus167/rustframe)
+
+This guide walks through the main building blocks of the library. Each chapter
+contains runnable snippets so you can follow along:
+
+1. [Data manipulation](./data-manipulation.md) for loading and transforming data
+2. [Compute features](./compute.md) for statistics and analytics
+3. [Machine learning](./machine-learning.md) for predictive models
+4. [Utilities](./utilities.md) for supporting helpers and upcoming modules
--- a/docs/src/machine-learning.md
+++ b/docs/src/machine-learning.md
@@ -0,0 +1,282 @@
+# Machine Learning
+
+The `compute::models` module bundles several learning algorithms that operate on
+`Matrix` structures. These examples highlight the basic training and prediction
+APIs. For more end‑to‑end walkthroughs see the examples directory in the
+repository.
+
+Currently implemented models include:
+
+- Linear and logistic regression
+- K‑means clustering
+- Principal component analysis (PCA)
+- Gaussian Naive Bayes
+- Dense neural networks
+
+## Linear Regression
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::models::linreg::LinReg;
+use rustframe::matrix::Matrix;
+
+let x = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 4, 1);
+let y = Matrix::from_vec(vec![2.0, 3.0, 4.0, 5.0], 4, 1);
+let mut model = LinReg::new(1);
+model.fit(&x, &y, 0.01, 100);
+let preds = model.predict(&x);
+assert_eq!(preds.rows(), 4);
+```
+
+## K-means Walkthrough
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::models::k_means::KMeans;
+use rustframe::matrix::Matrix;
+
+let data = Matrix::from_vec(vec![1.0, 1.0, 5.0, 5.0], 2, 2);
+let (model, _labels) = KMeans::fit(&data, 2, 10, 1e-4);
+let new_point = Matrix::from_vec(vec![0.0, 0.0], 1, 2);
+let cluster = model.predict(&new_point)[0];
+```
+
+## Logistic Regression
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::models::logreg::LogReg;
+use rustframe::matrix::Matrix;
+
+let x = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 4, 1);
+let y = Matrix::from_vec(vec![0.0, 0.0, 1.0, 1.0], 4, 1);
+let mut model = LogReg::new(1);
+model.fit(&x, &y, 0.1, 200);
+let preds = model.predict_proba(&x);
+assert_eq!(preds.rows(), 4);
+```
+
+## Principal Component Analysis
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::models::pca::PCA;
+use rustframe::matrix::Matrix;
+
+let data = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+let pca = PCA::fit(&data, 1, 0);
+let transformed = pca.transform(&data);
+assert_eq!(transformed.cols(), 1);
+```
+
+## Gaussian Naive Bayes
+
+Gaussian Naive Bayes classifier for continuous features:
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::models::gaussian_nb::GaussianNB;
+use rustframe::matrix::Matrix;
+
+// Training data with 2 features
+let x = Matrix::from_rows_vec(vec![
+    1.0, 2.0,
+    2.0, 3.0,
+    3.0, 4.0,
+    4.0, 5.0
+], 4, 2);
+
+// Class labels (0 or 1)
+let y = Matrix::from_vec(vec![0.0, 0.0, 1.0, 1.0], 4, 1);
+
+// Train the model
+let mut model = GaussianNB::new(1e-9, true);
+model.fit(&x, &y);
+
+// Make predictions
+let predictions = model.predict(&x);
+assert_eq!(predictions.rows(), 4);
+```
+
+## Dense Neural Networks
+
+Simple fully connected neural network:
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::models::dense_nn::{DenseNN, DenseNNConfig, ActivationKind, InitializerKind, LossKind};
+use rustframe::matrix::Matrix;
+
+// Training data with 2 features
+let x = Matrix::from_rows_vec(vec![
+    0.0, 0.0,
+    0.0, 1.0,
+    1.0, 0.0,
+    1.0, 1.0
+], 4, 2);
+
+// XOR target outputs
+let y = Matrix::from_vec(vec![0.0, 1.0, 1.0, 0.0], 4, 1);
+
+// Create a neural network with 2 hidden layers
+let config = DenseNNConfig {
+    input_size: 2,
+    hidden_layers: vec![4, 4],
+    output_size: 1,
+    activations: vec![ActivationKind::Sigmoid, ActivationKind::Sigmoid, ActivationKind::Sigmoid],
+    initializer: InitializerKind::Uniform(0.5),
+    loss: LossKind::MSE,
+    learning_rate: 0.1,
+    epochs: 1000,
+};
+let mut model = DenseNN::new(config);
+
+// Train the model
+model.train(&x, &y);
+
+// Make predictions
+let predictions = model.predict(&x);
+assert_eq!(predictions.rows(), 4);
+```
+
+## Real-world Examples
+
+### Housing Price Prediction
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::models::linreg::LinReg;
+use rustframe::matrix::Matrix;
+
+// Features: square feet and bedrooms
+let features = Matrix::from_rows_vec(vec![
+    2100.0, 3.0,
+    1600.0, 2.0,
+    2400.0, 4.0,
+    1400.0, 2.0,
+], 4, 2);
+
+// Sale prices
+let target = Matrix::from_vec(vec![400_000.0, 330_000.0, 369_000.0, 232_000.0], 4, 1);
+
+let mut model = LinReg::new(2);
+model.fit(&features, &target, 1e-8, 10_000);
+
+// Predict price of a new home
+let new_home = Matrix::from_vec(vec![2000.0, 3.0], 1, 2);
+let predicted_price = model.predict(&new_home);
+println!("Predicted price: ${}", predicted_price.data()[0]);
+```
+
+### Spam Detection
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::models::logreg::LogReg;
+use rustframe::matrix::Matrix;
+
+// 20 e-mails × 5 features = 100 numbers (row-major, spam first)
+let x = Matrix::from_rows_vec(
+    vec![
+        // ─────────── spam examples ───────────
+        2.0, 1.0, 1.0, 1.0, 1.0, // "You win a FREE offer - click for money-back bonus!"
+        1.0, 0.0, 1.0, 1.0, 0.0, // "FREE offer! Click now!"
+        0.0, 2.0, 0.0, 1.0, 1.0, // "Win win win - money inside, click…"
+        1.0, 1.0, 0.0, 0.0, 1.0, // "Limited offer to win easy money…"
+        1.0, 0.0, 1.0, 0.0, 1.0, // ...
+        0.0, 1.0, 1.0, 1.0, 0.0, // ...
+        2.0, 0.0, 0.0, 1.0, 1.0, // ...
+        0.0, 1.0, 1.0, 0.0, 1.0, // ...
+        1.0, 1.0, 1.0, 1.0, 0.0, // ...
+        1.0, 0.0, 0.0, 1.0, 1.0, // ...
+        // ─────────── ham examples ───────────
+        0.0, 0.0, 0.0, 0.0, 0.0, // "See you at the meeting tomorrow."
+        0.0, 0.0, 0.0, 1.0, 0.0, // "Here's the Zoom click-link."
+        0.0, 0.0, 0.0, 0.0, 1.0, // "Expense report: money attached."
+        0.0, 0.0, 0.0, 1.0, 1.0, // ...
+        0.0, 1.0, 0.0, 0.0, 0.0, // "Did we win the bid?"
+        0.0, 0.0, 0.0, 0.0, 0.0, // ...
+        0.0, 0.0, 0.0, 1.0, 0.0, // ...
+        1.0, 0.0, 0.0, 0.0, 0.0, // "Special offer for staff lunch."
+        0.0, 0.0, 0.0, 0.0, 0.0, // ...
+        0.0, 0.0, 0.0, 1.0, 0.0,
+    ],
+    20,
+    5,
+);
+
+// Labels: 1 = spam, 0 = ham
+let y = Matrix::from_vec(
+    vec![
+        1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, // 10 spam
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // 10 ham
+    ],
+    20,
+    1,
+);
+
+// Train
+let mut model = LogReg::new(5);
+model.fit(&x, &y, 0.01, 5000);
+
+// Predict
+// e.g. "free money offer"
+let email_data = vec![1.0, 0.0, 1.0, 0.0, 1.0];
+let email = Matrix::from_vec(email_data, 1, 5);
+let prob_spam = model.predict_proba(&email);
+println!("Probability of spam: {:.4}", prob_spam.data()[0]);
+```
+
+### Iris Flower Classification
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::models::gaussian_nb::GaussianNB;
+use rustframe::matrix::Matrix;
+
+// Features: sepal length and petal length
+let x = Matrix::from_rows_vec(vec![
+    5.1, 1.4, // setosa
+    4.9, 1.4, // setosa
+    6.2, 4.5, // versicolor
+    5.9, 5.1, // virginica
+], 4, 2);
+
+let y = Matrix::from_vec(vec![0.0, 0.0, 1.0, 2.0], 4, 1);
+let names = vec!["setosa", "versicolor", "virginica"];
+
+let mut model = GaussianNB::new(1e-9, true);
+model.fit(&x, &y);
+
+let sample = Matrix::from_vec(vec![5.0, 1.5], 1, 2);
+let predicted_class = model.predict(&sample);
+let class_name = names[predicted_class.data()[0] as usize];
+println!("Predicted class: {} ({:?})", class_name, predicted_class.data()[0]);
+```
+
+### Customer Segmentation
+
+```rust
+# extern crate rustframe;
+use rustframe::compute::models::k_means::KMeans;
+use rustframe::matrix::Matrix;
+
+// Each row: [age, annual_income]
+let customers = Matrix::from_rows_vec(
+    vec![
+        25.0, 40_000.0, 34.0, 52_000.0, 58.0, 95_000.0, 45.0, 70_000.0,
+    ],
+    4,
+    2,
+);
+
+let (model, labels) = KMeans::fit(&customers, 2, 20, 1e-4);
+
+let new_customer = Matrix::from_vec(vec![30.0, 50_000.0], 1, 2);
+let cluster = model.predict(&new_customer)[0];
+println!("New customer belongs to cluster: {}", cluster);
+println!("Cluster labels: {:?}", labels);
+```
+
+For helper functions and upcoming modules, visit the
+[utilities](./utilities.md) section.
--- a/docs/src/utilities.md
+++ b/docs/src/utilities.md
@@ -0,0 +1,63 @@
+# Utilities
+
+Utilities provide handy helpers around the core library. Existing tools
+include:
+
+- Date utilities for generating calendar sequences and business‑day sets
+- Random number generators for simulations and testing
+
+## Date Helpers
+
+```rust
+# extern crate rustframe;
+use rustframe::utils::dateutils::{BDatesList, BDateFreq, DatesList, DateFreq};
+
+// Calendar sequence
+let list = DatesList::new("2024-01-01".into(), "2024-01-03".into(), DateFreq::Daily);
+assert_eq!(list.count().unwrap(), 3);
+
+// Business days starting from 2024‑01‑02
+let bdates = BDatesList::from_n_periods("2024-01-02".into(), BDateFreq::Daily, 3).unwrap();
+assert_eq!(bdates.list().unwrap().len(), 3);
+```
+
+## Random Numbers
+
+The `random` module offers deterministic and cryptographically secure RNGs.
+
+```rust
+# extern crate rustframe;
+use rustframe::random::{Prng, Rng};
+
+let mut rng = Prng::new(42);
+let v1 = rng.next_u64();
+let v2 = rng.next_u64();
+assert_ne!(v1, v2);
+```
+
+## Stats Functions
+
+```rust
+# extern crate rustframe;
+use rustframe::matrix::Matrix;
+use rustframe::compute::stats::descriptive::{mean, median, stddev};
+
+let data = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0], 1, 5);
+
+let mean_value = mean(&data);
+assert_eq!(mean_value, 3.0);
+
+let median_value = median(&data);
+assert_eq!(median_value, 3.0);
+
+let std_value = stddev(&data);
+assert_eq!(std_value, 2.0_f64.sqrt());
+```
+
+Upcoming utilities will cover:
+
+- Data import/export helpers
+- Visualization adapters
+- Streaming data interfaces
+
+Contributions to these sections are welcome!
--- a/examples/correlation.rs
+++ b/examples/correlation.rs
@@ -0,0 +1,45 @@
+use rustframe::compute::stats::{covariance, covariance_matrix, pearson};
+use rustframe::matrix::{Axis, Matrix};
+
+/// Demonstrates covariance and correlation utilities.
+fn main() {
+    pairwise_cov();
+    println!("\n-----\n");
+    matrix_cov();
+}
+
+fn pairwise_cov() {
+    println!("Covariance & Pearson r\n----------------------");
+    let x = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+    let y = Matrix::from_vec(vec![1.0, 2.0, 3.0, 5.0], 2, 2);
+    println!("covariance : {:.2}", covariance(&x, &y));
+    println!("pearson r  : {:.3}", pearson(&x, &y));
+}
+
+fn matrix_cov() {
+    println!("Covariance matrix\n-----------------");
+    let data = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+    let cov = covariance_matrix(&data, Axis::Col);
+    println!("cov matrix : {:?}", cov.data());
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    const EPS: f64 = 1e-8;
+
+    #[test]
+    fn test_pairwise_cov() {
+        let x = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+        let y = Matrix::from_vec(vec![1.0, 2.0, 3.0, 5.0], 2, 2);
+        assert!((covariance(&x, &y) - 1.625).abs() < EPS);
+        assert!((pearson(&x, &y) - 0.9827076298239908).abs() < 1e-5,);
+    }
+
+    #[test]
+    fn test_matrix_cov() {
+        let data = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+        let cov = covariance_matrix(&data, Axis::Col);
+        assert_eq!(cov.data(), &[2.0, 2.0, 2.0, 2.0]);
+    }
+}
--- a/examples/descriptive_stats.rs
+++ b/examples/descriptive_stats.rs
@@ -0,0 +1,56 @@
+use rustframe::compute::stats::{mean, mean_horizontal, mean_vertical, median, percentile, stddev};
+use rustframe::matrix::Matrix;
+
+/// Demonstrates descriptive statistics utilities.
+///
+/// Part 1: simple mean/stddev/median/percentile on a vector.
+/// Part 2: mean across rows and columns.
+fn main() {
+    simple_stats();
+    println!("\n-----\n");
+    axis_stats();
+}
+
+fn simple_stats() {
+    println!("Basic stats\n-----------");
+    let data = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0], 1, 5);
+    println!("mean      : {:.2}", mean(&data));
+    println!("stddev    : {:.2}", stddev(&data));
+    println!("median    : {:.2}", median(&data));
+    println!("90th pct. : {:.2}", percentile(&data, 90.0));
+}
+
+fn axis_stats() {
+    println!("Row/column means\n----------------");
+    // 2x3 matrix
+    let data = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0], 2, 3);
+    let v = mean_vertical(&data); // 1x3
+    let h = mean_horizontal(&data); // 2x1
+    println!("vertical means  : {:?}", v.data());
+    println!("horizontal means: {:?}", h.data());
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    const EPS: f64 = 1e-8;
+
+    #[test]
+    fn test_simple_stats() {
+        let data = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0], 1, 5);
+        assert!((mean(&data) - 3.0).abs() < EPS);
+        assert!((stddev(&data) - 1.4142135623730951).abs() < EPS);
+        assert!((median(&data) - 3.0).abs() < EPS);
+        assert!((percentile(&data, 90.0) - 5.0).abs() < EPS);
+    }
+
+    #[test]
+    fn test_axis_stats() {
+        let data = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0], 2, 3);
+        let v = mean_vertical(&data);
+        assert_eq!(v.data(), &[2.5, 3.5, 4.5]);
+        let h = mean_horizontal(&data);
+        assert_eq!(h.data(), &[2.0, 5.0]);
+    }
+}
+
--- a/examples/distributions.rs
+++ b/examples/distributions.rs
@@ -0,0 +1,66 @@
+use rustframe::compute::stats::{binomial_cdf, binomial_pmf, normal_cdf, normal_pdf, poisson_pmf};
+use rustframe::matrix::Matrix;
+
+/// Demonstrates some probability distribution helpers.
+fn main() {
+    normal_example();
+    println!("\n-----\n");
+    binomial_example();
+    println!("\n-----\n");
+    poisson_example();
+}
+
+fn normal_example() {
+    println!("Normal distribution\n-------------------");
+    let x = Matrix::from_vec(vec![0.0, 1.0], 1, 2);
+    let pdf = normal_pdf(x.clone(), 0.0, 1.0);
+    let cdf = normal_cdf(x, 0.0, 1.0);
+    println!("pdf : {:?}", pdf.data());
+    println!("cdf : {:?}", cdf.data());
+}
+
+fn binomial_example() {
+    println!("Binomial distribution\n---------------------");
+    let k = Matrix::from_vec(vec![0_u64, 1, 2], 1, 3);
+    let pmf = binomial_pmf(4, k.clone(), 0.5);
+    let cdf = binomial_cdf(4, k, 0.5);
+    println!("pmf : {:?}", pmf.data());
+    println!("cdf : {:?}", cdf.data());
+}
+
+fn poisson_example() {
+    println!("Poisson distribution\n--------------------");
+    let k = Matrix::from_vec(vec![0_u64, 1, 2], 1, 3);
+    let pmf = poisson_pmf(3.0, k);
+    println!("pmf : {:?}", pmf.data());
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_normal_example() {
+        let x = Matrix::from_vec(vec![0.0, 1.0], 1, 2);
+        let pdf = normal_pdf(x.clone(), 0.0, 1.0);
+        let cdf = normal_cdf(x, 0.0, 1.0);
+        assert!((pdf.get(0, 0) - 0.39894228).abs() < 1e-6);
+        assert!((cdf.get(0, 1) - 0.8413447).abs() < 1e-6);
+    }
+
+    #[test]
+    fn test_binomial_example() {
+        let k = Matrix::from_vec(vec![0_u64, 1, 2], 1, 3);
+        let pmf = binomial_pmf(4, k.clone(), 0.5);
+        let cdf = binomial_cdf(4, k, 0.5);
+        assert!((pmf.get(0, 2) - 0.375).abs() < 1e-6);
+        assert!((cdf.get(0, 2) - 0.6875).abs() < 1e-6);
+    }
+
+    #[test]
+    fn test_poisson_example() {
+        let k = Matrix::from_vec(vec![0_u64, 1, 2], 1, 3);
+        let pmf = poisson_pmf(3.0, k);
+        assert!((pmf.get(0, 1) - 3.0_f64 * (-3.0_f64).exp()).abs() < 1e-6);
+    }
+}
--- a/examples/game_of_life.rs
+++ b/examples/game_of_life.rs
@@ -1,13 +1,26 @@
-use rand::{self, Rng};
+//! Conway's Game of Life Example
+//! This example implements Conway's Game of Life using a `BoolMatrix` to represent the game board.
+//! It demonstrates matrix operations like shifting, counting neighbors, and applying game rules.
+//! The game runs in a loop, updating the board state and printing it to the console.
+//! To modify the behaviour of the example, please change the constants at the top of this file.
+
+
 use rustframe::matrix::{BoolMatrix, BoolOps, IntMatrix, Matrix};
+use rustframe::random::{rng, Rng};
 use std::{thread, time};

-const BOARD_SIZE: usize = 50; // Size of the board (50x50)
-const TICK_DURATION_MS: u64 = 10; // Milliseconds per frame
+const BOARD_SIZE: usize = 20; // Size of the board (50x50)
+const MAX_FRAMES: u32 = 1000;
+
+const TICK_DURATION_MS: u64 = 0; // Milliseconds per frame
+const SKIP_FRAMES: u32 = 1;
+const PRINT_BOARD: bool = true; // Set to false to disable printing the board

 fn main() {
-    // Initialize the game board.
-    // This demonstrates `BoolMatrix::from_vec`.
+    let args = std::env::args().collect::<Vec<String>>();
+    let debug_mode = args.contains(&"--debug".to_string());
+    let print_mode = if debug_mode { false } else { PRINT_BOARD };
+
    let mut current_board =
        BoolMatrix::from_vec(vec![false; BOARD_SIZE * BOARD_SIZE], BOARD_SIZE, BOARD_SIZE);

@@ -16,31 +29,18 @@ fn main() {
    add_simulated_activity(&mut current_board, BOARD_SIZE);

    let mut generation_count: u32 = 0;
-    // `previous_board_state` will store a clone of the board.
-    // This demonstrates `Matrix::clone()` and later `PartialEq` for `Matrix`.
    let mut previous_board_state: Option<BoolMatrix> = None;
    let mut board_hashes = Vec::new();
-    // let mut print_board_bool = true;
    let mut print_bool_int = 0;

    loop {
-        // print!("{}[2J", 27 as char); // Clear screen and move cursor to top-left
+        if print_bool_int % SKIP_FRAMES == 0 {
+            print_board(&current_board, generation_count, print_mode);

-        // if print_board_bool {
-        if print_bool_int % 10 == 0 {
-            print!("{}[2J", 27 as char);
-            println!("Conway's Game of Life - Generation: {}", generation_count);
-
-            print_board(&current_board);
-            println!("Alive cells: {}", &current_board.count());
-
-            // print_board_bool = false;
            print_bool_int = 0;
        } else {
-            // print_board_bool = true;
            print_bool_int += 1;
        }
-        // `current_board.count()` demonstrates a method from `BoolOps`.
        board_hashes.push(hash_board(&current_board, primes.clone()));
        if detect_stable_state(&current_board, &previous_board_state) {
            println!(
@@ -61,20 +61,18 @@ fn main() {
            add_simulated_activity(&mut current_board, BOARD_SIZE);
        }

-        // `current_board.clone()` demonstrates `Clone` for `Matrix`.
        previous_board_state = Some(current_board.clone());

-        // This is the core call to your game logic.
        let next_board = game_of_life_next_frame(&current_board);
        current_board = next_board;

        generation_count += 1;
        thread::sleep(time::Duration::from_millis(TICK_DURATION_MS));

-        // if generation_count > 500 { // Optional limit
-        //     println!("\nReached generation limit.");
-        //     break;
-        // }
+        if (MAX_FRAMES > 0) && (generation_count > MAX_FRAMES) {
+            println!("\nReached generation limit.");
+            break;
+        }
    }
 }

@@ -82,7 +80,13 @@ fn main() {
 ///
 /// - `board`: A reference to the `BoolMatrix` representing the current game state.
 /// This function demonstrates `board.rows()`, `board.cols()`, and `board[(r, c)]` (Index trait).
-fn print_board(board: &BoolMatrix) {
+fn print_board(board: &BoolMatrix, generation_count: u32, print_mode: bool) {
+    if !print_mode {
+        return;
+    }
+
+    print!("{}[2J", 27 as char);
+    println!("Conway's Game of Life - Generation: {}", generation_count);
    let mut print_str = String::new();
    print_str.push_str("+");
    for _ in 0..board.cols() {
@@ -93,7 +97,6 @@ fn print_board(board: &BoolMatrix) {
        print_str.push_str("| ");
        for c in 0..board.cols() {
            if board[(r, c)] {
-                // Using Index trait for Matrix<bool>
                print_str.push_str("██");
            } else {
                print_str.push_str("  ");
@@ -107,6 +110,8 @@ fn print_board(board: &BoolMatrix) {
    }
    print_str.push_str("+\n\n");
    print!("{}", print_str);
+
+    println!("Alive cells: {}", board.count());
 }

 /// Helper function to create a shifted version of the game board.
@@ -173,74 +178,38 @@ pub fn game_of_life_next_frame(current_game: &BoolMatrix) -> BoolMatrix {
    if rows == 0 && cols == 0 {
        return BoolMatrix::from_vec(vec![], 0, 0); // Return an empty BoolMatrix
    }
-    // Assuming valid non-empty dimensions (e.g., 25x25) as per typical GOL.
-    // Your Matrix::from_vec would panic for other invalid 0-dim cases.
-
    // Define the 8 neighbor offsets (row_delta, col_delta)
    let neighbor_offsets: [(isize, isize); 8] = [
        (-1, -1),
        (-1, 0),
-        (-1, 1), // Top row (NW, N, NE)
+        (-1, 1),
        (0, -1),
-        (0, 1), // Middle row (W, E)
+        (0, 1),
        (1, -1),
        (1, 0),
-        (1, 1), // Bottom row (SW, S, SE)
+        (1, 1),
    ];

-    // 1. Initialize `neighbor_counts` with the first shifted layer.
-    //    This demonstrates creating an IntMatrix from a function and using it as a base.
    let (first_dr, first_dc) = neighbor_offsets[0];
    let mut neighbor_counts = get_shifted_neighbor_layer(current_game, first_dr, first_dc);

-    // 2. Add the remaining 7 neighbor layers.
-    //    This demonstrates element-wise addition of matrices (`Matrix + Matrix`).
    for i in 1..neighbor_offsets.len() {
        let (dr, dc) = neighbor_offsets[i];
        let next_neighbor_layer = get_shifted_neighbor_layer(current_game, dr, dc);
-        // `neighbor_counts` (owned IntMatrix) + `next_neighbor_layer` (owned IntMatrix)
-        // uses `impl Add for Matrix`, consumes both, returns new owned `IntMatrix`.
        neighbor_counts = neighbor_counts + next_neighbor_layer;
    }

-    // 3. Apply Game of Life rules using element-wise operations.
-
-    // Rule: Survival or Birth based on neighbor counts.
-    // A cell is alive in the next generation if:
-    //   (it's currently alive AND has 2 or 3 neighbors) OR
-    //   (it's currently dead AND has exactly 3 neighbors)
-
-    // `neighbor_counts.eq_elem(scalar)`:
-    //   Demonstrates element-wise comparison of a Matrix with a scalar (broadcast).
-    //   Returns an owned `BoolMatrix`.
    let has_2_neighbors = neighbor_counts.eq_elem(2);
-    let has_3_neighbors = neighbor_counts.eq_elem(3); // This will be reused
+    let has_3_neighbors = neighbor_counts.eq_elem(3);

-    // `has_2_neighbors | has_3_neighbors`:
-    //   Demonstrates element-wise OR (`Matrix<bool> | Matrix<bool>`).
-    //   Consumes both operands, returns an owned `BoolMatrix`.
-    let has_2_or_3_neighbors = has_2_neighbors | has_3_neighbors.clone(); // Clone has_3_neighbors as it's used again
+    let has_2_or_3_neighbors = has_2_neighbors | has_3_neighbors.clone();

-    // `current_game & &has_2_or_3_neighbors`:
-    //   `current_game` is `&BoolMatrix`. `has_2_or_3_neighbors` is owned.
-    //   Demonstrates element-wise AND (`&Matrix<bool> & &Matrix<bool>`).
-    //   Borrows both operands, returns an owned `BoolMatrix`.
    let survives = current_game & &has_2_or_3_neighbors;

-    // `!current_game`:
-    //   Demonstrates element-wise NOT (`!&Matrix<bool>`).
-    //   Borrows operand, returns an owned `BoolMatrix`.
    let is_dead = !current_game;

-    // `is_dead & &has_3_neighbors`:
-    //   `is_dead` is owned. `has_3_neighbors` is owned.
-    //   Demonstrates element-wise AND (`Matrix<bool> & &Matrix<bool>`).
-    //   Consumes `is_dead`, borrows `has_3_neighbors`, returns an owned `BoolMatrix`.
    let births = is_dead & &has_3_neighbors;

-    // `survives | births`:
-    //   Demonstrates element-wise OR (`Matrix<bool> | Matrix<bool>`).
-    //   Consumes both operands, returns an owned `BoolMatrix`.
    let next_frame_game = survives | births;

    next_frame_game
@@ -250,7 +219,7 @@ pub fn generate_glider(board: &mut BoolMatrix, board_size: usize) {
    // Initialize with a Glider pattern.
    // It demonstrates how to set specific cells in the matrix.
    // This demonstrates `IndexMut` for `current_board[(r, c)] = true;`.
-    let mut rng = rand::rng();
+    let mut rng = rng();
    let r_offset = rng.random_range(0..(board_size - 3));
    let c_offset = rng.random_range(0..(board_size - 3));
    if board.rows() >= r_offset + 3 && board.cols() >= c_offset + 3 {
@@ -266,7 +235,7 @@ pub fn generate_pulsar(board: &mut BoolMatrix, board_size: usize) {
    // Initialize with a Pulsar pattern.
    // This demonstrates how to set specific cells in the matrix.
    // This demonstrates `IndexMut` for `current_board[(r, c)] = true;`.
-    let mut rng = rand::rng();
+    let mut rng = rng();
    let r_offset = rng.random_range(0..(board_size - 17));
    let c_offset = rng.random_range(0..(board_size - 17));
    if board.rows() >= r_offset + 17 && board.cols() >= c_offset + 17 {
--- a/examples/inferential_stats.rs
+++ b/examples/inferential_stats.rs
@@ -0,0 +1,66 @@
+use rustframe::compute::stats::{anova, chi2_test, t_test};
+use rustframe::matrix::Matrix;
+
+/// Demonstrates simple inferential statistics tests.
+fn main() {
+    t_test_demo();
+    println!("\n-----\n");
+    chi2_demo();
+    println!("\n-----\n");
+    anova_demo();
+}
+
+fn t_test_demo() {
+    println!("Two-sample t-test\n-----------------");
+    let a = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0], 1, 5);
+    let b = Matrix::from_vec(vec![6.0, 7.0, 8.0, 9.0, 10.0], 1, 5);
+    let (t, p) = t_test(&a, &b);
+    println!("t statistic: {:.2}, p-value: {:.4}", t, p);
+}
+
+fn chi2_demo() {
+    println!("Chi-square test\n---------------");
+    let observed = Matrix::from_vec(vec![12.0, 5.0, 8.0, 10.0], 2, 2);
+    let (chi2, p) = chi2_test(&observed);
+    println!("chi^2: {:.2}, p-value: {:.4}", chi2, p);
+}
+
+fn anova_demo() {
+    println!("One-way ANOVA\n-------------");
+    let g1 = Matrix::from_vec(vec![1.0, 2.0, 3.0], 1, 3);
+    let g2 = Matrix::from_vec(vec![2.0, 3.0, 4.0], 1, 3);
+    let g3 = Matrix::from_vec(vec![3.0, 4.0, 5.0], 1, 3);
+    let (f, p) = anova(vec![&g1, &g2, &g3]);
+    println!("F statistic: {:.2}, p-value: {:.4}", f, p);
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_t_test_demo() {
+        let a = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0], 1, 5);
+        let b = Matrix::from_vec(vec![6.0, 7.0, 8.0, 9.0, 10.0], 1, 5);
+        let (t, _p) = t_test(&a, &b);
+        assert!((t + 5.0).abs() < 1e-5);
+    }
+
+    #[test]
+    fn test_chi2_demo() {
+        let observed = Matrix::from_vec(vec![12.0, 5.0, 8.0, 10.0], 2, 2);
+        let (chi2, p) = chi2_test(&observed);
+        assert!(chi2 > 0.0);
+        assert!(p > 0.0 && p < 1.0);
+    }
+
+    #[test]
+    fn test_anova_demo() {
+        let g1 = Matrix::from_vec(vec![1.0, 2.0, 3.0], 1, 3);
+        let g2 = Matrix::from_vec(vec![2.0, 3.0, 4.0], 1, 3);
+        let g3 = Matrix::from_vec(vec![3.0, 4.0, 5.0], 1, 3);
+        let (f, p) = anova(vec![&g1, &g2, &g3]);
+        assert!(f > 0.0);
+        assert!(p > 0.0 && p < 1.0);
+    }
+}
--- a/examples/k_means.rs
+++ b/examples/k_means.rs
@@ -0,0 +1,65 @@
+use rustframe::compute::models::k_means::KMeans;
+use rustframe::matrix::Matrix;
+
+/// Two quick K-Means clustering demos.
+///
+/// Example 1 groups store locations on a city map.
+/// Example 2 segments customers by annual spending habits.
+fn main() {
+    city_store_example();
+    println!("\n-----\n");
+    customer_spend_example();
+}
+
+fn city_store_example() {
+    println!("Example 1: store locations");
+
+    // (x, y) coordinates of stores around a city
+    let raw = vec![
+        1.0, 2.0, 1.5, 1.8, 5.0, 8.0, 8.0, 8.0, 1.0, 0.6, 9.0, 11.0, 8.0, 2.0, 10.0, 2.0, 9.0, 3.0,
+    ];
+    let x = Matrix::from_rows_vec(raw, 9, 2);
+
+    // Group stores into two areas
+    let (model, labels) = KMeans::fit(&x, 2, 100, 1e-4);
+
+    println!("Centres: {:?}", model.centroids.data());
+    println!("Labels: {:?}", labels);
+
+    let new_points = Matrix::from_rows_vec(vec![0.0, 0.0, 8.0, 3.0], 2, 2);
+    let pred = model.predict(&new_points);
+    println!("New store assignments: {:?}", pred);
+}
+
+fn customer_spend_example() {
+    println!("Example 2: customer spending");
+
+    // (grocery spend, electronics spend) in dollars
+    let raw = vec![
+        200.0, 150.0, 220.0, 170.0, 250.0, 160.0, 800.0, 750.0, 820.0, 760.0, 790.0, 770.0,
+    ];
+    let x = Matrix::from_rows_vec(raw, 6, 2);
+
+    let (model, labels) = KMeans::fit(&x, 2, 100, 1e-4);
+
+    println!("Centres: {:?}", model.centroids.data());
+    println!("Labels: {:?}", labels);
+
+    let new_customers = Matrix::from_rows_vec(vec![230.0, 155.0, 810.0, 760.0], 2, 2);
+    let pred = model.predict(&new_customers);
+    println!("Cluster of new customers: {:?}", pred);
+}
+
+#[test]
+fn k_means_store_locations() {
+    let raw = vec![
+        1.0, 2.0, 1.5, 1.8, 5.0, 8.0, 8.0, 8.0, 1.0, 0.6, 9.0, 11.0, 8.0, 2.0, 10.0, 2.0, 9.0, 3.0,
+    ];
+    let x = Matrix::from_rows_vec(raw, 9, 2);
+    let (model, labels) = KMeans::fit(&x, 2, 100, 1e-4);
+    assert_eq!(labels.len(), 9);
+    assert_eq!(model.centroids.rows(), 2);
+    let new_points = Matrix::from_rows_vec(vec![0.0, 0.0, 8.0, 3.0], 2, 2);
+    let pred = model.predict(&new_points);
+    assert_eq!(pred.len(), 2);
+}
--- a/examples/linear_regression.rs
+++ b/examples/linear_regression.rs
@@ -0,0 +1,118 @@
+use rustframe::compute::models::linreg::LinReg;
+use rustframe::matrix::Matrix;
+
+/// Two quick linear regression demonstrations.
+///
+/// Example 1 fits a model to predict house price from floor area.
+/// Example 2 adds number of bedrooms as a second feature.
+fn main() {
+    example_one_feature();
+    println!("\n-----\n");
+    example_two_features();
+}
+
+/// Price ~ floor area
+fn example_one_feature() {
+    println!("Example 1: predict price from floor area only");
+
+    // Square meters of floor area for a few houses
+    let sizes = vec![50.0, 60.0, 70.0, 80.0, 90.0, 100.0];
+    // Thousands of dollars in sale price
+    let prices = vec![150.0, 180.0, 210.0, 240.0, 270.0, 300.0];
+
+    // Each row is a sample with one feature
+    let x = Matrix::from_vec(sizes.clone(), sizes.len(), 1);
+    let y = Matrix::from_vec(prices.clone(), prices.len(), 1);
+
+    // Train with a small learning rate
+    let mut model = LinReg::new(1);
+    model.fit(&x, &y, 0.0005, 20000);
+
+    let preds = model.predict(&x);
+    println!("Size (m^2) -> predicted price (k) vs actual");
+    for i in 0..x.rows() {
+        println!(
+            "{:>3} -> {:>6.1} | {:>6.1}",
+            sizes[i],
+            preds[(i, 0)],
+            prices[i]
+        );
+    }
+
+    let new_house = Matrix::from_vec(vec![120.0], 1, 1);
+    let pred = model.predict(&new_house);
+    println!("Predicted price for 120 m^2: {:.1}k", pred[(0, 0)]);
+}
+
+/// Price ~ floor area + bedrooms
+fn example_two_features() {
+    println!("Example 2: price from area and bedrooms");
+
+    // (size m^2, bedrooms) for each house
+    let raw_x = vec![
+        50.0, 2.0, 70.0, 2.0, 90.0, 3.0, 110.0, 3.0, 130.0, 4.0, 150.0, 4.0,
+    ];
+    let prices = vec![160.0, 195.0, 250.0, 285.0, 320.0, 350.0];
+
+    let x = Matrix::from_rows_vec(raw_x, 6, 2);
+    let y = Matrix::from_vec(prices.clone(), prices.len(), 1);
+
+    let mut model = LinReg::new(2);
+    model.fit(&x, &y, 0.0001, 50000);
+
+    let preds = model.predict(&x);
+    println!("size, beds -> predicted | actual (k)");
+    for i in 0..x.rows() {
+        let size = x[(i, 0)];
+        let beds = x[(i, 1)];
+        println!(
+            "{:>3} m^2, {:>1} -> {:>6.1} | {:>6.1}",
+            size,
+            beds,
+            preds[(i, 0)],
+            prices[i]
+        );
+    }
+
+    let new_home = Matrix::from_rows_vec(vec![120.0, 3.0], 1, 2);
+    let pred = model.predict(&new_home);
+    println!(
+        "Predicted price for 120 m^2 with 3 bedrooms: {:.1}k",
+        pred[(0, 0)]
+    );
+}
+
+#[test]
+fn test_linear_regression_one_feature() {
+    let sizes = vec![50.0, 60.0, 70.0, 80.0, 90.0, 100.0];
+    let prices = vec![150.0, 180.0, 210.0, 240.0, 270.0, 300.0];
+    let scaled: Vec<f64> = sizes.iter().map(|s| s / 100.0).collect();
+    let x = Matrix::from_vec(scaled, sizes.len(), 1);
+    let y = Matrix::from_vec(prices.clone(), prices.len(), 1);
+    let mut model = LinReg::new(1);
+    model.fit(&x, &y, 0.1, 2000);
+    let preds = model.predict(&x);
+    for i in 0..y.rows() {
+        assert!((preds[(i, 0)] - prices[i]).abs() < 1.0);
+    }
+}
+
+#[test]
+fn test_linear_regression_two_features() {
+    let raw_x = vec![
+        50.0, 2.0, 70.0, 2.0, 90.0, 3.0, 110.0, 3.0, 130.0, 4.0, 150.0, 4.0,
+    ];
+    let prices = vec![170.0, 210.0, 270.0, 310.0, 370.0, 410.0];
+    let scaled_x: Vec<f64> = raw_x
+        .chunks(2)
+        .flat_map(|pair| vec![pair[0] / 100.0, pair[1]])
+        .collect();
+    let x = Matrix::from_rows_vec(scaled_x, 6, 2);
+    let y = Matrix::from_vec(prices.clone(), prices.len(), 1);
+    let mut model = LinReg::new(2);
+    model.fit(&x, &y, 0.01, 50000);
+    let preds = model.predict(&x);
+    for i in 0..y.rows() {
+        assert!((preds[(i, 0)] - prices[i]).abs() < 1.0);
+    }
+}
--- a/examples/logistic_regression.rs
+++ b/examples/logistic_regression.rs
@@ -0,0 +1,101 @@
+use rustframe::compute::models::logreg::LogReg;
+use rustframe::matrix::Matrix;
+
+/// Two binary classification demos using logistic regression.
+///
+/// Example 1 predicts exam success from hours studied.
+/// Example 2 predicts whether an online shopper will make a purchase.
+fn main() {
+    student_passing_example();
+    println!("\n-----\n");
+    purchase_prediction_example();
+}
+
+fn student_passing_example() {
+    println!("Example 1: exam pass prediction");
+
+    // Hours studied for each student
+    let hours = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0];
+    // Label: 0 denotes failure and 1 denotes success
+    let passed = vec![0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0];
+
+    let x = Matrix::from_vec(hours.clone(), hours.len(), 1);
+    let y = Matrix::from_vec(passed.clone(), passed.len(), 1);
+
+    let mut model = LogReg::new(1);
+    model.fit(&x, &y, 0.1, 10000);
+
+    let preds = model.predict(&x);
+    println!("Hours -> pred | actual");
+    for i in 0..x.rows() {
+        println!(
+            "{:>2} -> {} | {}",
+            hours[i] as i32,
+            preds[(i, 0)] as i32,
+            passed[i] as i32
+        );
+    }
+
+    // Probability estimate for a new student
+    let new_student = Matrix::from_vec(vec![5.5], 1, 1);
+    let p = model.predict_proba(&new_student);
+    println!("Probability of passing with 5.5h study: {:.2}", p[(0, 0)]);
+}
+
+fn purchase_prediction_example() {
+    println!("Example 2: purchase likelihood");
+
+    // minutes on site, pages viewed -> made a purchase?
+    let raw_x = vec![1.0, 2.0, 3.0, 1.0, 2.0, 4.0, 5.0, 5.0, 3.5, 2.0, 6.0, 6.0];
+    let bought = vec![0.0, 0.0, 0.0, 1.0, 0.0, 1.0];
+
+    let x = Matrix::from_rows_vec(raw_x, 6, 2);
+    let y = Matrix::from_vec(bought.clone(), bought.len(), 1);
+
+    let mut model = LogReg::new(2);
+    model.fit(&x, &y, 0.05, 20000);
+
+    let preds = model.predict(&x);
+    println!("time, pages -> pred | actual");
+    for i in 0..x.rows() {
+        println!(
+            "{:>4}m, {:>2} -> {} | {}",
+            x[(i, 0)],
+            x[(i, 1)] as i32,
+            preds[(i, 0)] as i32,
+            bought[i] as i32
+        );
+    }
+
+    let new_visit = Matrix::from_rows_vec(vec![4.0, 4.0], 1, 2);
+    let p = model.predict_proba(&new_visit);
+    println!("Prob of purchase for 4min/4pages: {:.2}", p[(0, 0)]);
+}
+
+#[test]
+fn test_student_passing_example() {
+    let hours = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0];
+    let passed = vec![0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0];
+    let x = Matrix::from_vec(hours.clone(), hours.len(), 1);
+    let y = Matrix::from_vec(passed.clone(), passed.len(), 1);
+    let mut model = LogReg::new(1);
+    model.fit(&x, &y, 0.1, 10000);
+    let preds = model.predict(&x);
+    for i in 0..y.rows() {
+        assert_eq!(preds[(i, 0)], passed[i]);
+    }
+}
+
+#[test]
+fn test_purchase_prediction_example() {
+    let raw_x = vec![1.0, 2.0, 3.0, 1.0, 2.0, 4.0, 5.0, 5.0, 3.5, 2.0, 6.0, 6.0];
+    let bought = vec![0.0, 0.0, 0.0, 1.0, 0.0, 1.0];
+    let x = Matrix::from_rows_vec(raw_x, 6, 2);
+    let y = Matrix::from_vec(bought.clone(), bought.len(), 1);
+    let mut model = LogReg::new(2);
+    model.fit(&x, &y, 0.05, 20000);
+    let preds = model.predict(&x);
+    for i in 0..y.rows() {
+        assert_eq!(preds[(i, 0)], bought[i]);
+    }
+}
--- a/examples/pca.rs
+++ b/examples/pca.rs
@@ -0,0 +1,60 @@
+use rustframe::compute::models::pca::PCA;
+use rustframe::matrix::Matrix;
+
+/// Two dimensionality reduction examples using PCA.
+///
+/// Example 1 reduces 3D sensor readings to two components.
+/// Example 2 compresses a small four-feature dataset.
+fn main() {
+    sensor_demo();
+    println!("\n-----\n");
+    finance_demo();
+}
+
+fn sensor_demo() {
+    println!("Example 1: 3D sensor data");
+
+    // Ten 3D observations from an accelerometer
+    let raw = vec![
+        2.5, 2.4, 0.5, 0.5, 0.7, 1.5, 2.2, 2.9, 0.7, 1.9, 2.2, 1.0, 3.1, 3.0, 0.6, 2.3, 2.7, 0.9,
+        2.0, 1.6, 1.1, 1.0, 1.1, 1.9, 1.5, 1.6, 2.2, 1.1, 0.9, 2.1,
+    ];
+    let x = Matrix::from_rows_vec(raw, 10, 3);
+
+    let pca = PCA::fit(&x, 2, 0);
+    let reduced = pca.transform(&x);
+
+    println!("Components: {:?}", pca.components.data());
+    println!("First row -> {:.2?}", [reduced[(0, 0)], reduced[(0, 1)]]);
+}
+
+fn finance_demo() {
+    println!("Example 2: 4D finance data");
+
+    // Four daily percentage returns of different stocks
+    let raw = vec![
+        0.2, 0.1, -0.1, 0.0, 0.3, 0.2, -0.2, 0.1, 0.1, 0.0, -0.1, -0.1, 0.4, 0.3, -0.3, 0.2, 0.0,
+        -0.1, 0.1, -0.1,
+    ];
+    let x = Matrix::from_rows_vec(raw, 5, 4);
+
+    // Keep two principal components
+    let pca = PCA::fit(&x, 2, 0);
+    let reduced = pca.transform(&x);
+
+    println!("Reduced shape: {:?}", reduced.shape());
+    println!("First row -> {:.2?}", [reduced[(0, 0)], reduced[(0, 1)]]);
+}
+
+#[test]
+fn test_sensor_demo() {
+    let raw = vec![
+        2.5, 2.4, 0.5, 0.5, 0.7, 1.5, 2.2, 2.9, 0.7, 1.9, 2.2, 1.0, 3.1, 3.0, 0.6, 2.3, 2.7, 0.9,
+        2.0, 1.6, 1.1, 1.0, 1.1, 1.9, 1.5, 1.6, 2.2, 1.1, 0.9, 2.1,
+    ];
+    let x = Matrix::from_rows_vec(raw, 10, 3);
+    let pca = PCA::fit(&x, 2, 0);
+    let reduced = pca.transform(&x);
+    assert_eq!(reduced.rows(), 10);
+    assert_eq!(reduced.cols(), 2);
+}
--- a/examples/random_demo.rs
+++ b/examples/random_demo.rs
@@ -0,0 +1,67 @@
+use rustframe::random::{crypto_rng, rng, Rng, SliceRandom};
+
+/// Demonstrates basic usage of the random number generators.
+///
+/// It showcases uniform ranges, booleans, normal distribution,
+/// shuffling and the cryptographically secure generator.
+fn main() {
+    basic_usage();
+    println!("\n-----\n");
+    normal_demo();
+    println!("\n-----\n");
+    shuffle_demo();
+}
+
+fn basic_usage() {
+    println!("Basic PRNG usage\n----------------");
+    let mut prng = rng();
+    println!("random u64   : {}", prng.next_u64());
+    println!("range [10,20): {}", prng.random_range(10..20));
+    println!("bool         : {}", prng.gen_bool());
+}
+
+fn normal_demo() {
+    println!("Normal distribution\n-------------------");
+    let mut prng = rng();
+    for _ in 0..3 {
+        let v = prng.normal(0.0, 1.0);
+        println!("sample: {:.3}", v);
+    }
+}
+
+fn shuffle_demo() {
+    println!("Slice shuffling\n----------------");
+    let mut prng = rng();
+    let mut data = [1, 2, 3, 4, 5];
+    data.shuffle(&mut prng);
+    println!("shuffled: {:?}", data);
+
+    let mut secure = crypto_rng();
+    let byte = secure.random_range(0..256usize);
+    println!("crypto byte: {}", byte);
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use rustframe::random::{CryptoRng, Prng};
+
+    #[test]
+    fn test_basic_usage_range_bounds() {
+        let mut rng = Prng::new(1);
+        for _ in 0..50 {
+            let v = rng.random_range(5..10);
+            assert!(v >= 5 && v < 10);
+        }
+    }
+
+    #[test]
+    fn test_crypto_byte_bounds() {
+        let mut rng = CryptoRng::new();
+        for _ in 0..50 {
+            let v = rng.random_range(0..256usize);
+            assert!(v < 256);
+        }
+    }
+}
+
--- a/examples/random_stats.rs
+++ b/examples/random_stats.rs
@@ -0,0 +1,57 @@
+use rustframe::random::{crypto_rng, rng, Rng};
+
+/// Demonstrates simple statistical checks on random number generators.
+fn main() {
+    chi_square_demo();
+    println!("\n-----\n");
+    monobit_demo();
+}
+
+fn chi_square_demo() {
+    println!("Chi-square test on PRNG");
+    let mut rng = rng();
+    let mut counts = [0usize; 10];
+    let samples = 10000;
+    for _ in 0..samples {
+        let v = rng.random_range(0..10usize);
+        counts[v] += 1;
+    }
+    let expected = samples as f64 / 10.0;
+    let chi2: f64 = counts
+        .iter()
+        .map(|&c| {
+            let diff = c as f64 - expected;
+            diff * diff / expected
+        })
+        .sum();
+    println!("counts: {:?}", counts);
+    println!("chi-square: {:.3}", chi2);
+}
+
+fn monobit_demo() {
+    println!("Monobit test on crypto RNG");
+    let mut rng = crypto_rng();
+    let mut ones = 0usize;
+    let samples = 1000;
+    for _ in 0..samples {
+        ones += rng.next_u64().count_ones() as usize;
+    }
+    let ratio = ones as f64 / (samples as f64 * 64.0);
+    println!("ones ratio: {:.4}", ratio);
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_chi_square_demo_runs() {
+        chi_square_demo();
+    }
+
+    #[test]
+    fn test_monobit_demo_runs() {
+        monobit_demo();
+    }
+}
+
--- a/examples/stats_overview.rs
+++ b/examples/stats_overview.rs
@@ -0,0 +1,93 @@
+use rustframe::compute::stats::{
+    chi2_test, covariance, covariance_matrix, mean, median, pearson, percentile, stddev, t_test,
+};
+use rustframe::matrix::{Axis, Matrix};
+
+/// Demonstrates some of the statistics utilities in Rustframe.
+///
+/// The example is split into three parts:
+///   - Basic descriptive statistics on a small data set
+///   - Covariance and correlation calculations
+///   - Simple inferential tests (t-test and chi-square)
+fn main() {
+    descriptive_demo();
+    println!("\n-----\n");
+    correlation_demo();
+    println!("\n-----\n");
+    inferential_demo();
+}
+
+fn descriptive_demo() {
+    println!("Descriptive statistics\n----------------------");
+    let data = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0], 1, 5);
+    println!("mean           : {:.2}", mean(&data));
+    println!("std dev        : {:.2}", stddev(&data));
+    println!("median         : {:.2}", median(&data));
+    println!("25th percentile: {:.2}", percentile(&data, 25.0));
+}
+
+fn correlation_demo() {
+    println!("Covariance and Correlation\n--------------------------");
+    let x = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+    let y = Matrix::from_vec(vec![1.0, 2.0, 3.0, 5.0], 2, 2);
+    let cov = covariance(&x, &y);
+    let cov_mat = covariance_matrix(&x, Axis::Col);
+    let corr = pearson(&x, &y);
+    println!("covariance : {:.2}", cov);
+    println!("cov matrix : {:?}", cov_mat.data());
+    println!("pearson r  : {:.2}", corr);
+}
+
+fn inferential_demo() {
+    println!("Inferential statistics\n----------------------");
+    let s1 = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0], 1, 5);
+    let s2 = Matrix::from_vec(vec![6.0, 7.0, 8.0, 9.0, 10.0], 1, 5);
+    let (t_stat, t_p) = t_test(&s1, &s2);
+    println!("t statistic : {:.2}, p-value: {:.4}", t_stat, t_p);
+
+    let observed = Matrix::from_vec(vec![12.0, 5.0, 8.0, 10.0], 2, 2);
+    let (chi2, chi_p) = chi2_test(&observed);
+    println!("chi^2       : {:.2}, p-value: {:.4}", chi2, chi_p);
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    const EPS: f64 = 1e-8;
+
+    #[test]
+    fn test_descriptive_demo() {
+        let data = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0], 1, 5);
+        assert!((mean(&data) - 3.0).abs() < EPS);
+        assert!((stddev(&data) - 1.4142135623730951).abs() < EPS);
+        assert!((median(&data) - 3.0).abs() < EPS);
+        assert!((percentile(&data, 25.0) - 2.0).abs() < EPS);
+    }
+
+    #[test]
+    fn test_correlation_demo() {
+        let x = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+        let y = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 5.0], 2, 2);
+        let cov = covariance(&x, &y);
+        assert!((cov - 1.625).abs() < EPS);
+        let cov_mat = covariance_matrix(&x, Axis::Col);
+        assert!((cov_mat.get(0, 0) - 2.0).abs() < EPS);
+        assert!((cov_mat.get(1, 1) - 2.0).abs() < EPS);
+        let corr = pearson(&x, &y);
+        assert!((corr - 0.9827076298239908).abs() < 1e-6);
+    }
+
+    #[test]
+    fn test_inferential_demo() {
+        let s1 = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0], 1, 5);
+        let s2 = Matrix::from_rows_vec(vec![6.0, 7.0, 8.0, 9.0, 10.0], 1, 5);
+        let (t_stat, p_value) = t_test(&s1, &s2);
+        assert!((t_stat + 5.0).abs() < 1e-5);
+        assert!(p_value > 0.0 && p_value < 1.0);
+
+        let observed = Matrix::from_rows_vec(vec![12.0, 5.0, 8.0, 10.0], 2, 2);
+        let (chi2, p) = chi2_test(&observed);
+        assert!(chi2 > 0.0);
+        assert!(p > 0.0 && p < 1.0);
+    }
+}
--- a/src/compute/mod.rs
+++ b/src/compute/mod.rs
@@ -1,3 +1,16 @@
+//! Algorithms and statistical utilities built on top of the core matrices.
+//!
+//! This module groups together machine‑learning models and statistical helper
+//! functions. For quick access to basic statistics see [`stats`](crate::compute::stats), while
+//! [`models`](crate::compute::models) contains small learning algorithms.
+//!
+//! ```
+//! use rustframe::compute::stats;
+//! use rustframe::matrix::Matrix;
+//!
+//! let m = Matrix::from_vec(vec![1.0, 2.0, 3.0], 3, 1);
+//! assert_eq!(stats::mean(&m), 2.0);
+//! ```
 pub mod models;

 pub mod stats;
--- a/src/compute/models/activations.rs
+++ b/src/compute/models/activations.rs
@@ -1,3 +1,15 @@
+//! Common activation functions used in neural networks.
+//!
+//! Functions operate element-wise on [`Matrix`] values.
+//!
+//! ```
+//! use rustframe::compute::models::activations::sigmoid;
+//! use rustframe::matrix::Matrix;
+//!
+//! let x = Matrix::from_vec(vec![0.0], 1, 1);
+//! let y = sigmoid(&x);
+//! assert!((y.get(0,0) - 0.5).abs() < 1e-6);
+//! ```
 use crate::matrix::{Matrix, SeriesOps};

 pub fn sigmoid(x: &Matrix<f64>) -> Matrix<f64> {
@@ -25,6 +37,7 @@ pub fn dleaky_relu(x: &Matrix<f64>) -> Matrix<f64> {
    x.map(|v| if v > 0.0 { 1.0 } else { 0.01 })
 }

+#[cfg(test)]
 mod tests {
    use super::*;

--- a/src/compute/models/dense_nn.rs
+++ b/src/compute/models/dense_nn.rs
@@ -1,6 +1,33 @@
+//! A minimal dense neural network implementation for educational purposes.
+//!
+//! Layers operate on [`Matrix`] values and support ReLU and Sigmoid
+//! activations. This is not meant to be a performant deep‑learning framework
+//! but rather a small example of how the surrounding matrix utilities can be
+//! composed.
+//!
+//! ```
+//! use rustframe::compute::models::dense_nn::{ActivationKind, DenseNN, DenseNNConfig, InitializerKind, LossKind};
+//! use rustframe::matrix::Matrix;
+//!
+//! // Tiny network with one input and one output neuron.
+//! let config = DenseNNConfig {
+//!     input_size: 1,
+//!     hidden_layers: vec![],
+//!     output_size: 1,
+//!     activations: vec![ActivationKind::Relu],
+//!     initializer: InitializerKind::Uniform(0.5),
+//!     loss: LossKind::MSE,
+//!     learning_rate: 0.1,
+//!     epochs: 1,
+//! };
+//! let mut nn = DenseNN::new(config);
+//! let x = Matrix::from_vec(vec![1.0, 2.0], 2, 1);
+//! let y = Matrix::from_vec(vec![2.0, 3.0], 2, 1);
+//! nn.train(&x, &y);
+//! ```
 use crate::compute::models::activations::{drelu, relu, sigmoid};
 use crate::matrix::{Matrix, SeriesOps};
-use rand::prelude::*;
+use crate::random::prelude::*;

 /// Supported activation functions
 #[derive(Clone)]
@@ -46,7 +73,7 @@ pub enum InitializerKind {

 impl InitializerKind {
    pub fn initialize(&self, rows: usize, cols: usize) -> Matrix<f64> {
-        let mut rng = rand::rng();
+        let mut rng = rng();
        let fan_in = rows;
        let fan_out = cols;
        let limit = match self {
--- a/src/compute/models/gaussian_nb.rs
+++ b/src/compute/models/gaussian_nb.rs
@@ -1,3 +1,16 @@
+//! Gaussian Naive Bayes classifier for dense matrices.
+//!
+//! ```
+//! use rustframe::compute::models::gaussian_nb::GaussianNB;
+//! use rustframe::matrix::Matrix;
+//!
+//! let x = Matrix::from_vec(vec![1.0, 2.0, 1.0, 2.0], 2, 2); // two samples
+//! let y = Matrix::from_vec(vec![0.0, 1.0], 2, 1);
+//! let mut model = GaussianNB::new(1e-9, false);
+//! model.fit(&x, &y);
+//! let preds = model.predict(&x);
+//! assert_eq!(preds.rows(), 2);
+//! ```
 use crate::matrix::Matrix;
 use std::collections::HashMap;

--- a/src/compute/models/k_means.rs
+++ b/src/compute/models/k_means.rs
@@ -1,7 +1,17 @@
+//! Simple k-means clustering working on [`Matrix`] data.
+//!
+//! ```
+//! use rustframe::compute::models::k_means::KMeans;
+//! use rustframe::matrix::Matrix;
+//!
+//! let data = Matrix::from_vec(vec![1.0, 1.0, 5.0, 5.0], 2, 2);
+//! let (model, labels) = KMeans::fit(&data, 2, 10, 1e-4);
+//! assert_eq!(model.centroids.rows(), 2);
+//! assert_eq!(labels.len(), 2);
+//! ```
 use crate::compute::stats::mean_vertical;
 use crate::matrix::Matrix;
-use rand::rng;
-use rand::seq::SliceRandom;
+use crate::random::prelude::*;

 pub struct KMeans {
    pub centroids: Matrix<f64>, // (k, n_features)
@@ -193,7 +203,8 @@ mod tests {
                                break;
                            }
                        }
-                        assert!(matches_data_point, "Centroid {} (empty cluster) does not match any data point", c);
+                        // "Centroid {} (empty cluster) does not match any data point",c
+                        assert!(matches_data_point);
                    }
                }
                break;
@@ -360,5 +371,4 @@ mod tests {
        assert_eq!(predicted_label.len(), 1);
        assert!(predicted_label[0] < k);
    }
-
 }
--- a/src/compute/models/linreg.rs
+++ b/src/compute/models/linreg.rs
@@ -1,3 +1,16 @@
+//! Ordinary least squares linear regression.
+//!
+//! ```
+//! use rustframe::compute::models::linreg::LinReg;
+//! use rustframe::matrix::Matrix;
+//!
+//! let x = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 4, 1);
+//! let y = Matrix::from_vec(vec![2.0, 3.0, 4.0, 5.0], 4, 1);
+//! let mut model = LinReg::new(1);
+//! model.fit(&x, &y, 0.01, 100);
+//! let preds = model.predict(&x);
+//! assert_eq!(preds.rows(), 4);
+//! ```
 use crate::matrix::{Matrix, SeriesOps};

 pub struct LinReg {
--- a/src/compute/models/logreg.rs
+++ b/src/compute/models/logreg.rs
@@ -1,3 +1,16 @@
+//! Binary logistic regression classifier.
+//!
+//! ```
+//! use rustframe::compute::models::logreg::LogReg;
+//! use rustframe::matrix::Matrix;
+//!
+//! let x = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 4, 1);
+//! let y = Matrix::from_vec(vec![0.0, 0.0, 1.0, 1.0], 4, 1);
+//! let mut model = LogReg::new(1);
+//! model.fit(&x, &y, 0.1, 100);
+//! let preds = model.predict(&x);
+//! assert_eq!(preds[(0,0)], 0.0);
+//! ```
 use crate::compute::models::activations::sigmoid;
 use crate::matrix::{Matrix, SeriesOps};

--- a/src/compute/models/mod.rs
+++ b/src/compute/models/mod.rs
@@ -1,3 +1,19 @@
+//! Lightweight machine‑learning models built on matrices.
+//!
+//! Models are intentionally minimal and operate on the [`Matrix`](crate::matrix::Matrix) type for
+//! inputs and parameters.
+//!
+//! ```
+//! use rustframe::compute::models::linreg::LinReg;
+//! use rustframe::matrix::Matrix;
+//!
+//! let x = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 4, 1);
+//! let y = Matrix::from_vec(vec![2.0, 3.0, 4.0, 5.0], 4, 1);
+//! let mut model = LinReg::new(1);
+//! model.fit(&x, &y, 0.01, 1000);
+//! let preds = model.predict(&x);
+//! assert_eq!(preds.rows(), 4);
+//! ```
 pub mod activations;
 pub mod dense_nn;
 pub mod gaussian_nb;
--- a/src/compute/models/pca.rs
+++ b/src/compute/models/pca.rs
@@ -1,3 +1,14 @@
+//! Principal Component Analysis using covariance matrices.
+//!
+//! ```
+//! use rustframe::compute::models::pca::PCA;
+//! use rustframe::matrix::Matrix;
+//!
+//! let data = Matrix::from_rows_vec(vec![1.0, 1.0, 2.0, 2.0], 2, 2);
+//! let pca = PCA::fit(&data, 1, 0);
+//! let projected = pca.transform(&data);
+//! assert_eq!(projected.cols(), 1);
+//! ```
 use crate::compute::stats::correlation::covariance_matrix;
 use crate::compute::stats::descriptive::mean_vertical;
 use crate::matrix::{Axis, Matrix, SeriesOps};
@@ -44,11 +55,7 @@ mod tests {

    #[test]
    fn test_pca_basic() {
-        // Simple 2D data, points along y=x line
-        // Data:
-        // 1.0, 1.0
-        // 2.0, 2.0
-        // 3.0, 3.0
+        // Simple 2D data with points along the y = x line
        let data = Matrix::from_rows_vec(vec![1.0, 1.0, 2.0, 2.0, 3.0, 3.0], 3, 2);
        let (_n_samples, _n_features) = data.shape();

@@ -71,15 +78,7 @@ mod tests {
        assert!((pca.components.get(0, 0) - 1.0).abs() < EPSILON);
        assert!((pca.components.get(0, 1) - 1.0).abs() < EPSILON);

-        // Test transform
-        // Centered data:
-        // -1.0, -1.0
-        //  0.0,  0.0
-        //  1.0,  1.0
-        // Projected: (centered_data * components.transpose())
-        // (-1.0 * 1.0 + -1.0 * 1.0) = -2.0
-        // ( 0.0 * 1.0 +  0.0 * 1.0) =  0.0
-        // ( 1.0 * 1.0 +  1.0 * 1.0) =  2.0
+        // Test transform: centered data projects to [-2.0, 0.0, 2.0]
        let transformed_data = pca.transform(&data);
        assert_eq!(transformed_data.rows(), 3);
        assert_eq!(transformed_data.cols(), 1);
--- a/src/compute/stats/correlation.rs
+++ b/src/compute/stats/correlation.rs
@@ -1,3 +1,16 @@
+//! Covariance and correlation helpers.
+//!
+//! This module provides routines for measuring the relationship between
+//! columns or rows of matrices.
+//!
+//! ```
+//! use rustframe::compute::stats::correlation;
+//! use rustframe::matrix::Matrix;
+//!
+//! let x = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+//! let cov = correlation::covariance(&x, &x);
+//! assert!((cov - 1.25).abs() < 1e-8);
+//! ```
 use crate::compute::stats::{mean, mean_horizontal, mean_vertical, stddev};
 use crate::matrix::{Axis, Matrix, SeriesOps};

@@ -137,10 +150,7 @@ mod tests {

    #[test]
    fn test_covariance_scalar_same_matrix() {
-        // M =
-        // 1,2
-        // 3,4
-        // mean = 2.5
+        // Matrix with rows [1, 2] and [3, 4]; mean is 2.5
        let data = vec![1.0, 2.0, 3.0, 4.0];
        let m = Matrix::from_vec(data.clone(), 2, 2);

@@ -152,10 +162,7 @@ mod tests {

    #[test]
    fn test_covariance_scalar_diff_matrix() {
-        // x =
-        // 1,2
-        // 3,4
-        // y = 2*x
+        // Matrix x has rows [1, 2] and [3, 4]; y is two times x
        let x = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
        let y = Matrix::from_vec(vec![2.0, 4.0, 6.0, 8.0], 2, 2);

@@ -167,10 +174,7 @@ mod tests {

    #[test]
    fn test_covariance_vertical() {
-        // M =
-        // 1,2
-        // 3,4
-        // cols are [1,3] and [2,4], each var=1, cov=1
+        // Matrix with rows [1, 2] and [3, 4]; columns are [1,3] and [2,4], each var=1, cov=1
        let m = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
        let cov_mat = covariance_vertical(&m);

@@ -184,10 +188,7 @@ mod tests {

    #[test]
    fn test_covariance_horizontal() {
-        // M =
-        // 1,2
-        // 3,4
-        // rows are [1,2] and [3,4], each var=0.25, cov=0.25
+        // Matrix with rows [1,2] and [3,4], each var=0.25, cov=0.25
        let m = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
        let cov_mat = covariance_horizontal(&m);

@@ -201,10 +202,7 @@ mod tests {

    #[test]
    fn test_covariance_matrix_vertical() {
-        // Test with a simple 2x2 matrix
-        // M =
-        // 1, 2
-        // 3, 4
+        // Test with a simple 2x2 matrix with rows [1, 2] and [3, 4]
        // Expected covariance matrix (vertical, i.e., between columns):
        // Col1: [1, 3], mean = 2
        // Col2: [2, 4], mean = 3
@@ -212,9 +210,7 @@ mod tests {
        // Cov(Col2, Col2) = ((2-3)^2 + (4-3)^2) / (2-1) = (1+1)/1 = 2
        // Cov(Col1, Col2) = ((1-2)*(2-3) + (3-2)*(4-3)) / (2-1) = ((-1)*(-1) + (1)*(1))/1 = (1+1)/1 = 2
        // Cov(Col2, Col1) = 2
-        // Expected:
-        // 2, 2
-        // 2, 2
+        // Expected matrix filled with 2
        let m = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
        let cov_mat = covariance_matrix(&m, Axis::Col);

@@ -226,10 +222,7 @@ mod tests {

    #[test]
    fn test_covariance_matrix_horizontal() {
-        // Test with a simple 2x2 matrix
-        // M =
-        // 1, 2
-        // 3, 4
+        // Test with a simple 2x2 matrix with rows [1, 2] and [3, 4]
        // Expected covariance matrix (horizontal, i.e., between rows):
        // Row1: [1, 2], mean = 1.5
        // Row2: [3, 4], mean = 3.5
@@ -237,9 +230,7 @@ mod tests {
        // Cov(Row2, Row2) = ((3-3.5)^2 + (4-3.5)^2) / (2-1) = (0.25+0.25)/1 = 0.5
        // Cov(Row1, Row2) = ((1-1.5)*(3-3.5) + (2-1.5)*(4-3.5)) / (2-1) = ((-0.5)*(-0.5) + (0.5)*(0.5))/1 = (0.25+0.25)/1 = 0.5
        // Cov(Row2, Row1) = 0.5
-        // Expected:
-        // 0.5, -0.5
-        // -0.5, 0.5
+        // Expected matrix: [[0.5, -0.5], [-0.5, 0.5]]
        let m = Matrix::from_rows_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
        let cov_mat = covariance_matrix(&m, Axis::Row);

--- a/src/compute/stats/descriptive.rs
+++ b/src/compute/stats/descriptive.rs
@@ -1,3 +1,15 @@
+//! Descriptive statistics for matrices.
+//!
+//! Provides means, variances, medians and other aggregations computed either
+//! across the whole matrix or along a specific axis.
+//!
+//! ```
+//! use rustframe::compute::stats::descriptive;
+//! use rustframe::matrix::Matrix;
+//!
+//! let m = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+//! assert_eq!(descriptive::mean(&m), 2.5);
+//! ```
 use crate::matrix::{Axis, Matrix, SeriesOps};

 pub fn mean(x: &Matrix<f64>) -> f64 {
@@ -350,11 +362,7 @@ mod tests {
        let data: Vec<f64> = (1..=24).map(|x| x as f64).collect();
        let x = Matrix::from_vec(data, 4, 6);

-        // columns:
-        //  1,  5,  9, 13, 17, 21
-        //  2,  6, 10, 14, 18, 22
-        //  3,  7, 11, 15, 19, 23
-        //  4,  8, 12, 16, 20, 24
+        // columns contain sequences increasing by four starting at 1 through 4

        let er0 = vec![1., 5., 9., 13., 17., 21.];
        let er50 = vec![3., 7., 11., 15., 19., 23.];
--- a/src/compute/stats/distributions.rs
+++ b/src/compute/stats/distributions.rs
@@ -1,3 +1,16 @@
+//! Probability distribution functions applied element-wise to matrices.
+//!
+//! Includes approximations for the normal, uniform and gamma distributions as
+//! well as the error function.
+//!
+//! ```
+//! use rustframe::compute::stats::distributions;
+//! use rustframe::matrix::Matrix;
+//!
+//! let x = Matrix::from_vec(vec![0.0], 1, 1);
+//! let pdf = distributions::normal_pdf(x.clone(), 0.0, 1.0);
+//! assert!((pdf.get(0,0) - 0.3989).abs() < 1e-3);
+//! ```
 use crate::matrix::{Matrix, SeriesOps};

 use std::f64::consts::PI;
--- a/src/compute/stats/inferential.rs
+++ b/src/compute/stats/inferential.rs
@@ -1,3 +1,14 @@
+//! Basic inferential statistics such as t‑tests and chi‑square tests.
+//!
+//! ```
+//! use rustframe::compute::stats::inferential;
+//! use rustframe::matrix::Matrix;
+//!
+//! let a = Matrix::from_vec(vec![1.0, 2.0], 2, 1);
+//! let b = Matrix::from_vec(vec![1.1, 1.9], 2, 1);
+//! let (t, _p) = inferential::t_test(&a, &b);
+//! assert!(t.abs() < 1.0);
+//! ```
 use crate::matrix::{Matrix, SeriesOps};

 use crate::compute::stats::{gamma_cdf, mean, sample_variance};
--- a/src/compute/stats/mod.rs
+++ b/src/compute/stats/mod.rs
@@ -1,3 +1,16 @@
+//! Statistical routines for matrices.
+//!
+//! Functions are grouped into submodules for descriptive statistics,
+//! correlations, probability distributions and basic inferential tests.
+//!
+//! ```
+//! use rustframe::compute::stats;
+//! use rustframe::matrix::Matrix;
+//!
+//! let m = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 2, 2);
+//! let cov = stats::covariance(&m, &m);
+//! assert!((cov - 1.25).abs() < 1e-8);
+//! ```
 pub mod correlation;
 pub mod descriptive;
 pub mod distributions;
--- a/src/frame/base.rs
+++ b/src/frame/base.rs
@@ -1,3 +1,19 @@
+//! Core data-frame structures such as [`Frame`] and [`RowIndex`].
+//!
+//! The [`Frame`] type stores column-labelled data with an optional row index
+//! and builds upon the [`crate::matrix::Matrix`] type.
+//!
+//! # Examples
+//!
+//! ```
+//! use rustframe::frame::{Frame, RowIndex};
+//! use rustframe::matrix::Matrix;
+//!
+//! let data = Matrix::from_cols(vec![vec![1, 2], vec![3, 4]]);
+//! let frame = Frame::new(data, vec!["L", "R"], Some(RowIndex::Int(vec![10, 20])));
+//! assert_eq!(frame.columns(), &["L", "R"]);
+//! assert_eq!(frame.index(), &RowIndex::Int(vec![10, 20]));
+//! ```
 use crate::matrix::Matrix;
 use chrono::NaiveDate;
 use std::collections::HashMap;
--- a/src/frame/mod.rs
+++ b/src/frame/mod.rs
@@ -1,3 +1,21 @@
+//! High-level interface for working with columnar data and row indices.
+//!
+//! The [`Frame`](crate::frame::Frame) type combines a matrix with column labels and a typed row
+//! index, similar to data frames in other data-analysis libraries.
+//!
+//! # Examples
+//!
+//! ```
+//! use rustframe::frame::{Frame, RowIndex};
+//! use rustframe::matrix::Matrix;
+//!
+//! // Build a frame from two columns labelled "A" and "B".
+//! let data = Matrix::from_cols(vec![vec![1.0, 2.0], vec![3.0, 4.0]]);
+//! let frame = Frame::new(data, vec!["A", "B"], None);
+//!
+//! assert_eq!(frame["A"], vec![1.0, 2.0]);
+//! assert_eq!(frame.index(), &RowIndex::Range(0..2));
+//! ```
 pub mod base;
 pub mod ops;

--- a/src/frame/ops.rs
+++ b/src/frame/ops.rs
@@ -1,3 +1,16 @@
+//! Trait implementations that allow [`Frame`] to reuse matrix operations.
+//!
+//! These modules forward numeric and boolean aggregation methods from the
+//! underlying [`Matrix`](crate::matrix::Matrix) type so that they can be called
+//! directly on a [`Frame`].
+//!
+//! ```
+//! use rustframe::frame::Frame;
+//! use rustframe::matrix::{Matrix, SeriesOps};
+//!
+//! let frame = Frame::new(Matrix::from_cols(vec![vec![1.0, 2.0]]), vec!["A"], None);
+//! assert_eq!(frame.sum_vertical(), vec![3.0]);
+//! ```
 use crate::frame::Frame;
 use crate::matrix::{Axis, BoolMatrix, BoolOps, FloatMatrix, SeriesOps};

--- a/src/lib.rs
+++ b/src/lib.rs
@@ -14,3 +14,6 @@ pub mod utils;

 /// Documentation for the [`crate::compute`] module.
 pub mod compute;
+
+/// Documentation for the [`crate::random`] module.
+pub mod random;
--- a/src/matrix/boolops.rs
+++ b/src/matrix/boolops.rs
@@ -1,3 +1,14 @@
+//! Logical reductions for boolean matrices.
+//!
+//! The [`BoolOps`] trait mirrors common boolean aggregations such as `any` and
+//! `all` over rows or columns of a [`BoolMatrix`].
+//!
+//! ```
+//! use rustframe::matrix::{BoolMatrix, BoolOps};
+//!
+//! let m = BoolMatrix::from_vec(vec![true, false], 2, 1);
+//! assert!(m.any());
+//! ```
 use crate::matrix::{Axis, BoolMatrix};

 /// Boolean operations on `Matrix<bool>`
--- a/src/matrix/mat.rs
+++ b/src/matrix/mat.rs
@@ -1028,9 +1028,7 @@ mod tests {

    #[test]
    fn test_from_rows_vec() {
-        // Representing:
-        // 1 2 3
-        // 4 5 6
+        // Matrix with rows [1, 2, 3] and [4, 5, 6]
        let rows_data = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
        let matrix = Matrix::from_rows_vec(rows_data, 2, 3);

@@ -1042,19 +1040,14 @@ mod tests {

    // Helper function to create a basic Matrix for testing
    fn static_test_matrix() -> Matrix<i32> {
-        // Column-major data:
-        // 1 4 7
-        // 2 5 8
-        // 3 6 9
+        // Column-major data representing a 3x3 matrix of sequential integers
        let data = vec![1, 2, 3, 4, 5, 6, 7, 8, 9];
        Matrix::from_vec(data, 3, 3)
    }

    // Another helper for a different size
    fn static_test_matrix_2x4() -> Matrix<i32> {
-        // Column-major data:
-        // 1 3 5 7
-        // 2 4 6 8
+        // Column-major data representing a 2x4 matrix of sequential integers
        let data = vec![1, 2, 3, 4, 5, 6, 7, 8];
        Matrix::from_vec(data, 2, 4)
    }
@@ -1132,10 +1125,7 @@ mod tests {

    #[test]
    fn test_from_cols_basic() {
-        // Representing:
-        // 1 4 7
-        // 2 5 8
-        // 3 6 9
+        // Matrix with columns forming a 3x3 sequence
        let cols_data = vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]];
        let matrix = Matrix::from_cols(cols_data);

@@ -1512,8 +1502,7 @@ mod tests {

        // Delete the first row
        matrix.delete_row(0);
-        // Should be:
-        // 3 6 9
+        // Resulting data should be [3, 6, 9]
        assert_eq!(matrix.rows(), 1);
        assert_eq!(matrix.cols(), 3);
        assert_eq!(matrix.data(), &[3, 6, 9]);
--- a/src/matrix/mod.rs
+++ b/src/matrix/mod.rs
@@ -1,3 +1,18 @@
+//! Core matrix types and operations.
+//!
+//! The [`Matrix`](crate::matrix::Matrix) struct provides a simple column‑major 2D array with a
+//! suite of numeric helpers. Additional traits like [`SeriesOps`](crate::matrix::SeriesOps) and
+//! [`BoolOps`](crate::matrix::BoolOps) extend functionality for common statistics and logical reductions.
+//!
+//! # Examples
+//!
+//! ```
+//! use rustframe::matrix::Matrix;
+//!
+//! let m = Matrix::from_cols(vec![vec![1, 2], vec![3, 4]]);
+//! assert_eq!(m.shape(), (2, 2));
+//! assert_eq!(m[(0,1)], 3);
+//! ```
 pub mod boolops;
 pub mod mat;
 pub mod seriesops;
--- a/src/matrix/seriesops.rs
+++ b/src/matrix/seriesops.rs
@@ -1,3 +1,14 @@
+//! Numeric reductions and transformations over matrix axes.
+//!
+//! [`SeriesOps`] provides methods like [`SeriesOps::sum_vertical`] or
+//! [`SeriesOps::map`] that operate on [`FloatMatrix`] values.
+//!
+//! ```
+//! use rustframe::matrix::{Matrix, SeriesOps};
+//!
+//! let m = Matrix::from_cols(vec![vec![1.0, 2.0], vec![3.0, 4.0]]);
+//! assert_eq!(m.sum_horizontal(), vec![4.0, 6.0]);
+//! ```
 use crate::matrix::{Axis, BoolMatrix, FloatMatrix};

 /// "Series-like" helpers that work along a single axis.
@@ -215,20 +226,13 @@ mod tests {

    // Helper function to create a FloatMatrix for SeriesOps testing
    fn create_float_test_matrix() -> FloatMatrix {
-        // 3x3 matrix (column-major) with some NaNs
-        // 1.0  4.0  7.0
-        // 2.0  NaN  8.0
-        // 3.0  6.0  NaN
+        // 3x3 column-major matrix containing a few NaN values
        let data = vec![1.0, 2.0, 3.0, 4.0, f64::NAN, 6.0, 7.0, 8.0, f64::NAN];
        FloatMatrix::from_vec(data, 3, 3)
    }

    fn create_float_test_matrix_4x4() -> FloatMatrix {
-        // 4x4 matrix (column-major) with some NaNs
-        // 1.0  5.0  9.0  13.0
-        // 2.0  NaN  10.0 NaN
-        // 3.0  6.0  NaN  14.0
-        // NaN  7.0  11.0 NaN
+        // 4x4 column-major matrix with NaNs inserted at positions where index % 5 == 0
        // first make array with 16 elements
        FloatMatrix::from_vec(
            (0..16)
--- a/src/random/crypto.rs
+++ b/src/random/crypto.rs
@@ -0,0 +1,237 @@
+//! Cryptographically secure random number generator.
+//!
+//! On Unix systems this reads from `/dev/urandom`; on Windows it uses the
+//! system's preferred CNG provider.
+//!
+//! ```
+//! use rustframe::random::{crypto_rng, Rng};
+//! let mut rng = crypto_rng();
+//! let _v = rng.next_u64();
+//! ```
+#[cfg(unix)]
+use std::{fs::File, io::Read};
+
+use crate::random::Rng;
+
+#[cfg(unix)]
+pub struct CryptoRng {
+    file: File,
+}
+
+#[cfg(unix)]
+impl CryptoRng {
+    /// Open `/dev/urandom`.
+    pub fn new() -> Self {
+        let file = File::open("/dev/urandom").expect("failed to open /dev/urandom");
+        Self { file }
+    }
+}
+
+#[cfg(unix)]
+impl Rng for CryptoRng {
+    fn next_u64(&mut self) -> u64 {
+        let mut buf = [0u8; 8];
+        self.file
+            .read_exact(&mut buf)
+            .expect("failed reading from /dev/urandom");
+        u64::from_ne_bytes(buf)
+    }
+}
+
+#[cfg(windows)]
+pub struct CryptoRng;
+
+#[cfg(windows)]
+impl CryptoRng {
+    /// No handle is needed on Windows.
+    pub fn new() -> Self {
+        Self
+    }
+}
+
+#[cfg(windows)]
+impl Rng for CryptoRng {
+    fn next_u64(&mut self) -> u64 {
+        let mut buf = [0u8; 8];
+        win_fill(&mut buf).expect("BCryptGenRandom failed");
+        u64::from_ne_bytes(buf)
+    }
+}
+
+/// Fill `buf` with cryptographically secure random bytes using CNG.
+///
+/// * `BCryptGenRandom(NULL, buf, len, BCRYPT_USE_SYSTEM_PREFERRED_RNG)`
+///   asks the OS for its system‑preferred DRBG (CTR_DRBG on modern
+///   Windows).
+#[cfg(windows)]
+fn win_fill(buf: &mut [u8]) -> Result<(), ()> {
+    use core::ffi::c_void;
+
+    type BcryptAlgHandle = *mut c_void;
+    type NTSTATUS = i32;
+
+    const BCRYPT_USE_SYSTEM_PREFERRED_RNG: u32 = 0x0000_0002;
+
+    #[link(name = "bcrypt")]
+    extern "system" {
+        fn BCryptGenRandom(
+            hAlgorithm: BcryptAlgHandle,
+            pbBuffer: *mut u8,
+            cbBuffer: u32,
+            dwFlags: u32,
+        ) -> NTSTATUS;
+    }
+
+    // NT_SUCCESS(status) == status >= 0
+    let status = unsafe {
+        BCryptGenRandom(
+            core::ptr::null_mut(),
+            buf.as_mut_ptr(),
+            buf.len() as u32,
+            BCRYPT_USE_SYSTEM_PREFERRED_RNG,
+        )
+    };
+
+    if status >= 0 {
+        Ok(())
+    } else {
+        Err(())
+    }
+}
+
+/// Convenience constructor for [`CryptoRng`].
+pub fn crypto_rng() -> CryptoRng {
+    CryptoRng::new()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::random::Rng;
+    use std::collections::HashSet;
+
+    #[test]
+    fn test_crypto_rng_nonzero() {
+        let mut rng = CryptoRng::new();
+        let mut all_same = true;
+        let mut prev = rng.next_u64();
+        for _ in 0..5 {
+            let val = rng.next_u64();
+            if val != prev {
+                all_same = false;
+            }
+            prev = val;
+        }
+        assert!(!all_same, "CryptoRng produced identical values");
+    }
+
+    #[test]
+    fn test_crypto_rng_variation_large() {
+        let mut rng = CryptoRng::new();
+        let mut values = HashSet::new();
+        for _ in 0..100 {
+            values.insert(rng.next_u64());
+        }
+        assert!(values.len() > 90, "CryptoRng output not varied enough");
+    }
+
+    #[test]
+    fn test_crypto_rng_random_range_uniform() {
+        let mut rng = CryptoRng::new();
+        let mut counts = [0usize; 10];
+        for _ in 0..1000 {
+            let v = rng.random_range(0..10usize);
+            counts[v] += 1;
+        }
+        for &c in &counts {
+            // "Crypto RNG counts far from uniform: {c}"
+            assert!((c as isize - 100).abs() < 50);
+        }
+    }
+
+    #[test]
+    fn test_crypto_normal_distribution() {
+        let mut rng = CryptoRng::new();
+        let mean = 0.0;
+        let sd = 1.0;
+        let n = 2000;
+        let mut sum = 0.0;
+        let mut sum_sq = 0.0;
+        for _ in 0..n {
+            let val = rng.normal(mean, sd);
+            sum += val;
+            sum_sq += val * val;
+        }
+        let sample_mean = sum / n as f64;
+        let sample_var = sum_sq / n as f64 - sample_mean * sample_mean;
+        assert!(sample_mean.abs() < 0.1);
+        assert!((sample_var - 1.0).abs() < 0.2);
+    }
+
+    #[test]
+    fn test_two_instances_different_values() {
+        let mut a = CryptoRng::new();
+        let mut b = CryptoRng::new();
+        let va = a.next_u64();
+        let vb = b.next_u64();
+        assert_ne!(va, vb);
+    }
+
+    #[test]
+    fn test_crypto_rng_helper_function() {
+        let mut rng = crypto_rng();
+        let _ = rng.next_u64();
+    }
+
+    #[test]
+    fn test_crypto_normal_zero_sd() {
+        let mut rng = CryptoRng::new();
+        for _ in 0..5 {
+            let v = rng.normal(10.0, 0.0);
+            assert_eq!(v, 10.0);
+        }
+    }
+
+    #[test]
+    fn test_crypto_shuffle_empty_slice() {
+        use crate::random::SliceRandom;
+        let mut rng = CryptoRng::new();
+        let mut arr: [u8; 0] = [];
+        arr.shuffle(&mut rng);
+        assert!(arr.is_empty());
+    }
+
+    #[test]
+    fn test_crypto_chi_square_uniform() {
+        let mut rng = CryptoRng::new();
+        let mut counts = [0usize; 10];
+        let samples = 10000;
+        for _ in 0..samples {
+            let v = rng.random_range(0..10usize);
+            counts[v] += 1;
+        }
+        let expected = samples as f64 / 10.0;
+        let chi2: f64 = counts
+            .iter()
+            .map(|&c| {
+                let diff = c as f64 - expected;
+                diff * diff / expected
+            })
+            .sum();
+        assert!(chi2 < 40.0, "chi-square statistic too high: {chi2}");
+    }
+
+    #[test]
+    fn test_crypto_monobit() {
+        let mut rng = CryptoRng::new();
+        let mut ones = 0usize;
+        let samples = 1000;
+        for _ in 0..samples {
+            ones += rng.next_u64().count_ones() as usize;
+        }
+        let total_bits = samples * 64;
+        let ratio = ones as f64 / total_bits as f64;
+        // "bit ratio far from 0.5: {ratio}"
+        assert!((ratio - 0.5).abs() < 0.02);
+    }
+}
--- a/src/random/mod.rs
+++ b/src/random/mod.rs
@@ -0,0 +1,29 @@
+//! Random number generation utilities.
+//!
+//! Provides both a simple pseudo-random generator [`Prng`](crate::random::Prng) and a
+//! cryptographically secure alternative [`CryptoRng`](crate::random::CryptoRng). The
+//! [`SliceRandom`](crate::random::SliceRandom) trait offers shuffling of slices using any RNG
+//! implementing [`Rng`](crate::random::Rng).
+//!
+//! ```
+//! use rustframe::random::{rng, SliceRandom};
+//!
+//! let mut rng = rng();
+//! let mut data = [1, 2, 3, 4];
+//! data.shuffle(&mut rng);
+//! assert_eq!(data.len(), 4);
+//! ```
+pub mod crypto;
+pub mod prng;
+pub mod random_core;
+pub mod seq;
+
+pub use crypto::{crypto_rng, CryptoRng};
+pub use prng::{rng, Prng};
+pub use random_core::{RangeSample, Rng};
+pub use seq::SliceRandom;
+
+pub mod prelude {
+    pub use super::seq::SliceRandom;
+    pub use super::{crypto_rng, rng, CryptoRng, Prng, RangeSample, Rng};
+}
--- a/src/random/prng.rs
+++ b/src/random/prng.rs
@@ -0,0 +1,235 @@
+//! A tiny XorShift64-based pseudo random number generator.
+//!
+//! ```
+//! use rustframe::random::{rng, Rng};
+//! let mut rng = rng();
+//! let x = rng.next_u64();
+//! assert!(x >= 0);
+//! ```
+use std::time::{SystemTime, UNIX_EPOCH};
+
+use crate::random::Rng;
+
+/// Simple XorShift64-based pseudo random number generator.
+#[derive(Clone)]
+pub struct Prng {
+    state: u64,
+}
+
+impl Prng {
+    /// Create a new generator from the given seed.
+    pub fn new(seed: u64) -> Self {
+        Self { state: seed }
+    }
+
+    /// Create a generator seeded from the current time.
+    pub fn from_entropy() -> Self {
+        let nanos = SystemTime::now()
+            .duration_since(UNIX_EPOCH)
+            .unwrap()
+            .as_nanos() as u64;
+        Self::new(nanos)
+    }
+}
+
+impl Rng for Prng {
+    fn next_u64(&mut self) -> u64 {
+        let mut x = self.state;
+        x ^= x << 13;
+        x ^= x >> 7;
+        x ^= x << 17;
+        self.state = x;
+        x
+    }
+}
+
+/// Convenience constructor using system entropy.
+pub fn rng() -> Prng {
+    Prng::from_entropy()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::random::Rng;
+
+    #[test]
+    fn test_prng_determinism() {
+        let mut a = Prng::new(42);
+        let mut b = Prng::new(42);
+        for _ in 0..5 {
+            assert_eq!(a.next_u64(), b.next_u64());
+        }
+    }
+
+    #[test]
+    fn test_random_range_f64() {
+        let mut rng = Prng::new(1);
+        for _ in 0..10 {
+            let v = rng.random_range(-1.0..1.0);
+            assert!(v >= -1.0 && v < 1.0);
+        }
+    }
+
+    #[test]
+    fn test_random_range_usize() {
+        let mut rng = Prng::new(9);
+        for _ in 0..100 {
+            let v = rng.random_range(10..20);
+            assert!(v >= 10 && v < 20);
+        }
+    }
+
+    #[test]
+    fn test_gen_bool_balance() {
+        let mut rng = Prng::new(123);
+        let mut trues = 0;
+        for _ in 0..1000 {
+            if rng.gen_bool() {
+                trues += 1;
+            }
+        }
+        let ratio = trues as f64 / 1000.0;
+        assert!(ratio > 0.4 && ratio < 0.6);
+    }
+
+    #[test]
+    fn test_normal_distribution() {
+        let mut rng = Prng::new(7);
+        let mut sum = 0.0;
+        let mut sum_sq = 0.0;
+        let mean = 5.0;
+        let sd = 2.0;
+        let n = 5000;
+        for _ in 0..n {
+            let val = rng.normal(mean, sd);
+            sum += val;
+            sum_sq += val * val;
+        }
+        let sample_mean = sum / n as f64;
+        let sample_var = sum_sq / n as f64 - sample_mean * sample_mean;
+        assert!((sample_mean - mean).abs() < 0.1);
+        assert!((sample_var - sd * sd).abs() < 0.2 * sd * sd);
+    }
+
+    #[test]
+    fn test_prng_from_entropy_unique() {
+        use std::{collections::HashSet, thread, time::Duration};
+        let mut seen = HashSet::new();
+        for _ in 0..5 {
+            let mut rng = Prng::from_entropy();
+            seen.insert(rng.next_u64());
+            thread::sleep(Duration::from_micros(1));
+        }
+        assert!(seen.len() > 1, "Entropy seeds produced identical outputs");
+    }
+
+    #[test]
+    fn test_prng_uniform_distribution() {
+        let mut rng = Prng::new(12345);
+        let mut counts = [0usize; 10];
+        for _ in 0..10000 {
+            let v = rng.random_range(0..10usize);
+            counts[v] += 1;
+        }
+        for &c in &counts {
+            // "PRNG counts far from uniform: {c}"
+            assert!((c as isize - 1000).abs() < 150);
+        }
+    }
+
+    #[test]
+    fn test_prng_different_seeds_different_output() {
+        let mut a = Prng::new(1);
+        let mut b = Prng::new(2);
+        let va = a.next_u64();
+        let vb = b.next_u64();
+        assert_ne!(va, vb);
+    }
+
+    #[test]
+    fn test_prng_gen_bool_varies() {
+        let mut rng = Prng::new(99);
+        let mut seen_true = false;
+        let mut seen_false = false;
+        for _ in 0..100 {
+            if rng.gen_bool() {
+                seen_true = true;
+            } else {
+                seen_false = true;
+            }
+        }
+        assert!(seen_true && seen_false);
+    }
+
+    #[test]
+    fn test_random_range_single_usize() {
+        let mut rng = Prng::new(42);
+        for _ in 0..10 {
+            let v = rng.random_range(5..6);
+            assert_eq!(v, 5);
+        }
+    }
+
+    #[test]
+    fn test_random_range_single_f64() {
+        let mut rng = Prng::new(42);
+        for _ in 0..10 {
+            let v = rng.random_range(1.234..1.235);
+            assert!(v >= 1.234 && v < 1.235);
+        }
+    }
+
+    #[test]
+    fn test_prng_normal_zero_sd() {
+        let mut rng = Prng::new(7);
+        for _ in 0..5 {
+            let v = rng.normal(3.0, 0.0);
+            assert_eq!(v, 3.0);
+        }
+    }
+
+    #[test]
+    fn test_random_range_extreme_usize() {
+        let mut rng = Prng::new(5);
+        for _ in 0..10 {
+            let v = rng.random_range(0..usize::MAX);
+            assert!(v < usize::MAX);
+        }
+    }
+
+    #[test]
+    fn test_prng_chi_square_uniform() {
+        let mut rng = Prng::new(12345);
+        let mut counts = [0usize; 10];
+        let samples = 10000;
+        for _ in 0..samples {
+            let v = rng.random_range(0..10usize);
+            counts[v] += 1;
+        }
+        let expected = samples as f64 / 10.0;
+        let chi2: f64 = counts
+            .iter()
+            .map(|&c| {
+                let diff = c as f64 - expected;
+                diff * diff / expected
+            })
+            .sum();
+        //  "chi-square statistic too high: {chi2}"
+        assert!(chi2 < 20.0);
+    }
+
+    #[test]
+    fn test_prng_monobit() {
+        let mut rng = Prng::new(42);
+        let mut ones = 0usize;
+        let samples = 1000;
+        for _ in 0..samples {
+            ones += rng.next_u64().count_ones() as usize;
+        }
+        let total_bits = samples * 64;
+        let ratio = ones as f64 / total_bits as f64;
+        // "bit ratio far from 0.5: {ratio}"
+        assert!((ratio - 0.5).abs() < 0.01);
+    }
+}
--- a/src/random/random_core.rs
+++ b/src/random/random_core.rs
@@ -0,0 +1,106 @@
+//! Core traits for random number generators and sampling ranges.
+//!
+//! ```
+//! use rustframe::random::{rng, Rng};
+//! let mut r = rng();
+//! let value: f64 = r.random_range(0.0..1.0);
+//! assert!(value >= 0.0 && value < 1.0);
+//! ```
+use std::f64::consts::PI;
+use std::ops::Range;
+
+/// Trait implemented by random number generators.
+pub trait Rng {
+    /// Generate the next random `u64` value.
+    fn next_u64(&mut self) -> u64;
+
+    /// Generate a value uniformly in the given range.
+    fn random_range<T>(&mut self, range: Range<T>) -> T
+    where
+        T: RangeSample,
+    {
+        T::from_u64(self.next_u64(), &range)
+    }
+
+    /// Generate a boolean with probability 0.5 of being `true`.
+    fn gen_bool(&mut self) -> bool {
+        self.random_range(0..2usize) == 1
+    }
+
+    /// Sample from a normal distribution using the Box-Muller transform.
+    fn normal(&mut self, mean: f64, sd: f64) -> f64 {
+        let u1 = self.random_range(0.0..1.0);
+        let u2 = self.random_range(0.0..1.0);
+        mean + sd * (-2.0 * u1.ln()).sqrt() * (2.0 * PI * u2).cos()
+    }
+}
+
+/// Conversion from a raw `u64` into a type within a range.
+pub trait RangeSample: Sized {
+    fn from_u64(value: u64, range: &Range<Self>) -> Self;
+}
+
+impl RangeSample for usize {
+    fn from_u64(value: u64, range: &Range<Self>) -> Self {
+        let span = range.end - range.start;
+        (value as usize % span) + range.start
+    }
+}
+
+impl RangeSample for f64 {
+    fn from_u64(value: u64, range: &Range<Self>) -> Self {
+        let span = range.end - range.start;
+        range.start + (value as f64 / u64::MAX as f64) * span
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_range_sample_usize_boundary() {
+        assert_eq!(<usize as RangeSample>::from_u64(0, &(0..1)), 0);
+        assert_eq!(<usize as RangeSample>::from_u64(u64::MAX, &(0..1)), 0);
+    }
+
+    #[test]
+    fn test_range_sample_f64_boundary() {
+        let v0 = <f64 as RangeSample>::from_u64(0, &(0.0..1.0));
+        let vmax = <f64 as RangeSample>::from_u64(u64::MAX, &(0.0..1.0));
+        assert!(v0 >= 0.0 && v0 < 1.0);
+        assert!(vmax > 0.999999999999 && vmax <= 1.0);
+    }
+
+    #[test]
+    fn test_range_sample_usize_varied() {
+        for i in 0..5 {
+            let v = <usize as RangeSample>::from_u64(i, &(10..15));
+            assert!(v >= 10 && v < 15);
+        }
+    }
+
+    #[test]
+    fn test_range_sample_f64_span() {
+        for val in [0, u64::MAX / 2, u64::MAX] {
+            let f = <f64 as RangeSample>::from_u64(val, &(2.0..4.0));
+            assert!(f >= 2.0 && f <= 4.0);
+        }
+    }
+
+    #[test]
+    fn test_range_sample_usize_single_value() {
+        for val in [0, 1, u64::MAX] {
+            let n = <usize as RangeSample>::from_u64(val, &(5..6));
+            assert_eq!(n, 5);
+        }
+    }
+
+    #[test]
+    fn test_range_sample_f64_negative_range() {
+        for val in [0, u64::MAX / 3, u64::MAX] {
+            let f = <f64 as RangeSample>::from_u64(val, &(-2.0..2.0));
+            assert!(f >= -2.0 && f <= 2.0);
+        }
+    }
+}
--- a/src/random/seq.rs
+++ b/src/random/seq.rs
@@ -0,0 +1,113 @@
+//! Extensions for shuffling slices with a random number generator.
+//!
+//! ```
+//! use rustframe::random::{rng, SliceRandom};
+//! let mut data = [1, 2, 3];
+//! data.shuffle(&mut rng());
+//! assert_eq!(data.len(), 3);
+//! ```
+use crate::random::Rng;
+
+/// Trait for randomizing slices.
+pub trait SliceRandom {
+    /// Shuffle the slice in place using the provided RNG.
+    fn shuffle<R: Rng>(&mut self, rng: &mut R);
+}
+
+impl<T> SliceRandom for [T] {
+    fn shuffle<R: Rng>(&mut self, rng: &mut R) {
+        for i in (1..self.len()).rev() {
+            let j = rng.random_range(0..(i + 1));
+            self.swap(i, j);
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::random::{CryptoRng, Prng};
+
+    #[test]
+    fn test_shuffle_slice() {
+        let mut rng = Prng::new(3);
+        let mut arr = [1, 2, 3, 4, 5];
+        let orig = arr.clone();
+        arr.shuffle(&mut rng);
+        assert_eq!(arr.len(), orig.len());
+        let mut sorted = arr.to_vec();
+        sorted.sort();
+        assert_eq!(sorted, orig.to_vec());
+    }
+
+    #[test]
+    fn test_slice_shuffle_deterministic_with_prng() {
+        let mut rng1 = Prng::new(11);
+        let mut rng2 = Prng::new(11);
+        let mut a = [1u8, 2, 3, 4, 5, 6, 7, 8, 9];
+        let mut b = a.clone();
+        a.shuffle(&mut rng1);
+        b.shuffle(&mut rng2);
+        assert_eq!(a, b);
+    }
+
+    #[test]
+    fn test_slice_shuffle_crypto_random_changes() {
+        let mut rng1 = CryptoRng::new();
+        let mut rng2 = CryptoRng::new();
+        let orig = [1u8, 2, 3, 4, 5, 6, 7, 8, 9];
+        let mut a = orig.clone();
+        let mut b = orig.clone();
+        a.shuffle(&mut rng1);
+        b.shuffle(&mut rng2);
+        assert!(a != orig || b != orig, "Shuffles did not change order");
+        assert_ne!(a, b, "Two Crypto RNG shuffles produced same order");
+    }
+
+    #[test]
+    fn test_shuffle_single_element_no_change() {
+        let mut rng = Prng::new(1);
+        let mut arr = [42];
+        arr.shuffle(&mut rng);
+        assert_eq!(arr, [42]);
+    }
+
+    #[test]
+    fn test_multiple_shuffles_different_results() {
+        let mut rng = Prng::new(5);
+        let mut arr1 = [1, 2, 3, 4];
+        let mut arr2 = [1, 2, 3, 4];
+        arr1.shuffle(&mut rng);
+        arr2.shuffle(&mut rng);
+        assert_ne!(arr1, arr2);
+    }
+
+    #[test]
+    fn test_shuffle_empty_slice() {
+        let mut rng = Prng::new(1);
+        let mut arr: [i32; 0] = [];
+        arr.shuffle(&mut rng);
+        assert!(arr.is_empty());
+    }
+
+    #[test]
+    fn test_shuffle_three_uniform() {
+        use std::collections::HashMap;
+        let mut rng = Prng::new(123);
+        let mut counts: HashMap<[u8; 3], usize> = HashMap::new();
+        for _ in 0..6000 {
+            let mut arr = [1u8, 2, 3];
+            arr.shuffle(&mut rng);
+            *counts.entry(arr).or_insert(0) += 1;
+        }
+        let expected = 1000.0;
+        let chi2: f64 = counts
+            .values()
+            .map(|&c| {
+                let diff = c as f64 - expected;
+                diff * diff / expected
+            })
+            .sum();
+        assert!(chi2 < 30.0, "shuffle chi-square too high: {chi2}");
+    }
+}
--- a/src/utils/dateutils/dates.rs
+++ b/src/utils/dateutils/dates.rs
@@ -1,3 +1,10 @@
+//! Generation and manipulation of calendar date sequences.
+//!
+//! ```
+//! use rustframe::utils::dateutils::dates::{DateFreq, DatesList};
+//! let list = DatesList::new("2024-01-01".into(), "2024-01-03".into(), DateFreq::Daily);
+//! assert_eq!(list.count().unwrap(), 3);
+//! ```
 use chrono::{Datelike, Duration, NaiveDate, Weekday};
 use std::collections::HashMap;
 use std::error::Error;
--- a/src/utils/dateutils/mod.rs
+++ b/src/utils/dateutils/mod.rs
@@ -1,3 +1,13 @@
+//! Generators for sequences of calendar and business dates.
+//!
+//! See [`dates`] for all-day calendars and [`bdates`] for business-day aware
+//! variants.
+//!
+//! ```
+//! use rustframe::utils::dateutils::{DatesList, DateFreq};
+//! let list = DatesList::new("2024-01-01".into(), "2024-01-02".into(), DateFreq::Daily);
+//! assert_eq!(list.count().unwrap(), 2);
+//! ```
 pub mod bdates;
 pub mod dates;

--- a/src/utils/mod.rs
+++ b/src/utils/mod.rs
@@ -1,3 +1,14 @@
+//! Assorted helper utilities.
+//!
+//! Currently this module exposes date generation utilities in [`dateutils`](crate::utils::dateutils),
+//! including calendar and business date sequences.
+//!
+//! ```
+//! use rustframe::utils::DatesList;
+//! use rustframe::utils::DateFreq;
+//! let dates = DatesList::new("2024-01-01".into(), "2024-01-03".into(), DateFreq::Daily);
+//! assert_eq!(dates.count().unwrap(), 3);
+//! ```
 pub mod dateutils;

 pub use dateutils::{BDateFreq, BDatesGenerator, BDatesList};