7

私は興味深いがかなり厄介な問題に出くわしました。

データセットから計算された関数を統合しようとしています。データはここにあります：sample.txtへのリンク。

まず、データに線を当てます。approxfunこれは、で線形または非線形で行うことができますsplinefun。以下の例では、後者を使用しています。さて、フィットした関数を統合しようとすると、エラーが発生します

  • maximum number of subdivisions reached

しかし、細分化を増やすと、

  • roundoff error

サンプルコードの値から、この特定のデータセットのしきい値は754->755であることがわかります。

私の同僚は、このデータセットをMatlabに統合するのに問題はありません。データを操作して統合する方法はありますか？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
4

1 に答える 1

7

Rには多くの統合ルーチンがあり、それらのいくつかは「RSiteSearch」または「sos」パッケージを使用して見つけることができます。

たとえば、パッケージpracmaにはいくつかの実装があります。

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

これらは純粋なRスクリプトであるため、たとえば、このintegrateような振動関数を使用してコンパイルされたルーチンよりも低速であることに注意してください。

于 2012-05-04T19:04:16.210 に答える