R: Integrar: Número máximo de subdivisiones alcanzadas, error de redondeo

Question

Me encontré con un problema interesante pero bastante molesto.

Estoy tratando de integrar una función que se ha calculado a partir de un conjunto de datos. Los datos se pueden encontrar aquí: Link to sample.txt.

Empiezo por ajustar una línea a mis datos. esto se puede hacer de forma lineal con approxfun o no lineal con splinefun. En mi ejemplo a continuación, uso este último. Ahora, cuando intento de integrar la función ajustada corro en el error

• maximum number of subdivisions reached

pero cuando puedo aumentar la subdivisión consigo

• roundoff error

De los valores en mi código de muestra, puede ver eso para esta información específica establecer el umbral es 754-> 755.

Mi colega no tiene problemas para integrar este conjunto de datos en Matlab. ¿Hay alguna manera de manipular mis datos para integrarlos? ¿Hay otro método para la integración numérica en R?

data<-read.table('sample.txt',sep=',') 
colnames(data)<-c('wave','trans') 
plot(data$wave,data$trans,type='l') 

trans<- -1 * log(data$trans) 
plot(data$wave,trans,type='l') 

fx.spline<-splinefun(data$wave,trans) 

#Try either 
Fx.spline<-integrate(fx.spline,min(data$wave),max(data$wave)) 
#Above: Number of subdivision reached 
Fx.spline<-integrate(fx.spline,min(data$wave),max(data$wave),subdivisions=754) 
#Above: Number of subdivision reached 
Fx.spline<-integrate(fx.spline,min(data$wave),max(data$wave),subdivisions=755) 
#Above: Roundoff error 

Answer 1

Hay muchas rutinas de integración en I, y se pueden encontrar algunos de ellos por 'RSiteSearch'ing o utilizando el 'paquete SOS'.

Por ejemplo, el paquete de pracma tiene varias implementaciones, por ejemplo

quad(fx.spline,min(data$wave),max(data$wave)) # adaptive Simpson 
# [1] 2.170449         # 2.5 sec 
quadgk(fx.spline,min(data$wave),max(data$wave)) # adaptive Gauss-Kronrod 
# [1] 2.170449         # 0.9 sec 
quadl(fx.spline,min(data$wave),max(data$wave)) # adaptive Lobatto 
# [1] 2.170449         # 0.8 sec 

no Por favor, que estos son secuencias de comandos puros R y por lo tanto más lento que, por ejemplo, el compilado integrate rutina con una función de este tipo oscilante.