Repository: scipy/scipy
Pull request: MAINT: stats: chisquare check sum of observed/expected frequencies
#### Reference issue
Closes gh-12282

#### What does this implement/fix?
There was no check in `scipy.stats.chisquare` that the sum of observed frequencies matched the sum of expected frequencies. This PR raises an error if the two do not agree within a certain tolerance. 

#### Additional information
The side effect is that it is no longer possible to express the expected frequencies in terms of probabilities; we can add automatic normalization of the observed frequencies for this case if desired.
A question is whether this should be specific to `chisquare` or is more generally applicable in `power_divergence`. Currently the check is performed in `power_divergence` whenever `lambda_ == 1` (Pearson).

@sam-data-guy would this address the concern?

Problem statement:
API: stats.chisquare: must observed and expected frequencies add to same total?
<!-- 
Thank you for taking the time to file a bug report. 
Please fill in the fields below, deleting the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
Note: This is a comment, and won't appear in the output.
-->

The function scipy.stats.chisquare performs the chi-square test on a vector of observed and expected frequencies.  For the test to make sens, the observed and expected frequency vectors must sum to the same total, otherwise the result is nonsense since the inputs are incompatible.  So there are two options: either an error should be thrown if they do not, or f_exp need to be rescaled so they sum to the total of f_obs.  The first is probably better, or the second with a warning thrown.

#### Reproducing code example:


```
# In the following example, one vector is an exact multiple of the other.
# this means that the observed and expected frequencies are exactly
# proportional.  This should give a p-value of 1 (not significant at all).  
# Instead, you get the following:
from scipy.stats import chisquare
chisquare(f_obs=[10,20], f_exp=[30,60])

# Power_divergenceResult(statistic=40.0, pvalue=2.5396285894708634e-10)
# the statistic of 40 is calculated as following, which directly follows the 
# formula.  This is then plugged into a chi squared distribution to get a 
#p-value of close to 0, which is the opposite of the significance you should
# get.
# ((10-30)**2 /30) + ((20-60)** 2 / 60) = 40

# Instead, here is what should happen
import numpy as np
fobs = np.array([10,20])
fexp = np.array([30,60])
# adjust the totals
# gives array([10., 20.]), the same as observed
fexp = fexp * (np.sum(fobs)/np.sum(fexp)) 
chisquare(f_obs=fobs, f_exp=fexp)
# the correct result
# Power_divergenceResult(statistic=0.0, pvalue=1.0)
```


#### Scipy/Numpy/Python version information:
```
1.4.1 1.18.4 sys.version_info(major=3, minor=7, micro=4, releaselevel='final', serial=0)
```

Hints from issue discussion:
The two proposed changes sound reasonable to me and I can implement either. A third option, if this is a "feature", would be to change nothing about the way the test functions but to document why the test is performed the way it is. I don't have an informed opinion. 
As I said, the test output is nonsensical since the inputs are incompatible because the test requires them to sum to the same total.  It seems to me that the fact it has been programmed this way is a defect from the beginning.   The third option you note amounts to simply documenting the fact that the test does not work as expected, and leaving the user to discover the mistake.  I mean, it does do a chi square calculation but if the inputs are not compatible it is not a real test.  I think having a warning or note thrown is best, if you are concerned about backwards compatibility.  For instance, the anderson_darling test throws a note about p-value truncation (which I frankly think is silly as well, but what the heck).  But I really think you should think about correcting it in some way (perhaps deprecating the function and using a new one), otherwise you run the risk of perpetuating the error in the future.  Furthermore, it shouldn't be up to the user to make sure that their inputs are legal, particularly since not everyone knows better.
AFAIR, what I do in statsmodels in a similar case is raise an exception if sums differs by more than a numeric small number, or rescale if sums differ but `expected` has values between 0 and 1.
The second case allows that expected values can be defined as probabilities or relative frequencies.

Different sums of integer values is ambiguous and rescaling cannot be done in a safe way, i.e. we don't know which to rescale, unless we want to second guess which mistake the user made or what the user really meant.

Disclaimer: I haven't looked at the scipy functions in a long time.

Commit message: MAINT: stats: redefine relative difference in a more standard way

Patch to review:
diff --git a/scipy/stats/stats.py b/scipy/stats/stats.py
index 79e0930bad11..531724cf8d80 100644
--- a/scipy/stats/stats.py
+++ b/scipy/stats/stats.py
@@ -5966,6 +5966,13 @@ def _count(a, axis=None):
     return num
 
 
+def _m_broadcast_to(a, shape):
+    if np.ma.isMaskedArray(a):
+        return np.ma.masked_array(np.broadcast_to(a, shape),
+                                  mask=np.broadcast_to(a.mask, shape))
+    return np.broadcast_to(a, shape, subok=True)
+
+
 Power_divergenceResult = namedtuple('Power_divergenceResult',
                                     ('statistic', 'pvalue'))
 
@@ -6030,6 +6037,10 @@ def power_divergence(f_obs, f_exp=None, ddof=0, axis=0, lambda_=None):
     category are too small.  A typical rule is that all of the observed
     and expected frequencies should be at least 5.
 
+    Also, the sum of the observed and expected frequencies must be the same
+    for the test to be valid; `power_divergence` raises an error if the sums
+    do not agree within a relative tolerance of ``1e-8``.
+
     When `lambda_` is less than zero, the formula for the statistic involves
     dividing by `f_obs`, so a warning or error may be generated if any value
     in `f_obs` is 0.
@@ -6134,9 +6145,25 @@ def power_divergence(f_obs, f_exp=None, ddof=0, axis=0, lambda_=None):
         lambda_ = 1
 
     f_obs = np.asanyarray(f_obs)
+    f_obs_float = f_obs.astype(np.float64)
 
     if f_exp is not None:
-        f_exp = np.asanyarray(f_exp)
+        f_exp = _m_broadcast_to(np.asanyarray(f_exp), f_obs_float.shape)
+        rtol = 1e-8  # to pass existing tests
+        with np.errstate(invalid='ignore'):
+            f_obs_sum = f_obs_float.sum(axis=axis)
+            f_exp_sum = f_exp.sum(axis=axis)
+            relative_diff = (np.abs(f_obs_sum - f_exp_sum) /
+                             np.minimum(f_obs_sum, f_exp_sum))
+            diff_gt_tol = (relative_diff > rtol).any()
+        if diff_gt_tol:
+            msg = (f"For each axis slice, the sum of the observed "
+                   f"frequencies must agree with the sum of the "
+                   f"expected frequencies to a relative tolerance "
+                   f"of {rtol}, but the percent differences are:\n"
+                   f"{relative_diff}")
+            raise ValueError(msg)
+
     else:
         # Ignore 'invalid' errors so the edge case of a data set with length 0
         # is handled without spurious warnings.
@@ -6148,7 +6175,7 @@ def power_divergence(f_obs, f_exp=None, ddof=0, axis=0, lambda_=None):
     # cases of lambda_.
     if lambda_ == 1:
         # Pearson's chi-squared statistic
-        terms = (f_obs.astype(np.float64) - f_exp)**2 / f_exp
+        terms = (f_obs_float - f_exp)**2 / f_exp
     elif lambda_ == 0:
         # Log-likelihood ratio (i.e. G-test)
         terms = 2.0 * special.xlogy(f_obs, f_obs / f_exp)
@@ -6213,6 +6240,10 @@ def chisquare(f_obs, f_exp=None, ddof=0, axis=0):
     category are too small.  A typical rule is that all of the observed
     and expected frequencies should be at least 5.
 
+    Also, the sum of the observed and expected frequencies must be the same
+    for the test to be valid; `chisquare` raises an error if the sums do not
+    agree within a relative tolerance of ``1e-8``.
+
     The default degrees of freedom, k-1, are for the case when no parameters
     of the distribution are estimated. If p parameters are estimated by
     efficient maximum likelihood then the correct degrees of freedom are
diff --git a/scipy/stats/tests/test_stats.py b/scipy/stats/tests/test_stats.py
index 5bdd7d37ccd3..d5133a1fa270 100644
--- a/scipy/stats/tests/test_stats.py
+++ b/scipy/stats/tests/test_stats.py
@@ -2979,6 +2979,26 @@ def test_power_divergence_result_attributes(self):
         attributes = ('statistic', 'pvalue')
         check_named_results(res, attributes)
 
+    def test_power_divergence_gh_12282(self):
+        # The sums of observed and expected frequencies must match
+        f_obs = np.array([[10, 20], [30, 20]])
+        f_exp = np.array([[5, 15], [35, 25]])
+        with assert_raises(ValueError, match='For each axis slice...'):
+            stats.power_divergence(f_obs=[10, 20], f_exp=[30, 60])
+        with assert_raises(ValueError, match='For each axis slice...'):
+            stats.power_divergence(f_obs=f_obs, f_exp=f_exp, axis=1)
+        stat, pval = stats.power_divergence(f_obs=f_obs, f_exp=f_exp)
+        assert_allclose(stat, [5.71428571, 2.66666667])
+        assert_allclose(pval, [0.01682741, 0.10247043])
+
+
+def test_gh_chisquare_12282():
+    # Currently `chisquare` is implemented via power_divergence
+    # in case that ever changes, perform a basic test like
+    # test_power_divergence_gh_12282
+    with assert_raises(ValueError, match='For each axis slice...'):
+        stats.chisquare(f_obs=[10, 20], f_exp=[30, 60])
+
 
 @pytest.mark.parametrize("n, dtype", [(200, np.uint8), (1000000, np.int32)])
 def test_chiquare_data_types(n, dtype):
@@ -3070,28 +3090,18 @@ def test_power_divergence_against_cressie_read_data():
     # J. R. Statist. Soc. B (1984), Vol 46, No. 3, pp. 440-464.
     # This tests the calculation for several values of lambda.
 
+    # Table 4 data recalculated for greater precision according to:
+    # Shelby J. Haberman, Analysis of Qualitative Data: Volume 1
+    # Introductory Topics, Academic Press, New York, USA (1978).
+    obs = np.array([15, 11, 14, 17, 5, 11, 10, 4, 8, 10, 7, 9, 11, 3, 6, 1, 1, 4])
+    beta = -0.083769  # Haberman (1978), p. 15
+    i = np.arange(1, len(obs) + 1)
+    alpha = np.log(obs.sum() / np.exp(beta*i).sum())
+    expected_counts = np.exp(alpha + beta*i)
+
     # `table4` holds just the second and third columns from Table 4.
-    table4 = np.array([
-        # observed, expected,
-        15, 15.171,
-        11, 13.952,
-        14, 12.831,
-        17, 11.800,
-        5, 10.852,
-        11, 9.9796,
-        10, 9.1777,
-        4, 8.4402,
-        8, 7.7620,
-        10, 7.1383,
-        7, 6.5647,
-        9, 6.0371,
-        11, 5.5520,
-        3, 5.1059,
-        6, 4.6956,
-        1, 4.3183,
-        1, 3.9713,
-        4, 3.6522,
-        ]).reshape(-1, 2)
+    table4 = np.vstack((obs, expected_counts)).T
+
     table5 = np.array([
         # lambda, statistic
         -10.0, 72.2e3,


Code review comments:
diff_hunk: @@ -6134,9 +6145,25 @@ def power_divergence(f_obs, f_exp=None, ddof=0, axis=0, lambda_=None):
         lambda_ = 1
 
     f_obs = np.asanyarray(f_obs)
+    f_obs_float = f_obs.astype(np.float64)
 
     if f_exp is not None:
-        f_exp = np.asanyarray(f_exp)
+        f_exp = _m_broadcast_to(np.asanyarray(f_exp), f_obs_float.shape)
path: scipy/stats/stats.py
text: @user1:
With this change, the test and the refguide-check pass:
```suggestion
        f_exp = np.asanyarray(f_exp)
        bshape = _broadcast_shapes(f_obs_float.shape, f_exp.shape)
        f_obs_float = _m_broadcast_to(f_obs_float, bshape)
        f_exp = _m_broadcast_to(f_exp, bshape)
```
This figures out the common broadcast shape, and then uses `_m_broadcast_to()` on both inputs.

`_broadcast_shapes()` is already defined here in `stats.py`.

What do you think?

@author:
It makes me cry a little.

> _broadcast_shapes() is already defined here in stats.py.
```
np.broadcast(x, y).shape
```

But I'll commit it.


@user1:
Looks like that would work fine too, if you prefer it.

diff_hunk: @@ -3070,28 +3090,18 @@ def test_power_divergence_against_cressie_read_data():
     # J. R. Statist. Soc. B (1984), Vol 46, No. 3, pp. 440-464.
     # This tests the calculation for several values of lambda.
 
+    # Table 4 data recalculated for greater precision according to:
+    # Shelby J. Haberman, Analysis of Qualitative Data: Volume 1
+    # Introductory Topics, Academic Press, New York, USA (1978).
+    obs = np.array([15, 11, 14, 17, 5, 11, 10, 4, 8, 10, 7, 9, 11, 3, 6, 1, 1, 4])
path: scipy/stats/tests/test_stats.py
text: @user1:
```suggestion
    obs = np.array([15, 11, 14, 17, 5, 11, 10, 4, 8,
                    10, 7, 9, 11, 3, 6, 1, 1, 4])
```