Рассчитать стандартизированные остатки Пирсона в Python

Question

Я хочу вычислить стандартизированные остатки Пирсона в Python (3.7.1), используя вывод scipy.stats.chi2_contingency.Я уже наткнулся на эту запись stackoverflow, и это именно то, что мне нужно, однако я получаю ошибочные результаты.Я могу только догадываться, что это может быть связано с моей более новой версией Python (ссылка с 2013 года)?

Я уже сломал формулу вычисления

v = csum * rsum * (n - rsum) * (n - csum) / n**3

в терминах cr_sum = csum * rsum и n_rcsum = (n - rsum) * (n - csum).Оба выходных массива имеют форму (2,5).Кажется, здесь необходимо вычислить произведение Адамара cr_sum и n_rcsum.Когда я сделал это вручную для первой ячейки (со значением частоты 33), я получил правильный остаток (-2.62309082).Однако я не могу заставить этот продукт Адамара работать в Python.Вместо этого Python, похоже, передает некоторые вещания и выводит:

array([[-1125512208, -267063340, -274153780, -1725637260, 691228240], [-1125512208, -267063340, -274153780, -1725637260, 691228240]]).

Более того, я обычно не понимаю, когда использовать какой тип умножения.В сообщении stackoverflow комментатор использовал только звездочку, и все, казалось, работало нормально.Какие изменения необходимо внести в код и почему?

Это мой код:

from __future__ import division

import numpy as np
from scipy.stats.contingency import margins
from scipy.stats import chi2_contingency

def residuals(observed, expected):
    return (observed - expected) / np.sqrt(expected)

def stdres(observed, expected):
    n = observed.sum()
    rsum, csum = margins(observed)
    v = csum * rsum * (n - rsum) * (n - csum) / n**3
    return (observed - expected) / np.sqrt(v)

F = np.array([[33, 250, 196, 136, 32], [55, 293, 190, 71, 13]])
chi2, p, dof, expected = chi2_contingency(F)
stdres = stdres(F,expected)

Answer 1

В Windows целочисленный тип по умолчанию для массива NumPy - 32-разрядный.Когда код в Что эквивалентно R data.chisq $ residuals в python? выполняется в Windows с входным массивом F = np.array([[33, 250, 196, 136, 32], [55, 293, 190, 71, 13]]), промежуточный расчет выражения csum * rsum * (n - rsum) * (n - csum) в функции stdres приводит к переполнению целого числа.Переполнение помещает отрицательные отрицательные значения в переменную v, поэтому при вычислении sqrt(v) вы получите nan с и предупреждение.

Исправление заключается в преобразовании rsum и csum вс плавающей запятой, прежде чем делать этот промежуточный расчет.Попробуйте эту версию:

def stdres(observed, expected):
    n = observed.sum()
    rsum, csum = margins(observed)
    rsum = rsum.astype(np.float64)
    csum = csum.astype(np.float64)
    v = csum * rsum * (n - rsum) * (n - csum) / n**3
    return (observed - expected) / np.sqrt(v)