分布差异的平方，为什么会算出负数？

"""Exact finite enumeration for biased and unbiased MMD squared."""
import itertools
import json
import math
from pathlib import Path


def mmd2(x,y,sigma=1.0,unbiased=True):
    if not math.isfinite(sigma) or sigma <= 0:
        raise ValueError('positive finite bandwidth required')
    if len(x)<(2 if unbiased else 1) or len(y)<(2 if unbiased else 1):
        raise ValueError('not enough observations')
    if not all(math.isfinite(v) for v in list(x)+list(y)):
        raise ValueError('finite observations required')
    def k(a,b): return math.exp(-0.5*((a-b)/sigma)**2)
    def within(values):
        n=len(values)
        return sum(k(a,b) for i,a in enumerate(values) for j,b in enumerate(values)
                   if not unbiased or i!=j)/(n*(n-1) if unbiased else n*n)
    return within(x)+within(y)-2*sum(k(a,b) for a in x for b in y)/(len(x)*len(y))


def main():
    x=y=[0.0,1.0]
    u=mmd2(x,y); v=mmd2(x,y,unbiased=False)
    k=math.exp(-0.5)
    assert abs(u-(k-1))<1e-12 and abs(v)<1e-12
    # Each four-bit sequence has probability 1/16: independent iid Bernoulli samples.
    outcomes=[]
    for bits in itertools.product([0.0,1.0],repeat=4):
        a,b=bits[:2],bits[2:]
        outcomes.append(dict(x=a,y=b,u=mmd2(a,b),v=mmd2(a,b,unbiased=False)))
    mean_u=sum(r['u'] for r in outcomes)/16
    mean_v=sum(r['v'] for r in outcomes)/16
    clipped_mean=sum(max(0,r['u']) for r in outcomes)/16
    assert abs(mean_u)<1e-12 and abs(mean_v-(1-k)/2)<1e-12
    assert clipped_mean > 0
    assert all(r['v'] >= -1e-12 for r in outcomes)
    assert abs(mmd2([0,2],[1,3],2)-mmd2([0,20],[10,30],20))<1e-12
    assert abs(mmd2([0,2],[1,3])-mmd2([2,0],[3,1]))<1e-12
    rejected=0
    for a,b,s in [([0],[1],1),([0,1],[0,1],0),([0,1],[0,1],float('nan'))]:
        try: mmd2(a,b,s)
        except ValueError: rejected+=1
    assert rejected==3
    report=dict(example_u=u,example_v=v,population_mmd2=0,mean_u=mean_u,mean_v=mean_v,
        clipped_mean=clipped_mean,exact_outcomes=outcomes,invalid_rejections=rejected,checks_passed=True)
    Path('report.json').write_text(json.dumps(report,indent=2)+'\n',encoding='utf-8')
    print(f'same observed lists: U={u:.6f}; V={v:.6f}')
    print(f'16 equally likely iid outcomes: E[U]={mean_u:.6f}; E[V]={mean_v:.6f}')
    print(f'E[max(0,U)]={clipped_mean:.6f}')
    print('checks=PASS; report.json written')


if __name__=='__main__': main()