Correlation & Fusion

# ============================================================================
#  Sensor fusion & correlation — the whole idea in plain Python (stdlib only)
# ============================================================================
#  Two jobs:
#    1. fuse()       blend several readings of ONE number into one best guess.
#    2. similarity() score how likely two readings came from the SAME object.
#  Nothing to install. Run it with:  python3 fusion_demo.py
# ============================================================================

import math


# ---------------------------------------------------------------------------
#  1. INVERSE-VARIANCE FUSION  —  "trust the steadier sensor more"
#
#  Every reading comes with a `variance`: how jittery that sensor is.
#  Small variance = precise, big variance = sloppy. We give each reading a
#  weight of 1/variance, then take a weighted average — so a sensor with 4x
#  less variance ends up counting 4x as much.
# ---------------------------------------------------------------------------
def fuse(values, variances):
    # Step 1: turn each variance into a "trust" score (1 / variance).
    inverse_variances = [1.0 / v for v in variances]

    # Step 2: add up all the trust, so we can turn the scores into fractions.
    total_trust = math.fsum(inverse_variances)

    # Step 3: each weight is that sensor's share of the total trust.
    #         All the weights add up to 1.0.
    weights = [iv / total_trust for iv in inverse_variances]

    # Step 4: the fused value is the weighted average of the readings.
    mean = math.fsum(w * x for w, x in zip(weights, values))

    # Step 5: report how spread out the readings were around that average.
    #         This grows when the sensors DISAGREE (honest error bars).
    n = len(values)
    spread = math.fsum(w * (x - mean) ** 2 for w, x in zip(weights, values))
    variance = spread * n / (n - 1)          # small-sample (Bessel) correction

    return {
        "value": mean,
        "variance": variance,                    # empirical spread
        "std_dev": math.sqrt(variance),
        "std_error": math.sqrt(variance) / math.sqrt(n),
        "combined_variance": 1.0 / total_trust,  # classical, always tighter
        "weights": weights,
        "n": n,
    }


# ---------------------------------------------------------------------------
#  2. CORRELATION  —  "are these two readings the same object?"
#
#  Measure how many error-bars apart the two readings are (call it `t`), then
#  ask: if they really were the same thing, how often would pure chance push
#  them this far apart? That probability is the "similarity" (0..1).
#  High similarity = probably the same object.
#
#  Working out that probability needs the area under a bell-shaped curve. The
#  three helpers below (lgamma is built in; _betacf and _betai are written out)
#  are standard, well-tested numerical recipes for exactly that. Treat them as a
#  sealed calculator. The plain-English part is similarity() and combine().
# ---------------------------------------------------------------------------

# _betacf + _betai together compute the "regularised incomplete beta", which is
# the math behind the Student-t tail. Method: Lentz's continued fraction
# (Numerical Recipes section 6.4). You should rarely need to touch this.
def _betacf(a, b, x):
    tiny = 1e-30          # a stand-in for "basically zero", avoids divide-by-0
    c = 1.0
    d = 1.0 - (a + b) * x / (a + 1.0)
    d = 1.0 / (tiny if abs(d) < tiny else d)
    h = d
    for m in range(1, 200):           # loop until the fraction stops changing
        m2 = 2 * m
        # --- even step of the continued fraction ---
        aa = m * (b - m) * x / ((a - 1.0 + m2) * (a + m2))
        d = 1.0 / (1.0 + aa * d)
        c = 1.0 + aa / c
        h *= d * c
        # --- odd step of the continued fraction ---
        aa = -(a + m) * (a + b + m) * x / ((a + m2) * (a + 1.0 + m2))
        d = 1.0 / (1.0 + aa * d)
        c = 1.0 + aa / c
        delta = d * c
        h *= delta
        if abs(delta - 1.0) < 3e-11:  # converged: more steps wouldn't change it
            break
    return h


def _betai(a, b, x):
    if x <= 0.0:
        return 0.0
    if x >= 1.0:
        return 1.0
    # the log-scale prefactor keeps the numbers from overflowing
    front = math.exp(math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b)
                     + a * math.log(x) + b * math.log(1.0 - x))
    # use whichever side of the fraction converges faster
    if x < (a + 1.0) / (a + b + 2.0):
        return front * _betacf(a, b, x) / a
    return 1.0 - front * _betacf(b, a, 1.0 - x) / b


def similarity(mean_a, se_a, n_a, mean_b, se_b, n_b):
    # how far apart are the two readings, measured in combined error-bars?
    combined_error = math.hypot(se_a, se_b)        # = sqrt(se_a**2 + se_b**2)
    t = (mean_a - mean_b) / combined_error

    # with 20+ samples the curve is "normal"; its tail uses erfc.
    if min(n_a, n_b) >= 20:                         # z-test
        return math.erfc(abs(t) / math.sqrt(2.0))

    # with few samples, use the fatter-tailed Student-t curve instead.
    dof = n_a + n_b - 2                             # t-test "degrees of freedom"
    return _betai(dof / 2.0, 0.5, dof / (dof + t * t))


def combine(scores):
    # Fold several feature scores (length, width, ...) into one number by taking
    # their geometric mean. Each score is floored at 1e-6 so a single 0 can't
    # wipe out the whole product.
    logs = [math.log(max(1e-6, s)) for s in scores]
    return math.exp(math.fsum(logs) / len(scores))


# ---------------------------------------------------------------------------
#  Run it on the two real sensors (EO first, radar second) and print the
#  verified numbers.
# ---------------------------------------------------------------------------
print(fuse([10.0, 10.00008], [25.0, 100.0]))     # latitude
print(fuse([123.0, 121.0],   [100.0, 400.0]))    # altitude
sl = similarity(120, 0.87, 12, 121, 0.73, 30)    # length feature
sw = similarity(18,  0.29, 12, 17.5, 0.37, 30)   # width feature
print(sl, sw, combine([sl, sw]))

# latitude → {'value': 10.000016, 'variance': 2.048e-09, ..., 'combined_variance': 20.0, 'weights': [0.8, 0.2], 'n': 2}
# altitude → {'value': 122.6, 'variance': 1.28, 'std_dev': 1.1313708..., 'combined_variance': 80.0, 'weights': [0.8, 0.2], 'n': 2}
# 0.383838  0.293894  0.335868   →  combined ≥ 0.05  →  same object

// ============================================================================
//  Sensor fusion & correlation — the whole idea in plain TypeScript
// ============================================================================
//  Two jobs:
//    1. fuse()       blend several readings of ONE number into one best guess.
//    2. similarity() score how likely two readings came from the SAME object.
//  No dependencies. Run with:  npx tsx fusionDemo.ts   (or compile with tsc)
// ============================================================================

// ---------------------------------------------------------------------------
//  1. INVERSE-VARIANCE FUSION  —  "trust the steadier sensor more"
//     Each reading has a `variance` (small = precise, big = sloppy). We weight
//     every reading by 1/variance, then take the weighted average — so a sensor
//     with 4x less variance counts 4x as much.
// ---------------------------------------------------------------------------
function fuse(values: number[], variances: number[]) {
  // Step 1: turn each variance into a "trust" score (1 / variance).
  const inverseVariances = variances.map((v) => 1 / v);

  // Step 2: total trust, used to turn the scores into fractions.
  const totalTrust = inverseVariances.reduce((a, b) => a + b, 0);

  // Step 3: each weight is that sensor's share of the total (weights sum to 1).
  const weights = inverseVariances.map((iv) => iv / totalTrust);

  // Step 4: the fused value is the weighted average of the readings.
  const mean = weights.reduce((sum, w, i) => sum + w * values[i], 0);

  // Step 5: how spread out the readings were (grows when sensors disagree).
  const n = values.length;
  const spread = weights.reduce((sum, w, i) => sum + w * (values[i] - mean) ** 2, 0);
  const variance = (spread * n) / (n - 1);   // small-sample (Bessel) correction

  return {
    value: mean,
    variance,                                 // empirical spread
    stdDev: Math.sqrt(variance),
    stdError: Math.sqrt(variance) / Math.sqrt(n),
    combinedVariance: 1 / totalTrust,         // classical, always tighter
    weights,
    n,
  };
}

// ---------------------------------------------------------------------------
//  2. CORRELATION  —  "are these two readings the same object?"
//     Measure how many error-bars apart they are, then ask how often pure
//     chance would spread them that far. That probability is the "similarity".
//
//     lgamma / betacf / betai below are standard numerical recipes that compute
//     the area under the Student-t curve. Treat them as a sealed calculator;
//     the plain-English part is similarity() and combine() at the bottom.
// ---------------------------------------------------------------------------

// lgamma(z) = log of the Gamma function, via the Lanczos approximation.
// (JavaScript has no built-in lgamma, so we supply one.)
const LANCZOS = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
  771.32342877765313, -176.61502916214059, 12.507343278686905,
  -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7];
function lgamma(z: number): number {
  // reflection formula keeps the approximation valid for small z
  if (z < 0.5) return Math.log(Math.PI / Math.sin(Math.PI * z)) - lgamma(1 - z);
  z -= 1;
  let x = LANCZOS[0];
  for (let i = 1; i < 9; i++) x += LANCZOS[i] / (z + i);
  const t = z + 7.5;
  return 0.5 * Math.log(2 * Math.PI) + (z + 0.5) * Math.log(t) - t + Math.log(x);
}

// Continued fraction for the incomplete beta (Lentz's method, Num. Recipes 6.4).
function betacf(a: number, b: number, x: number): number {
  const tiny = 1e-30;         // stand-in for "basically zero", avoids /0
  let c = 1;
  let d = 1 - ((a + b) * x) / (a + 1);
  d = 1 / (Math.abs(d) < tiny ? tiny : d);
  let h = d;
  for (let m = 1; m < 200; m++) {     // loop until it stops changing
    const m2 = 2 * m;
    // even step
    let aa = (m * (b - m) * x) / ((a - 1 + m2) * (a + m2));
    d = 1 / (1 + aa * d);
    c = 1 + aa / c;
    h *= d * c;
    // odd step
    aa = (-(a + m) * (a + b + m) * x) / ((a + m2) * (a + 1 + m2));
    d = 1 / (1 + aa * d);
    c = 1 + aa / c;
    const delta = d * c;
    h *= delta;
    if (Math.abs(delta - 1) < 3e-11) break;   // converged
  }
  return h;
}

// Regularised incomplete beta I_x(a, b) — the Student-t tail builds on this.
function betai(a: number, b: number, x: number): number {
  if (x <= 0) return 0;
  if (x >= 1) return 1;
  const front = Math.exp(
    lgamma(a + b) - lgamma(a) - lgamma(b) + a * Math.log(x) + b * Math.log(1 - x));
  // use whichever side of the fraction converges faster
  return x < (a + 1) / (a + b + 2)
    ? (front * betacf(a, b, x)) / a
    : 1 - (front * betacf(b, a, 1 - x)) / b;
}

// erfc(x) = the normal-curve tail, via the Abramowitz & Stegun 7.1.26 polynomial.
function erfc(x: number): number {
  const z = Math.abs(x);
  const t = 1 / (1 + 0.3275911 * z);
  const y = t * (0.254829592 + t * (-0.284496736 + t * (1.421413741
          + t * (-1.453152027 + t * 1.061405429)))) * Math.exp(-z * z);
  return x >= 0 ? y : 2 - y;
}

function similarity(mA: number, seA: number, nA: number,
                    mB: number, seB: number, nB: number): number {
  // how far apart, measured in combined error-bars
  const combinedError = Math.hypot(seA, seB);    // = sqrt(seA^2 + seB^2)
  const t = (mA - mB) / combinedError;
  if (Math.min(nA, nB) >= 20) {                  // 20+ samples → normal curve
    return erfc(Math.abs(t) / Math.SQRT2);       // z-test
  }
  const dof = nA + nB - 2;                        // few samples → Student-t
  return betai(dof / 2, 0.5, dof / (dof + t * t)); // t-test
}

function combine(scores: number[]): number {
  // geometric mean; floor each score at 1e-6 so one 0 can't wipe out the rest
  const sumOfLogs = scores.reduce((s, v) => s + Math.log(Math.max(1e-6, v)), 0);
  return Math.exp(sumOfLogs / scores.length);
}

// Run it on the two real sensors (EO first, radar second).
console.log(fuse([10.0, 10.00008], [25, 100]));   // latitude
console.log(fuse([123, 121], [100, 400]));        // altitude
const sl = similarity(120, 0.87, 12, 121, 0.73, 30);   // length
const sw = similarity(18, 0.29, 12, 17.5, 0.37, 30);   // width
console.log(sl, sw, combine([sl, sw]));

// latitude → { value: 10.000016, variance: 2.048e-9, ..., combinedVariance: 20, weights: [0.8, 0.2], n: 2 }
// altitude → { value: 122.6, variance: 1.28, stdDev: 1.1313708..., combinedVariance: 80, weights: [0.8, 0.2], n: 2 }
// 0.383838  0.293894  0.335868   →  combined ≥ 0.05  →  same object

// ============================================================================
//  Sensor fusion & correlation — the whole idea in plain Rust (no crates)
// ============================================================================
//  Two jobs:
//    1. fuse()       blend several readings of ONE number into one best guess.
//    2. similarity() score how likely two readings came from the SAME object.
//  Run with:  rustc fusion_demo.rs && ./fusion_demo
// ============================================================================

// ---------------------------------------------------------------------------
//  1. INVERSE-VARIANCE FUSION  —  "trust the steadier sensor more"
//     Each reading has a `variance` (small = precise, big = sloppy). We weight
//     every reading by 1/variance and take the weighted average — so a sensor
//     with 4x less variance counts 4x as much.
// ---------------------------------------------------------------------------
struct Fused {
    value: f64,
    variance: f64,             // empirical spread (grows when sensors disagree)
    std_dev: f64,
    std_error: f64,
    combined_variance: f64,    // classical fused variance, always tighter
    weights: Vec<f64>,
    n: usize,
}

fn fuse(values: &[f64], variances: &[f64]) -> Fused {
    // Step 1: turn each variance into a "trust" score (1 / variance).
    let inverse_variances: Vec<f64> = variances.iter().map(|v| 1.0 / v).collect();

    // Step 2: total trust, used to turn the scores into fractions.
    let total_trust: f64 = inverse_variances.iter().sum();

    // Step 3: each weight is that sensor's share of the total (weights sum to 1).
    let weights: Vec<f64> = inverse_variances.iter().map(|iv| iv / total_trust).collect();

    // Step 4: the fused value is the weighted average of the readings.
    let mean: f64 = weights.iter().zip(values).map(|(w, x)| w * x).sum();

    // Step 5: how spread out the readings were around that average.
    let n = values.len();
    let spread: f64 = weights.iter().zip(values).map(|(w, x)| w * (x - mean).powi(2)).sum();
    let variance = spread * n as f64 / (n as f64 - 1.0);   // small-sample correction

    Fused {
        value: mean,
        variance,
        std_dev: variance.sqrt(),
        std_error: variance.sqrt() / (n as f64).sqrt(),
        combined_variance: 1.0 / total_trust,
        weights,
        n,
    }
}

// ---------------------------------------------------------------------------
//  2. CORRELATION  —  "are these two readings the same object?"
//     Measure how many error-bars apart they are, then ask how often pure
//     chance would spread them that far. That probability is the "similarity".
//
//     lgamma / betacf / betai are standard numerical recipes computing the area
//     under the Student-t curve — treat them as a sealed calculator. The plain-
//     English part is similarity() and combine().
// ---------------------------------------------------------------------------

// lgamma(z) = log of the Gamma function, via the Lanczos approximation.
fn lgamma(z: f64) -> f64 {
    const C: [f64; 9] = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
        771.32342877765313, -176.61502916214059, 12.507343278686905,
        -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7];
    use std::f64::consts::PI;
    if z < 0.5 {                                  // reflection formula for small z
        return (PI / (PI * z).sin()).ln() - lgamma(1.0 - z);
    }
    let z = z - 1.0;
    let mut x = C[0];
    for i in 1..9 {
        x += C[i] / (z + i as f64);
    }
    let t = z + 7.5;
    0.5 * (2.0 * PI).ln() + (z + 0.5) * t.ln() - t + x.ln()
}

// Continued fraction for the incomplete beta (Lentz's method, Num. Recipes 6.4).
fn betacf(a: f64, b: f64, x: f64) -> f64 {
    let tiny = 1e-30;                 // stand-in for "basically zero", avoids /0
    let mut c = 1.0;
    let mut d = 1.0 - (a + b) * x / (a + 1.0);
    if d.abs() < tiny {
        d = tiny;
    }
    d = 1.0 / d;
    let mut h = d;
    for m in 1..200 {                 // loop until it stops changing
        let (mf, m2) = (m as f64, (2 * m) as f64);
        // even step
        let mut aa = mf * (b - mf) * x / ((a - 1.0 + m2) * (a + m2));
        d = 1.0 / (1.0 + aa * d);
        c = 1.0 + aa / c;
        h *= d * c;
        // odd step
        aa = -(a + mf) * (a + b + mf) * x / ((a + m2) * (a + 1.0 + m2));
        d = 1.0 / (1.0 + aa * d);
        c = 1.0 + aa / c;
        let delta = d * c;
        h *= delta;
        if (delta - 1.0).abs() < 3e-11 {          // converged
            break;
        }
    }
    h
}

// Regularised incomplete beta I_x(a, b) — the Student-t tail builds on this.
fn betai(a: f64, b: f64, x: f64) -> f64 {
    if x <= 0.0 {
        return 0.0;
    }
    if x >= 1.0 {
        return 1.0;
    }
    let front = (lgamma(a + b) - lgamma(a) - lgamma(b)
        + a * x.ln() + b * (1.0 - x).ln()).exp();
    // use whichever side of the fraction converges faster
    if x < (a + 1.0) / (a + b + 2.0) {
        front * betacf(a, b, x) / a
    } else {
        1.0 - front * betacf(b, a, 1.0 - x) / b
    }
}

// erfc(x) = the normal-curve tail, via the Abramowitz & Stegun 7.1.26 polynomial.
fn erfc(x: f64) -> f64 {
    let z = x.abs();
    let t = 1.0 / (1.0 + 0.3275911 * z);
    let y = t * (0.254829592 + t * (-0.284496736 + t * (1.421413741
          + t * (-1.453152027 + t * 1.061405429)))) * (-z * z).exp();
    if x >= 0.0 { y } else { 2.0 - y }
}

fn similarity(ma: f64, sea: f64, na: u32, mb: f64, seb: f64, nb: u32) -> f64 {
    let combined_error = sea.hypot(seb);          // = sqrt(sea^2 + seb^2)
    let t = (ma - mb) / combined_error;
    if na.min(nb) >= 20 {                          // 20+ samples → normal curve
        erfc(t.abs() / std::f64::consts::SQRT_2)   // z-test
    } else {
        let dof = (na + nb - 2) as f64;            // few samples → Student-t
        betai(dof / 2.0, 0.5, dof / (dof + t * t)) // t-test
    }
}

fn combine(scores: &[f64]) -> f64 {
    // geometric mean; floor each score at 1e-6 so one 0 can't wipe out the rest
    let sum_of_logs: f64 = scores.iter().map(|s| s.max(1e-6).ln()).sum();
    (sum_of_logs / scores.len() as f64).exp()
}

fn main() {
    let lat = fuse(&[10.0, 10.00008], &[25.0, 100.0]);   // latitude
    let alt = fuse(&[123.0, 121.0], &[100.0, 400.0]);    // altitude
    println!("lat  value={} cvar={} weights={:?}", lat.value, lat.combined_variance, lat.weights);
    println!("alt  value={} variance={} std_dev={} cvar={}", alt.value, alt.variance, alt.std_dev, alt.combined_variance);
    let sl = similarity(120.0, 0.87, 12, 121.0, 0.73, 30);   // length
    let sw = similarity(18.0, 0.29, 12, 17.5, 0.37, 30);     // width
    println!("length={:.6} width={:.6} combined={:.6}", sl, sw, combine(&[sl, sw]));
    // lat  value=10.000016 cvar=20 weights=[0.8, 0.2]
    // alt  value=122.6 variance=1.28 std_dev=1.1313708498984762 cvar=80
    // length=0.383838 width=0.293894 combined=0.335868   →  same object
}

import java.util.Arrays;

// ============================================================================
//  Sensor fusion & correlation — the whole idea in plain Java (no libraries)
// ============================================================================
//  Two jobs:
//    1. fuse()       blend several readings of ONE number into one best guess.
//    2. similarity() score how likely two readings came from the SAME object.
//  Run with:  javac FusionDemo.java && java FusionDemo
// ============================================================================
public class FusionDemo {

    // -----------------------------------------------------------------------
    //  1. INVERSE-VARIANCE FUSION  —  "trust the steadier sensor more"
    //     Each reading has a variance (small = precise, big = sloppy). We weight
    //     each reading by 1/variance and take the weighted average — so a sensor
    //     with 4x less variance counts 4x as much.
    // -----------------------------------------------------------------------
    record Fused(double value, double variance, double stdDev, double stdError,
                 double combinedVariance, double[] weights, int n) {}

    static Fused fuse(double[] values, double[] variances) {
        int n = values.length;

        // Step 1 + 2: trust = 1/variance per reading; add up the total trust.
        double totalTrust = 0;
        for (double v : variances) totalTrust += 1.0 / v;

        // Step 3 + 4: weight = each reading's share of the trust (sums to 1);
        //             the fused value is the weighted average.
        double[] weights = new double[n];
        double mean = 0;
        for (int i = 0; i < n; i++) {
            weights[i] = (1.0 / variances[i]) / totalTrust;
            mean += weights[i] * values[i];
        }

        // Step 5: how spread out the readings were (grows when sensors disagree).
        double spread = 0;
        for (int i = 0; i < n; i++) spread += weights[i] * Math.pow(values[i] - mean, 2);
        double variance = spread * (double) n / (n - 1);   // small-sample correction

        return new Fused(mean, variance, Math.sqrt(variance),
                Math.sqrt(variance) / Math.sqrt(n), 1.0 / totalTrust, weights, n);
    }

    // -----------------------------------------------------------------------
    //  2. CORRELATION  —  "are these two readings the same object?"
    //     Measure how many error-bars apart they are, then ask how often pure
    //     chance would spread them that far. That probability is the similarity.
    //
    //     lgamma / betacf / betai are standard numerical recipes that compute the
    //     area under the Student-t curve — treat them as a sealed calculator. The
    //     plain-English part is similarity() and combine().
    // -----------------------------------------------------------------------

    // lgamma(z) = log of the Gamma function, via the Lanczos approximation.
    static final double[] LANCZOS = {0.99999999999980993, 676.5203681218851, -1259.1392167224028,
        771.32342877765313, -176.61502916214059, 12.507343278686905,
        -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
    static double lgamma(double z) {
        if (z < 0.5)                                  // reflection formula for small z
            return Math.log(Math.PI / Math.sin(Math.PI * z)) - lgamma(1 - z);
        z -= 1;
        double x = LANCZOS[0];
        for (int i = 1; i < 9; i++) x += LANCZOS[i] / (z + i);
        double t = z + 7.5;
        return 0.5 * Math.log(2 * Math.PI) + (z + 0.5) * Math.log(t) - t + Math.log(x);
    }

    // Continued fraction for the incomplete beta (Lentz's method, Num. Recipes 6.4).
    static double betacf(double a, double b, double x) {
        double tiny = 1e-30;          // stand-in for "basically zero", avoids /0
        double c = 1;
        double d = 1 - (a + b) * x / (a + 1);
        if (Math.abs(d) < tiny) d = tiny;
        d = 1 / d;
        double h = d;
        for (int m = 1; m < 200; m++) {     // loop until it stops changing
            int m2 = 2 * m;
            // even step
            double aa = m * (b - m) * x / ((a - 1 + m2) * (a + m2));
            d = 1 / (1 + aa * d);
            c = 1 + aa / c;
            h *= d * c;
            // odd step
            aa = -(a + m) * (a + b + m) * x / ((a + m2) * (a + 1 + m2));
            d = 1 / (1 + aa * d);
            c = 1 + aa / c;
            double delta = d * c;
            h *= delta;
            if (Math.abs(delta - 1) < 3e-11) break;   // converged
        }
        return h;
    }

    // Regularised incomplete beta I_x(a, b) — the Student-t tail builds on this.
    static double betai(double a, double b, double x) {
        if (x <= 0) return 0;
        if (x >= 1) return 1;
        double front = Math.exp(lgamma(a + b) - lgamma(a) - lgamma(b)
                + a * Math.log(x) + b * Math.log(1 - x));
        // use whichever side of the fraction converges faster
        return x < (a + 1) / (a + b + 2) ? front * betacf(a, b, x) / a
                                         : 1 - front * betacf(b, a, 1 - x) / b;
    }

    // erfc(x) = the normal-curve tail, via the Abramowitz & Stegun 7.1.26 polynomial.
    static double erfc(double x) {
        double z = Math.abs(x);
        double t = 1 / (1 + 0.3275911 * z);
        double y = t * (0.254829592 + t * (-0.284496736 + t * (1.421413741
                 + t * (-1.453152027 + t * 1.061405429)))) * Math.exp(-z * z);
        return x >= 0 ? y : 2 - y;
    }

    static double similarity(double mA, double seA, int nA, double mB, double seB, int nB) {
        double combinedError = Math.hypot(seA, seB);   // = sqrt(seA^2 + seB^2)
        double t = (mA - mB) / combinedError;
        if (Math.min(nA, nB) >= 20)                     // 20+ samples → normal curve
            return erfc(Math.abs(t) / Math.sqrt(2));    // z-test
        double dof = nA + nB - 2;                       // few samples → Student-t
        return betai(dof / 2, 0.5, dof / (dof + t * t)); // t-test
    }

    static double combine(double... scores) {
        // geometric mean; floor each score at 1e-6 so one 0 can't wipe out the rest
        double sumOfLogs = 0;
        for (double v : scores) sumOfLogs += Math.log(Math.max(1e-6, v));
        return Math.exp(sumOfLogs / scores.length);
    }

    public static void main(String[] args) {
        Fused lat = fuse(new double[]{10.0, 10.00008}, new double[]{25, 100});   // latitude
        Fused alt = fuse(new double[]{123, 121}, new double[]{100, 400});        // altitude
        System.out.printf("lat  value=%.6f cvar=%.1f weights=%s%n",
                lat.value(), lat.combinedVariance(), Arrays.toString(lat.weights()));
        System.out.printf("alt  value=%.1f variance=%.1f std_dev=%.6f cvar=%.1f%n",
                alt.value(), alt.variance(), alt.stdDev(), alt.combinedVariance());
        double sl = similarity(120, 0.87, 12, 121, 0.73, 30);   // length
        double sw = similarity(18, 0.29, 12, 17.5, 0.37, 30);   // width
        System.out.printf("length=%.6f width=%.6f combined=%.6f%n", sl, sw, combine(sl, sw));
        // lat  value=10.000016 cvar=20.0 weights=[0.8, 0.2]
        // alt  value=122.6 variance=1.28 std_dev=1.131371 cvar=80.0
        // length=0.383838 width=0.293894 combined=0.335868   →  same object
    }
}