PySR for Cumulative Distribution Functions

[ghost]

Hello @MilesCranmer and PySR community,

I am using PySR to fit cumulative distribution functions (CDFs) which requires me to impose some constraints on the fit. Firstly, predicted probabilities must lie between 0 and 1, and secondly, the function must always be increasing.

Could you please offer some advice on how to enforce these constraints within the algorithm? I think the best way to do this might be to create a custom loss function in Julia but I am less familiar with the Julia programming language.

Any advice would be greatly appreciated!

Many thanks,
Strath

[MilesCranmer]

Maybe like this? It only checks the points you pass for monotonicity though; maybe that's enough? I guess you could always pass more points if needed.

from pysr import jl, PySRRegressor

jl.seval("using Zygote")  # Need to load package for gradient calculations

loss_function = """
function cdf_loss(tree::Node, dataset::Dataset{T,L}, options::Options) where {T,L}
    X = dataset.X
    y = dataset.y

    y_pred, grad, completed = eval_grad_tree_array(tree, X, options; variable=true)

    if !completed
        return convert(L, 1e9)
    end

    loss = convert(L, sum((y_pred .- y).^2) / length(y))
    is_monotonic = all(grad .>= 0)
    min_val = minimum(y_pred)
    max_val = maximum(y_pred)

    if !is_monotonic
        # Add some penalty:
        loss += convert(L, 1e3)
    end

    if min_val < 0 || max_val > 1
        # Add some penalty for not being a CDF:
        loss += convert(L, 1e3)
    end

    return loss
end
"""

model = PySRRegressor(
    binary_operators=["+", "-", "*", "/"],
    unary_operators=["exp", "square"],
    loss_function=loss_function,
    enable_autodiff=True,
)

With Julia syntax highlighting:

function cdf_loss(tree::Node, dataset::Dataset{T,L}, options::Options) where {T,L}
    X = dataset.X
    y = dataset.y

    y_pred, grad, completed = eval_grad_tree_array(tree, X, options)

    if !completed
        return convert(L, 1e9)
    end

    loss = sum((y_pred .- y).^2) / length(y)
    is_monotonic = all(grad .>= 0)
    min_val = minimum(y_pred)
    max_val = maximum(y_pred)

    if !is_monotonic
        # Add some penalty:
        loss += convert(L, 1e3)
    end

    if min_val < 0 || max_val > 1
        # Add some penalty for not being a CDF:
        loss += convert(L, 1e3)
    end

    return convert(L, loss)
end

[ghost]

Thanks so much, @MilesCranmer!

I have tested this on several data and it seems to work well. Much appreciated!

Strath