Logistic Regression Mixture Model

[carterjgreen]

I posted this on the discourse but it hasn't gained any traction and I'm not sure if it is an issue or just my implementation.

I'm trying to implement robust logistic regression fromchapter 10.5 of PML. The model is
$p(y|\mathbf{x}) = \pi Ber(y|0.5) + (1-\pi) Ber(y|\sigma(\mathbf{w}^\top\mathbf{x}))$
where $\pi$ are the mixture weights drawn from a uniform distribution and the regression weights are drawn from a normal distribution with mean 0 and std 10.

The author has a pymc3 model here but trying to implement it with any other sampler but MH() blows up with numerical errors. Does it have to do with the logistic function? Any help would be appreciated.

using Turing, MCMCChains, StatsPlots, Random, StatsFuns
using RDatasets

iris = dataset("datasets", "iris")
select!(iris, :SepalLength, :Species)
filter!(row -> row[:Species] in ["setosa", "versicolor"], iris)
iris[!, :target] = iris[!, :Species] .!= "setosa"
select!(iris, :SepalLength, :target);

log_x = iris[!, :SepalLength]
log_y = iris[!, :target]
log_xoutliers = [4.2, 4.5, 4.0, 4.3, 4.2, 4.4]
log_youtliers = ones(length(log_xoutliers))

log_x_train = [log_x; log_xoutliers]
log_y_train = [log_y; log_youtliers];

@model function robust_logreg(x, y)
    N = length(x)
    alpha ~ Normal(0, 10)
    beta ~ Normal(0, 10)
    
    mu  = alpha .+ beta .* x
    theta = logistic.(mu)
    
    pie ~ filldist(Beta(1, 1), N) # probability of contamination of outliers
    p = pie .* 0.5 .+ (1 .- pie) .* theta
    
    y .~ Bernoulli.(p)
end

roblog_model = robust_logreg(log_x_train, log_x_train)
roblog_chain = sample(roblog_model, HMC(0.05, 10), 1000)

[luiarthur]

I'm a little rusty with pymc3, but if the model you described above is what you're trying to implement, it should be this:

@model function robust_logreg(x, y)
    N = length(x)
    alpha ~ Normal(0, 10)
    beta ~ Normal(0, 10)
    
    mu  = alpha .+ beta * x  # Note that a second dot is not needed
    theta = logistic.(mu)
    
    pie ~ filldist(Beta(1, 1), N) # probability of contamination of outliers
    
    # Can't broadcast `MixtureModel` the way we want to, so we specify the 
    # the distribution of each `y[i]` via a loop. This specification is clear, but 
    # not very efficient.
    for i in 1:N
        y[i] ~ MixtureModel([Bernoulli(.5), Bernoulli(theta[i])], [pie[i], 1 - pie[i]])
    end
end

When N is large, however, the loop might be slow. The remedy would be to directly increment the log unnormalized joint posterior using Turing.@addlogprob!. This would allow you to broadcast, which in this case would be more efficient. So the last three lines would be replaced with:

        @Turing.@addlogprob! logaddexp.(log.(pie) .+ N * log(0.5), log1p.(-pie) + binomlogpdf.(1, theta, y))

Note that logaddexp and binomlogpdf are from the StatsFuns julia package.

[devmotion]

but if the model you described above is what you're trying to implement, it should be this:

There is no need to use a (less efficient) MixtureModel here - the mixture model of two Bernoulli distributions is just another Bernoulli distribution. The definition of the distribution of y in the OP is correct.

    mu  = alpha .+ beta * x  # Note that a second dot is not needed

While the second dot is not needed, it ensures that no intermediate arrays are created and all computations are fused. I.e., the implementation in the OP with additional dot is more efficient.

---

One problem seems to be that pie should be an univariate random variable, i.e., it should be the same for all data points. This would match both the model description and the PyMC3 implementation. Additionally, instead of Beta(1, 1) one can use the (probably more efficient) Uniform(). The main issue, however, seems to be a typo in the definition of roblog_model: you used log_x_train instead of log_y_train - this leads to numerical errors in the adaptation phase (all samples are rejected) since the values of the log pdf are non-fnite (you can inspect them by adding @show Turing.getlogp(_varinfo) at the end of the model).

Taken together, the model can be implemented as:

using RDatasets
using StatsFuns
using Turing

iris = dataset("datasets", "iris")
filter!(:Species => !=("virginica"), iris)
iris[!, :target] = iris.Species .== "versicolor"

log_x = iris.SepalLength
log_y = iris.target
log_xoutliers = [4.2, 4.5, 4.0, 4.3, 4.2, 4.4]
log_youtliers = trues(length(log_xoutliers))

log_x_train = [log_x; log_xoutliers]
log_y_train = [log_y; log_youtliers]

@model function robust_logreg(x, y)
    alpha ~ Normal(0, 10)
    beta ~ Normal(0, 10)
    # probability of contamination of outliers
    pie ~ Uniform() # alternatively (but probably less efficient): Beta(1, 1)
    
    y .~ Bernoulli.(pie * 0.5 .+ (1 - pie) .* logistic.(alpha .+ beta .* x))
end

roblog_model = robust_logreg(log_x_train, log_y_train)
roblog_chain = sample(roblog_model, NUTS(), 1000)

[carterjgreen]

Thank you both for your replies! I've learned something from both. My main problem was that I was using the wrong data. I also incorrectly assumed that the original model assigned a contamination probability for each point. I'm now seeing the results that I expected.