How to make this AFT regression stable?

Specific Domains Statistics
, , ,
6

Fixed a bunch of bugs, thanks for all the suggestions.
The following seems to run fine, but has identifiability issues. While the relative magnitude of positive beta seems roughly correct, it can’t get the sign right.
I tried adding a ridge regularizer only on the linear coefficients (commented) but that led to stability issues as well

using Revise,Distributions,Optim,DistributedTopicModels,Survival,LineSearches,Plots,Random,LinearAlgebra

n = 10000
# beta = [1,10]
beta = [10,1]
# beta = [1,-1]
# beta = [-1,1]
s = Random.seed!(1)
X = rand(s,n,2)
delta = rand([0,1],n) # delta[i] = 1 iff Y[i] censored
μ = X * beta

histogram(μ)

λ = 1.1
α = .1
wd = Weibull(λ,α)
s = Random.seed!(1)
ϵ = rand(s,wd,n)

Y = log.(0.1(μ) + ϵ)
# Y = log.(0.1(μ) + ϵ)*0.1

histogram(Y)

Ycen = [delta[i] == 1 ? Y[i]-rand()*Y[i] : Y[i] for i in 1:n]
paramInit = fill(0.1,3)

function tAFT(pvec,data)
    function f(w)
        exp(w - exp(w))
    end
    function F(w)
        exp(-exp(w))
    end
    function l(theta)
        -sum(log(f((d[2]-d[3:end]'*theta[2:end])/theta[1])^d[1]) + log(F((d[2]-d[3:end]'*theta[2:end])/theta[1])^(1-d[1])) for d in data)
        # -sum(log(f((d[2]-d[3:end]'*theta[2:end])/theta[1])^d[1]) + log(F((d[2]-d[3:end]'*theta[2:end])/theta[1])^(1-d[1])) + 0.1*norm(theta[3:end])^2 for d in data)
    end
    optimize(l,pvec,BFGS(linesearch=LineSearches.BackTracking()),autodiff = :forward)
    # optimize(l,pvec,LBFGS(),autodiff = :forward)
end


resp = tAFT(paramInit, collect(eachrow([delta Y X])))
resp.minimizer

julia> resp = tAFT(paramInit, collect(eachrow([delta Y X])))
 * Status: success

 * Candidate solution
    Final objective value:     8.489688e+03

 * Found with
    Algorithm:     BFGS

 * Convergence measures
    |x - x'|               = 4.88e-15 ≰ 0.0e+00
    |x - x'|/|x'|          = 1.87e-15 ≰ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 7.78e-08 ≰ 1.0e-08

 * Work counters
    Seconds run:   1  (vs limit Inf)
    Iterations:    76
    f(x) calls:    150
    ∇f(x) calls:   77


julia> resp.minimizer
3-element Array{Float64,1}:
 2.6089785118820537
 2.0258941072666534
 0.362426844935407