CODEWITHSUNDEEP
juliajulia-jump
Vincent
Vincent

Reputation: 1371

Julia JuMP Multiavariate ML Estimation

I am trying to perform a ML-Estimation of a normally distributed variable in a linear regression setting in Julia using JuMP and the NLopt solver.

There exists a good working example here however if I try to estimate the regression parameters (slope) the code becomes quite tedious to write, in particular if the parameter space increases.

Maybe someone has an idea how to write it more concise. Here is my Code:

#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()

Any ideas?

EDIT

As noted by Reza a working example is:

@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) )

Upvotes: 3

Views: 270

Answers (2)

IainDunning
IainDunning

Reputation: 11644

The sum{} syntax is a special syntax that only works inside JuMP macros, and is the preferred syntax for sums.

So your example would be written as:

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) )

where I've expanded it across multiple lines to be as clear as possible.

Reza's answer isn't technically wrong, but isn't idiomatic JuMP and won't be as efficient for larger models.

Upvotes: 4

Reza Afzalan
Reza Afzalan

Reputation: 5746

I didn't trace your code but anywhere, I wish that the following works for you:

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

as @IainDunning mentioned, JuMP package has a special syntax for summation inside it's macros, so the more efficient and abstract way to do this is:

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

Upvotes: 1

Related Questions