Hi,
I have a density that I need to get MLEs from, which includes definite
integrals both in the denominator and in the numerator of the density function.
It looks like the outcome depends on the initial values given. My program is
shown below:
library(circular)
########################################
4 parameters
########################################
z<-rvonmises(100,0,1)
z<-as.vector(z)
x<-rvonmises(100,0,2+cos(z))
x<-as.vector(x)
f<-function(x){
k1<-x[1]
k2<-x[2]
k12<-x[3]
g<-x[4]
-2*sum(log(
exp(k1*cos(x))*integrate(function(y)exp((k2+k12*cos(x))*cos(y)),g,2*pi,subdivisions=1000000)$value
/(2*pi*integrate(function(y)besselI(k1+k12*cos(y),0)*exp(k2*cos(y)),g,2*pi,subdivisions=1000000)$value)
))
}
gr<-function(u){
k1<-u[1]
k2<-u[2]
k12<-u[3]
g<-u[4]
grad1<-sum(integrate(function(y)integrate(function(a)(cos(x)-cos(a))*exp(k2*cos(y)+(k1+k12*cos(y))*cos(a)),subdivisions=10000000,-pi,pi)$value,subdivisions=10000000,g,pi)$value)
grad2<-sum(integral(function(y)cos(y)*exp((k2+k12*cos(x))*cos(y)),subdivisions=10000000,g,pi)$value*integrate(function(y)besselI(k1+k12*cos(y),0)*exp(k2*cos(y)),subdivisions=10000000,g,pi)$value-integrate(function(y)exp((k2+k12*cos(x))*cos(y)),subdivisions=10000000,g,pi)$value*integrate(function(y)besselI(k1+k12*cos(y),1)*cos(y)*exp(k2*cos(y))+besselI(k1+k12*cos(y),0)*cos(y)*exp(k2*cos(y)),subdivisions=10000000,g,pi)$value)
grad3<-sum(integrate(function(y)cos(x)*cos(y)*exp((k2+k12*cos(x))*cos(y)),subdivisions=10000000,g,pi)$value*integrate(function(y)besselI(k1+k12*cos(y),0)*exp(k2*cos(y)),subdivisions=10000000,g,pi)$value-integrate(function(y)exp((k2+k12*cos(x))*cos(y)),subdivisions=10000000,g,pi)$value*integrate(function(y)besselI(k1+k12*cos(y),1)*cos(y)*exp(k2*cos(y)),subdivisions=10000000,g,pi)$value)
grad4<-sum(integrate(function(y)besselI(k1+k12*cos(y),0)*exp(k2*cos(y)),subdivisions=10000000,g,pi)$value*exp((k2+k12*cos(x))*cos(g))-integrate(function(y)exp((k2+k12*cos(x))*cos(y)),subdivisions=10000000,g,pi)$value*besselI(k1+k12*cos(g))*exp(k2*cos(g)))
c(grad1,grad2,grad3,grad4)
}
constrOptim(c(2,2,12,2),f,gr,method="CG",ui=rbind(c(1,0,0,0),c(0,1,0,0),c(0,0,1,0),c(0,0,0,1),c(0,0,0,-1)),ci=c(0,0,0,-pi,-pi),control=list(fnscale=-1))
I get 3 different error messages for different initial values given. They are
shown below:
Error: cannot allocate vector of size 305.2 Mb
or
Error in integrate(function(y) exp((k2 + k12 * cos(x - m1)) * cos(y - :
evaluation of function gave a result of wrong length
or
Error in integrate(function(y) exp((k2 + k12 * cos(x)) * cos(y)), g, 2 * :
the integral is probably divergent
I would appreciate your help.
Sincerely,
---------------------------------
[[alternative HTML version deleted]]
