两个预测峰都合理，平均值为什么落在低密度区？

"""Fixed Gaussian mixtures: density, squared error, likelihood and variance.

Python 3.10+; standard library only. No model training or real dataset.
"""
import json
import math
import platform
from decimal import Decimal, localcontext
from pathlib import Path


def validate(weights, means, sigmas):
    if not weights or not (len(weights) == len(means) == len(sigmas)):
        raise ValueError("nonempty component arrays must have equal lengths")
    if not all(math.isfinite(x) for row in (weights, means, sigmas) for x in row):
        raise ValueError("parameters must be finite")
    if any(w < 0 for w in weights) or not math.isclose(
        math.fsum(weights), 1.0, rel_tol=0.0, abs_tol=1e-12
    ):
        raise ValueError("weights must be nonnegative and sum to one")
    if any(s <= 0 for s in sigmas):
        raise ValueError("standard deviations must be positive")


def log_normal(x, mean, sigma):
    residual = (x - mean) / sigma
    return -0.5 * math.log(2.0 * math.pi) - math.log(sigma) - 0.5 * residual * residual


def logsumexp(values):
    top = max(values)
    return top + math.log(math.fsum(math.exp(v - top) for v in values))


def log_mixture(x, weights, means, sigmas):
    validate(weights, means, sigmas)
    if not math.isfinite(x):
        raise ValueError("observation must be finite")
    terms = [math.log(w) + log_normal(x, m, s)
             for w, m, s in zip(weights, means, sigmas) if w > 0]
    if not all(math.isfinite(t) for t in terms):
        raise ValueError("standardized residual exceeds this float implementation")
    return logsumexp(terms)


def moments(weights, means, sigmas):
    validate(weights, means, sigmas)
    mean = math.fsum(w * m for w, m in zip(weights, means))
    within = math.fsum(w * s * s for w, s in zip(weights, sigmas))
    between = math.fsum(w * (m - mean) ** 2 for w, m in zip(weights, means))
    return mean, within, between


def decimal_nll(x, weights, means, sigmas):
    """Independent reference: directly sum high-precision densities, then log.

    This intentionally does not call log_normal or logsumexp.
    """
    with localcontext() as ctx:
        ctx.prec = 70
        pi = Decimal("3.1415926535897932384626433832795028841971693993751058209749445923")
        xd = Decimal(str(x))
        total = Decimal(0)
        for w, m, s in zip(weights, means, sigmas):
            wd, md, sd = map(lambda v: Decimal(str(v)), (w, m, s))
            exponent = -((xd - md) ** 2) / (2 * sd * sd)
            total += wd * exponent.exp() / ((2 * pi).sqrt() * sd)
        return float(-total.ln())


def simpson(fn, lo, hi, intervals=16000):
    """Independent numerical integration over a specified finite interval."""
    if intervals <= 0 or intervals % 2:
        raise ValueError("Simpson integration needs a positive even count")
    step = (hi - lo) / intervals
    total = math.fsum((4 if i % 2 else 2) * fn(lo + i * step)
                      for i in range(1, intervals))
    return step * (fn(lo) + fn(hi) + total) / 3


def main():
    checks = []

    def check(name, condition):
        if not condition:
            raise AssertionError(name)
        checks.append(name)

    w, m, s = [0.5, 0.5], [-2.0, 2.0], [0.5, 0.5]
    mean, within, between = moments(w, m, s)
    variance = within + between
    density_mean = math.exp(log_mixture(mean, w, m, s))
    density_center = math.exp(log_mixture(2.0, w, m, s))
    mixture_nll = -log_mixture(2.0, w, m, s)
    component_nll = math.fsum(-wi * log_normal(2.0, mi, si)
                              for wi, mi, si in zip(w, m, s))
    check("density at the mean is lower than at a component center", density_mean < density_center)
    check("mixture NLL is below weighted component NLL here", mixture_nll < component_nll)
    check("within variance", math.isclose(within, 0.25))
    check("between variance", math.isclose(between, 4.0))

    # Integrate a separately written density formula, avoiding log_mixture.
    def density(x):
        return sum(wi * math.exp(-((x - mi) ** 2) / (2 * si * si))
                   / (si * math.sqrt(2 * math.pi)) for wi, mi, si in zip(w, m, s))

    mass = simpson(density, -8, 8)
    integrated_mean = simpson(lambda x: x * density(x), -8, 8)
    integrated_variance = simpson(lambda x: (x - mean) ** 2 * density(x), -8, 8)
    mse_at_mean = simpson(lambda x: (x - mean) ** 2 * density(x), -8, 8)
    mse_at_center = simpson(lambda x: (x - 2) ** 2 * density(x), -8, 8)
    check("quadrature mass", math.isclose(mass, 1.0, abs_tol=1e-10))
    check("quadrature mean", math.isclose(integrated_mean, mean, abs_tol=1e-10))
    check("quadrature total variance", math.isclose(integrated_variance, variance, abs_tol=1e-10))
    check("quadrature MSE at mean", math.isclose(mse_at_mean, 4.25, abs_tol=1e-10))
    check("quadrature MSE at component center", math.isclose(mse_at_center, 8.25, abs_tol=1e-10))

    # An asymmetric case prevents symmetry from hiding incorrect mixture weights.
    wa, ma, sa = [0.2, 0.3, 0.5], [-1.0, 0.5, 3.0], [0.25, 1.0, 0.75]
    amean, awithin, abetween = moments(wa, ma, sa)
    def asymmetric_density(x):
        return sum(wi * math.exp(-((x - mi) ** 2) / (2 * si * si))
                   / (si * math.sqrt(2 * math.pi)) for wi, mi, si in zip(wa, ma, sa))
    amass = simpson(asymmetric_density, -10, 12)
    amean_ref = simpson(lambda x: x * asymmetric_density(x), -10, 12)
    avar_ref = simpson(lambda x: (x - amean) ** 2 * asymmetric_density(x), -10, 12)
    check("asymmetric quadrature mass", math.isclose(amass, 1, abs_tol=1e-10))
    check("asymmetric quadrature mean", math.isclose(amean, amean_ref, abs_tol=1e-10))
    check("asymmetric quadrature variance", math.isclose(awithin + abetween, avar_ref, abs_tol=1e-10))
    check("asymmetric Decimal NLL", math.isclose(
        -log_mixture(1.7, wa, ma, sa), decimal_nll(1.7, wa, ma, sa), abs_tol=1e-10))

    references = [(x, -log_mixture(x, w, m, s), decimal_nll(x, w, m, s))
                  for x in [-8.0, -2.0, 0.0, 2.0, 8.0, 30.0]]
    for x, actual, reference in references:
        check("Decimal NLL x=" + str(x), math.isclose(actual, reference, rel_tol=0, abs_tol=1e-10))
    naive_far = density(30.0)
    stable_far_nll = -log_mixture(30.0, w, m, s)
    check("ordinary floating density underflows at x=30", naive_far == 0.0)
    check("log-domain NLL remains finite at x=30", math.isfinite(stable_far_nll))
    check("zero-weight component is skipped", math.isclose(
        log_mixture(1, [0, 1], [-99, 2], [1, 0.5]), log_normal(1, 2, 0.5)))
    check("single component reduces to Gaussian", math.isclose(
        log_mixture(0.3, [1], [0], [2]), log_normal(0.3, 0, 2)))
    invalid = [([], [], []), ([1], [0, 1], [1]), ([-0.1, 1.1], [0, 1], [1, 1]),
               ([0.4, 0.4], [0, 1], [1, 1]), ([0, 0], [0, 1], [1, 1]),
               ([1], [0], [0]), ([1], [0], [-1]), ([1], [math.inf], [1]),
               ([math.nan], [0], [1])]
    for i, params in enumerate(invalid):
        try:
            log_mixture(0, *params)
        except ValueError:
            check("reject invalid parameters " + str(i), True)
        else:
            check("reject invalid parameters " + str(i), False)
    for x in [math.nan, math.inf, -math.inf]:
        try:
            log_mixture(x, w, m, s)
        except ValueError:
            check("reject nonfinite observation " + str(x), True)
        else:
            check("reject nonfinite observation " + str(x), False)
    try:
        log_mixture(1e308, [1], [0], [1e-308])
    except ValueError:
        check("reject residual outside float range", True)
    else:
        check("reject residual outside float range", False)

    report = {
        "python_version": platform.python_version(), "checks": len(checks), "passed": checks,
        "mean": mean, "within_variance": within, "between_variance": between,
        "total_variance": variance, "density_at_mean": density_mean,
        "density_at_component_center": density_center, "mse_at_mean": mse_at_mean,
        "mse_at_component_center": mse_at_center, "mixture_nll_at_2": mixture_nll,
        "weighted_component_nll_at_2": component_nll, "naive_density_at_30": naive_far,
        "stable_nll_at_30": stable_far_nll, "decimal_references": references,
        "quadrature": {"range": [-8, 8], "intervals": 16000, "mass": mass},
        "asymmetric_reference": {"weights": wa, "means": ma, "sigmas": sa,
                                 "mean": amean, "variance": awithin + abetween,
                                 "integrated_mean": amean_ref, "integrated_variance": avar_ref},
        "scope": "Fixed one-dimensional Gaussian mixtures only; no training or real-data evaluation."
    }
    Path("report.json").write_text(json.dumps(report, indent=2) + "\n", encoding="utf-8")
    print(f"mean={mean:.6f} density_at_mean={density_mean:.9f} density_at_center={density_center:.9f}")
    print(f"MSE_at_mean={mse_at_mean:.6f} MSE_at_center={mse_at_center:.6f}")
    print(f"mixture_NLL={mixture_nll:.6f} weighted_component_NLL={component_nll:.6f}")
    print(f"variance: within={within:.6f} between={between:.6f} total={variance:.6f}")
    print(f"x=30 naive_density={naive_far:.1f} stable_NLL={stable_far_nll:.6f}")
    print(f"asymmetric_reference: mean={amean:.6f} variance={awithin + abetween:.6f}")
    print(f"checks={len(checks)} passed; references=Decimal+quadrature; report.json written")


if __name__ == "__main__":
    main()