JuMP中Cox比例風險的得分函數

Question

我在將方程轉換爲代碼方面經驗不足。我被卡在將Juul中部分可能性的分數函數轉換爲代碼以在JuMP中進行評估。

score function其中，當在0解決beta時，是最大值。

我做了一個簡單的小數據集。

Using DataFrames, DataFramesMeta, JuMP, Ipopt 
#build DataFrame 
times = [6,7,10,15,19,25] 
is_censored = [1,0,1,1,0,1] 
x= is_control = [1,1,0,1,0,0] 

m = Model(solver=IpoptSolver(print_level=0)) 
using DataFrames 

df = DataFrame(); 
df[:times]=times; 
df[:is_censored]= is_censored; 
df[:x]=x; 
df 

#sort df 
df_sorted = sort!(df, cols = [order(:times)]) 

#make df_risk and df_uncensored 
df_uncensored = @where(df_sorted, :is_censored .== 0) 
df_risk = df_sorted 

#use JuMP 

##convert df to array 

uncensored = convert(Array,df_uncensored[:x]) 
risk_set = convert(Array,df_risk[:x]) 
risk_index = convert(Array,find(is_censored .== 0)) 
x = convert(Array, x) 
@variable(m, β, start = 0.0) 

# score 
@NLobjective(m, Max, sum(uncensored[i] - (([sum(exp(risk_set[j]*β)*x[j]) for j = risk_index[i]:length(risk_set)])/([(sum(exp(risk_set[j]*β)*x[j])) for j=risk_index[i]:length(risk_set)])) for i = 1:length(uncensored))) 

我得到的錯誤是

ERROR: exp is not defined for type AffExpr. Are you trying to build a nonlinear problem? Make sure you use @NLconstraint/@NLobjective. 
Stacktrace: 
[1] exp(::JuMP.GenericAffExpr{Float64,JuMP.Variable}) at /home/icarus/.julia/v0.6/JuMP/src/operators.jl:630 
[2] collect(::Base.Generator{UnitRange{Int64},##58#60}) at ./array.jl:418 
[3] macro expansion at /home/icarus/.julia/v0.6/JuMP/src/parseExpr_staged.jl:489 [inlined] 
[4] macro expansion at /home/icarus/.julia/v0.6/JuMP/src/parsenlp.jl:226 [inlined] 
[5] macro expansion at /home/icarus/.julia/v0.6/JuMP/src/macros.jl:1086 [inlined] 
[6] anonymous at ./<missing>:? 

錯誤說，指數的問題，但我以前一樣具有指數在它的數似然和有沒有錯誤。

在R中的測試= -1.3261

如果上述代碼是工作我期望運行

solve(m) 

println("β = ", getvalue(β)) 

Answer 1

的代碼中的問題後相同的輸出是一個位繞口，但以下是試圖提取相關部件用於工作方法：

using JuMP, Ipopt 

times = [6,7,10,15,19,25]; 
is_censored = 1-[1,0,1,1,0,1]; 
is_control = 1-[1,1,0,1,0,0]; 

uncensored = find(is_censored .== 0) 

println("times = $times") 
println("is_censored = $is_censored") 
println("is_control = $is_control") 

m = Model(solver=IpoptSolver(print_level=0)) 
@variable(m, β, start = 0.0) 
@NLobjective(m, Max, sum(log(1+(-1)^is_control[uncensored[i]]* 
    sum((-1)^is_control[j]*exp(is_control[j]*β) for j=uncensored[i]:length(times))/ 
    sum(exp(is_control[j]*β) for j=uncensored[i]:length(times))) 
    for i=1:length(uncensored))) 

solve(m) 
println("β = ", getvalue(β)) 

此輸出：

times = [6,7,10,15,19,25] 
is_censored = [0,1,0,0,1,0] 
is_control = [0,0,1,0,1,1] 

β = -1.3261290591982942 

β是同一個問題，所以我猜想輸入的調整是正確的，公式是對數似然。日誌內部表達式的開始使用了根據0/1值bool和(-1)^bool選擇+1或-1符號的常用技巧。