It is this well implemented?

[PhyAMR]

Hi, I'm using PySR on some data and I want to check how the different models fit that data by calculating the reduced chi squared for each equation PySR returns and check how similar to 1 they are, and only get the best expresion of all.

def reduced_chi_squared(y_true, y_pred, uncertainty=None, dof=None):
    uncertainty = np.std(y_true) if uncertainty is None else uncertainty
    residuals = (y_true - y_pred) / uncertainty
    chi2 = np.sum(residuals ** 2)
    dof = len(y_true) - 1 if dof is None else dof
    return chi2 / dof

def symbolic_regression(df, input_cols, target_col):
    X = df[input_cols].values.reshape(-1, len(input_cols))
    y = df[target_col].values
    
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.3, random_state=12)
    
    model = PySRRegressor(
        niterations=30,
        binary_operators=["+", "-", "*", "/", "^"],
        unary_operators=["exp", "log", "sqrt"],
        model_selection="best",
        elementwise_loss="loss(x, y) = (x - y)^2",
        verbosity=0,
    )
    model.fit(X_train, y_train)
    
    best_expr = None
    best_chi2 = np.inf
    for i, row in model.equations_.iterrows():
        try:
            y_pred = model.predict(X_test, index=i)
            chi2 = reduced_chi_squared(y_test, y_pred)
            if abs(chi2 - 1) < abs(best_chi2 - 1):
                best_chi2 = chi2
                best_expr = model.sympy(i)
            
        except Exception as e:
            logger.warning(f"Error evaluating model {i}: {e}")
    if best_chi2 < 0.8 or chi2 > 1.2:
        logger.info(f"Model {i} has a chi2 value of {chi2}, which is outside the acceptable range.")
        return None, np.inf
    
    return best_expr, best_chi2

I want to know if this can be done by PySR internally (my guess is yes, maybe in the model selection variable I can define a function that makes this). Also I would like to have some outside opinion to know if this makes sense, mainly because I get strange things like this one

Here it has to be something linear, but slightly different, bu I still get value near 1 for the chi2 which is not possible (the prediction of PySR is the number above the chi2)

Thanks for reading