I can't find the information how $\chi^2$ are used to select numerical features for a model.

Fro instance, If I employ the sklearn library:

from sklearn.datasets import load_iris
from sklearn.feature_selection import chi2

iris = load_iris()
X, y = iris.data, iris.target # X contains 4 features, y does 3 classes

I can build the following table for the dataset

pivot_table = np.zeros((X.shape[1], len(np.unique(y))))
pivot_table[:, 0] = X[y==0].sum(axis=0)
pivot_table[:, 1] = X[y==1].sum(axis=0)
pivot_table[:, 2] = X[y==2].sum(axis=0)

print(pivot_table.T)
array([[250.3, 171.4,  73.1,  12.3],
       [296.8, 138.5, 213. ,  66.3],
       [329.4, 148.7, 277.6, 101.3]])

Then applying the contingency table to it:

from scipy.stats import chi2_contingency
chi2_contingency(pivot_table, correction=False)

Output:

125.58397740261773
1.0924329599765022e-24
6
[[213.82265358 301.31664021 361.36070621]
 [111.8757204  157.6540915  189.0701881 ]
 [137.51492279 193.78458652 232.40049069]
 [ 43.88670323  61.84468177  74.168615  ]]

But this only provides the information about dependence of variables. And they are dependent because p-value is quite low.

Then I use the other function and I get what I want:

chi2(X,y)
(array([ 10.81782088,   3.7107283 , 116.31261309,  67.0483602 ]),
 array([4.47651499e-03, 1.56395980e-01, 5.53397228e-26, 2.75824965e-15]))

But how does this work? What math lies at the heart of it?