Skip to content

Instantly share code, notes, and snippets.

@korenmiklos
Created February 7, 2025 16:48
Show Gist options
  • Select an option

  • Save korenmiklos/fb33be8cf16862d92edad3d397ce3e8b to your computer and use it in GitHub Desktop.

Select an option

Save korenmiklos/fb33be8cf16862d92edad3d397ce3e8b to your computer and use it in GitHub Desktop.
using LogExpFunctions
using Distributions
function generate_data(K=1_000_000)
X = collect(range(-6, -3, length=K))
Y = rand.(Bernoulli.(logistic.(X)))
return X, Y
end
# dispatch on different types of numerical methods
logL(X, Y, β, ::Val{:direct}) = sum(Y .* log.(logistic.(X * β)) .+ (1 .- Y) .* log.(1 .- logistic.(X * β)))
logL(X, Y, β, ::Val{:loglogistic}) = sum(Y .* loglogistic.(X * β) .+ (1 .- Y) .* loglogistic.(-X * β))
logL(X, Y, β, ::Val{:split}) = sum(loglogistic.(X[Y] * β)) + sum(loglogistic.(-X[.!Y] * β))
function main()
X, Y = generate_data()
β = range(0, 2, length=1_000)
@time [logL(X, Y, b, Val(:direct)) for b in β]
@time [logL(X, Y, b, Val(:loglogistic)) for b in β]
@time [logL(X, Y, b, Val(:split)) for b in β]
end
main()
@korenmiklos
Copy link
Author

Don't worry about the type dispatch, this is just to reuse the same function name, a way to hijack Julia's multiple dispatch system. The important thing to watch is L#13. (But only if Y is a Bool, not an Int. L#6 shows how to generate random Bools.)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment