Non-Linear Least Squares Sine Frequency Estimation in julia

Question

I am trying to estimate an artificial non-linear maximum likelihood problem which works fine as long as there is no frequency in a sine function to estimate.

I am using julia Optim and Distributions for this example, here is the code:

function ml(  )
    m = 10000
    x = linspace(0,1,m)

    #non linear function
    y = 2 + 3.4*sin(x) + exp(1.1*x) + rand(Normal(0,1), m)

    #use maximum likelihood of the normally distributed error terms
    function lglklhd( beta  )
        err = (y - beta[2] - beta[3]*sin(x) - exp(beta[4]*x)).^2
        -loglikelihood( Normal(0, beta[1]), err )
    end
    res = optimize( lglklhd, [0.1, 0.0, 0.0, 0.0] )
    println( res.minimum )
    println( "1, 2, 3.4, 1.1" )
end

which returns:

> [1.72, 2.00, 3.31, 1.12]

The first number is the estimation of the variance, followed by the estimation of the parameters. Note that the error in the variance estimation is relatively larger then the error of the other parameters.

1. Question: Why? It should be close to 1. right?

Now let me change:

#small change in the response function
y = 2 + 3.4*sin(2.1*x) + exp(1.1*x) + rand(Normal(0,1), m)

#and here:
err = (y - beta[2] - beta[3]*sin(beta[4]x) - exp(beta[5]*x)).^2

returns:

> [1.71,1.99,-4.79, -1.711, 0.14]
1, 2, 3.4, 2.1, 1.1 #original parameter

These are quite different from the ones without the frequency in the sine function.

First: both sine parameters have the wrong sign and are simple relatively far from their true value.

Second: the estimator of the exponential is now notably lower then it was before.

Overall I must conclude the additional parameter in the sine function distorts the estimation heavily.

Second question: Why?

I also tried other optimization algorithms but the results were almost always worse, if I did get any.

EDIT

In response to suggestions made here I used JuMP and the NLopt solver in an attempt to solve the problem. First from the NLopt Documentation I used the following Local derivative-free optimization algorithms:

Code:

...
#ceteris paribus (almost, variance estimation being left out.) 
...

#using a local (L) derivative free (N) COBYLA algorithm
model = Model( solver = NLoptSolver( algorithm = :LN_COBYLA ) )
@defVar( model, beta[1:4] )
@setNLObjective( model, Min, sum{ ( y[i] - beta[1] - beta[2] * sin(beta[3] * x[i]) - exp(beta[4] * x[i]) )^2, i=1:m  } )
solve( model )
println( getValue(beta) )

which results in:

> [2.00,-3.36, -2.06, 1.12]

While the non-sine related parameters are estimated quite alright, the sine related ones have the wrong sign and value.

Then this:

...
model = Model( solver = NLoptSolver( algorithm = :LN_BOBYQA ) ) 
...

> [2.07, -4.18, -1.77, 0.56]

...
model = Model( solver = NLoptSolver( algorithm = :LN_NELDERMEAD ) ) 
...

> [2.01, 3.26, 2.20, 1.16]
> [1.94, 5.54, 1.73, -0.99]

which is far better. So I assume there must be something with the Nelder-Mead implementation in the Optim Package.

I assume this is not a proper mathematical "way" of getting workable estimates but running this algorithm a few thousand times and taking the average should provide estimates that can be worked with.

Answer 1

The reason why you have trouble finding optimal values is that the sine function is NOT convex which is one of the basic requirement for optimization algorithms that is based on derivatives. Non convex function have multiple minima or multiple maxima so it is the impossible to find global optima. If you try derivative free optimizes like genetic algorithms or simulated annealing you will find close to global optimum and if you run multiple times you will be able to find almost all global optimum.

You could plot your function and see that it it has multiple optimum points.