Mathematical argument
Convention: a random variable is "normal" if it is N(μ, σ²) with σ² > 0, or a constant (σ² = 0), per proposition A.
Proposition B — TRUE
Let X ~ N(μ₁, σ₁²), Y ~ N(μ₂, σ₂²) be independent. Characteristic functions: φ_X(t) = exp(iμ₁t − σ₁²t²/2), φ_Y(t) = exp(iμ₂t − σ₂²t²/2). By independence, φ_{X+Y}(t) = φ_X(t)φ_Y(t) = exp(i(μ₁+μ₂)t − (σ₁²+σ₂²)t²/2), which is the characteristic function of N(μ₁+μ₂, σ₁²+σ₂²). By the uniqueness theorem for characteristic functions, X+Y ~ N(μ₁+μ₂, σ₁²+σ₂²) (degenerate cases included). ∎ (More generally the result holds whenever (X, Y) is jointly Gaussian; independence of normal marginals is a special case.)
Proposition A — FALSE
Normal marginals do not determine the joint law; only joint normality guarantees normal sums.
Note on the trivial candidate. Y = −X gives X+Y ≡ 0, which is a degenerate normal under the stated convention, so it is not a counterexample. We need a sum that is neither constant nor normal.
Construction. Let X ~ N(0,1) and let S be independent of X with P(S = 1) = P(S = −1) = 1/2. Define Y = S·X.
- X is normal: by definition.
- Y is normal: for any Borel B, P(Y ∈ B) = ½P(X ∈ B) + ½P(−X ∈ B) = P(X ∈ B), since −X has the same law as X (symmetry of N(0,1)). So Y ~ N(0,1).
- X+Y is not normal: X+Y = (1+S)X, which equals 0 when S = −1 and 2X when S = 1. Hence P(X+Y = 0) = P(S = −1) = 1/2 (plus P(X=0)·½ = 0).
- A non-degenerate normal is absolutely continuous, so it gives probability 0 to every point; X+Y has an atom of mass 1/2 at 0 — not a non-degenerate normal.
- X+Y is not constant: P(X+Y = 0) = 1/2 < 1, and on {S=1} it equals 2X, which is continuous.
- Hence X+Y is not normal under any convention. Its law is ½δ₀ + ½N(0,4); E[X+Y]=0, Var(X+Y)=2.
- Remark: Cov(X,Y) = E[S]E[X²] = 0, so X and Y are uncorrelated yet not independent (|Y| = |X|), and the pair (X, Y) is not jointly Gaussian (if it were, uncorrelated would imply independent and X+Y would be N(0,2)).
Numerical illustration (seed and n from data/design.csv; out/counterexample.csv, out/analysis.py): x and y have KS distances to fitted normals of about 0.002 (at the level of sampling noise), while x+y has KS distance ≈ 0.25 to its fitted normal and ≈ 50% exact zeros. See out/results.json. The simulation illustrates the argument; the proof is the argument above.