JuMPのCox比例ハザードのスコア関数

Question

私は、方程式をコードに変換する経験はあまりありません。私はJuuliaでコード化する部分尤度のスコア関数をJuMPで評価するように変換することに固執しました。

score function 0でベータのために解くと、

小さな単純なデータセットを作成しました。

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の符号を選択​​するという共通のトリックを使用します。