I'm trying to da an optimization for the followig function
Zwischenwert <- function (x)
{
lambda<-x[1];
mu<-x[2];
gammal<-x[3];
mud<-x[4];
gammad<-x[5];
Mittelwert <-0;
for(t in 0:(T-1))
{
for(i in 0:(n-1))
{
for(j in i:(n-1))
{
Mittelwert <- Mittelwert +(
phi[i+1,j+1,t+1]*((j-i)*log(lambda)-log(factorial(j-i))-lambda));
};
if(t>0){Mittelwert <- Mittelwert
+(phisum[i+1,t+1]*(-0.5*log(2*pi)-0.5*log(t*gammal+i*gammad)-((Y[t+1]-i*mud-t*mu-0.5*t*gammal+lambda*exp(mud+0.5*gammad)-t*lambda)^2
)/(2*(t*gammal+i*gammad))));};
};
};
};
and it's gradient
ZWGrad <- function(x)
{
lambda<-x[1];
mu<-x[2];
gammal<-x[3];
mud<-x[4];
gammad<-x[5];
dlambda<-0;
dmu<-0;
dgamma<-0;
dmud<-0;
dgammad<-0;
for (t in 0:(T-1))
{
for (i in 0:(n-1))
{
for (j in i:(n-1))
{
dlambda<-dlambda-phi[i+1,j+1,t+1]*(-1+(j-i)/lambda);
};
};
};
for (t in 1:(T-1))
{
for (i in 0:(n-1))
{
dlambda<-dlambda+phisum[i+1,t+1]*(t-t*lambda*exp(mud+0.5*gammad))*(Y[t+1]-i*mud-mu*t-0.5*gammal*t+lambda*t*exp(mud+0.5*gammad)-lambda*t)/(2*(gammal*t+gammad*i));
dmu<-dmu+phisum[i+1,t+1]*t*(Y[t+1]-i*mud-t*mu-0.5*gammal*t+lambda*t*exp(mud+0.5*gammad)-lambda*t)/(gammal*t+gammad*i);
dgamma<-dgamma+phisum[i+1,t+1]*(
t*(0.5*i*gammad+Y[t+1]-i*mud-t*mu+lambda*t*exp(mud+0.5*gammad)-lambda*t)^2/(2*(gammal*t+gammad*i)^2)-t/(2*(gammal*t+gammad*i))-0.125*t);
dmud<-dmud+phisum[i+1,t+1]*(i-lambda*t*exp(mud+0.5*gammad))*(Y[t+1]-i*mud-t*mu-0.5*gammal*t+lambda*t*exp(mud+0.5*gammad)-lambda*t)/(gammal*t+gammad*i);
dgammad<-dgammad+phisum[i+1,t+1]*(Y[t+1]-i*mud-t*mu-0.5*gammal*t+lambda*t*exp(mud+0.5*gammad)-lambda*t)*(i*(Y[t+1]-i*mud-t*mu-0.5*gammal*t-lambda*t)-lambda*t*exp(mud+0.5*gammad)*(gammal*t+gammad*i-i))/(2*(gammal*t+gammad*i)^2);
};
};
GradZW <- c(dlambda,dmu,dgamma,dmud,dgammad);
};
lambda, gammal und gammd need to be larger than 0, so I've tried different
constrained methods like L-BFGS-B via optim and constrOptim, but both seem to
have problems with my functions. The standard optim-method works, but sometimes
it returns negative values for lambda, gammal and gammad.
Those phi and phisums are known at runtime, n and T can be pretty large, so
maybe R simple can't handel it.
Many thanks
Matthias
--