Julia JuMP Многомерная оценка машинного обучения

Я пытаюсь выполнить ML-оценку нормально распределенной переменной в настройке линейной регрессии в Julia с использованием JuMP и решателя NLopt.

здесь существует хороший рабочий пример, однако, если я попробую для оценки параметров регрессии (наклона) код становится довольно утомительным для написания, особенно если пространство параметров увеличивается.

Может у кого-то есть идея, как это написать более лаконично. Вот мой код:

#type definition to store data
type data
    n::Int
    A::Matrix
    β::Vector
    y::Vector
    ls::Vector
    err::Vector
end

#generate regression data
function Data( n = 1000 )
    A = [ones(n) rand(n, 2)]
    β = [2.1, 12.9, 3.7]
    y = A*β + rand(Normal(), n)
    ls = inv(A'A)A'y
    err = y - A * ls
    data(n, A, β, y, ls, err)
end

#initialize data
d = Data()
println( var(d.y) )

function ml(  )
    m = Model( solver = NLoptSolver( algorithm = :LD_LBFGS ) )
    @defVar( m, b[1:3] )
    @defVar( m, σ >= 0, start = 1.0 )

    #this is the working example. 
    #As you can see it's quite tedious to write 
    #and becomes rather infeasible if there are more then, 
    #let's say 10, slope parameters to estimate 
    @setNLObjective( m, Max,-(d.n/2)*log(2π*σ^2) \\cont. next line
                            -sum{(d.y[i]-d.A[i,1]*b[1] \\
                                        -d.A[i,2]*b[2] \\
                                        -d.A[i,3]*b[3])^2, i=1:d.n}/(2σ^2) )

    #julia returns:
    > slope: [2.14,12.85,3.65], variance: 1.04

    #which is what is to be expected
    #however:

    #this is what I would like the code to look like:
    @setNLObjective( m, Max,-(d.n/2)*log(2π*σ^2) \\
                            -sum{(d.y[i]-(d.A[i,j]*b[j]))^2, \\
                            i=1:d.n, j=1:3}/(2σ^2) )

    #I also tried:
    @setNLObjective( m, Max,-(d.n/2)*log(2π*σ^2) \\
                            -sum{sum{(d.y[i]-(d.A[i,j]*b[j]))^2, \\
                            i=1:d.n}, j=1:3}/(2σ^2) )

    #but unfortunately it returns:
    > slope: [10.21,18.89,15.88], variance: 54.78

    solve(m)
    println( getValue(b), " ",  getValue(σ^2) )
end
ml()

Любые идеи?

ИЗМЕНИТЬ

Как отметил Реза, рабочий пример:

@setNLObjective( m, Max,-(d.n/2)*log(2π*σ^2) \\
                        -sum{(d.y[i]-sum{d.A[i,j]*b[j],j=1:3})^2,
                        i=1:d.n}/(2σ^2) )

julia julia-jump
person Vincent    schedule 17.10.2015    source источник

Ответы (2)


arrow_upward
1
arrow_downward

Я не отслеживал ваш код, но я хочу, чтобы у вас работало следующее:

sum([(d.y[i]-sum([d.A[i,j]*b[j] for j=1:3]))^2 for i=1:d.n])

как упоминал @IainDunning, пакет JuMP имеет специальный синтаксис для суммирования внутри макросов, поэтому более эффективный и абстрактный способ сделать это:

sum{sum{(d.y[i]-d.A[i,j]*b[j], j=1:3}^2,i=1:d.n}
person Reza Afzalan    schedule 17.10.2015
comment
Это было так. Большое спасибо. - person Vincent; 17.10.2015

arrow_upward
4
arrow_downward

Синтаксис sum{} - это специальный синтаксис, который работает только внутри макросов JuMP и является предпочтительным синтаксисом для сумм.

Итак, ваш пример будет записан как:

function ml(  )
    m = Model( solver = NLoptSolver( algorithm = :LD_LBFGS ) )
    @variable( m, b[1:3] )
    @variable( m, σ >= 0, start = 1.0 )

    @NLobjective(m, Max,
        -(d.n/2)*log(2π*σ^2)
        - sum{
            sum{(d.y[i]-d.A[i,j]*b[j], j=1:3}^2,
            i=1:d.n}/(2σ^2) )

где я расширил его на несколько строк, чтобы было как можно яснее.

Ответ Резы технически не ошибочен, но не является идиоматическим JuMP и не будет таким эффективным для более крупных моделей.

person IainDunning    schedule 17.10.2015
comment
Спасибо за ответ. Я знал, что это технически неправильно (из-за скобок, как вы правильно заметили, но он написал проблему именно так, как мне нужно. - person Vincent; 17.10.2015