So yes, your implementation was helpful for finding the two bugs. First is the Bessel’s correction(i.e. diving by n-1 in the variance formula), due to unbiased=True.
Second is the dimensions I am getting the mean and variance from. It should have been
x.mean([0, 2, 3])
and
x.var([0, 2, 3], unbiased=False)