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.
- 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.
sin(beta[4]x)
" shouldn't there be a*
betweenbeta[4]
andx
? $\endgroup$