1

R の積分関数を使用していますが、被積分関数が問題を引き起こしています (特異点があります)。このことを考慮、

distrib <- function(x){
  dnorm(x, mean=700,sd=50)
}

phi <- function(nu, epsilon=20, maxi=Inf){

  fun <- function(nup)  {
    lambdap <- 1e7/nup
    distrib(lambdap) / (nup*(nup - nu))
  }
  # first part: from epsilon to nu - epsilon
  I1 <- integrate(fun, epsilon, nu-epsilon, rel.tol = 1e-8,
                  subdivisions = 200L)
  # second part: from nu + epsilon to Infty
  I2 <- integrate(fun, nu+epsilon, maxi, rel.tol = 1e-8, 
                  subdivisions = 200L)

  I1$value + I2$value 

}


x <- seq(1e7/500, 1e7/1000, length=200)

phi1 <- sapply(x, phi)    
phi1 <- sapply(x, phi, epsilon=5)
#Error in integrate(fun, epsilon, nu - epsilon, rel.tol = 1e-08, subdivisions = #200L) : 
#  the integral is probably divergent

カットオフパラメータをいじって結果が得られるまで手動で調整するよりも、Rでそのような主要な値を計算するより良い方法はありますか(必ずしも正確ではありません)

4

1 に答える 1

1

特異点を中心とする区間の積分を検討する場合、変数変換を使用して被積分関数を対称化できます。一般化されたコーシーの主値の数値評価（A.NyiriおよびL.Bayanyi）で詳しく説明されているように、特異点は消えます。 、1999）。

# Principal value of \( \int_{x_0-a}^{x_0+a} \dfrac{ f(x) }{ h(x) - h(x0) } dx \)
pv_integrate <- function( f, h, x0, a, epsilon = 1e-6 ) {
  # Estimate the derivatives of f and h
  f1 <- ( f(x0+epsilon) - f(x0-epsilon) ) / ( 2 * epsilon )
  h1 <- ( h(x0+epsilon) - h(x0-epsilon) ) / ( 2 * epsilon )
  h2 <- ( h(x0+epsilon) + h(x0-epsilon) - 2 * h(x0) ) / epsilon^2 
  # Function to integrate. 
  # We just "symmetrize" the integrand, to get rid of the singularity,
  # but, to be able to integrate it numerically, 
  # we have to provide a value at the singularity.
  # Details: http://mat76.mat.uni-miskolc.hu/~mnotes/downloader.php?article=mmn_16.pdf
  g <- function(u) 
    ifelse( u != 0,
      f( x0 - u ) / ( h(x0 - u) - h(x0) ) + f( x0 + u ) / ( h(x0+u) - h(x0) ),
      2 * f1 / h1 - f(x0) * h2 / h1^2
    )
  integrate( g, 0, a )
}

# Compare with the numeric results in the article
pv_integrate( function(x) 1,      function(x) x,    0, 1  ) # 0
pv_integrate( function(x) 1,      function(x) x^3,  1, .5 ) # -0.3425633
pv_integrate( function(x) exp(x), function(x) x,    0, 1  ) # 2.114502
pv_integrate( function(x) x^2,    function(x) x^4,  1, .5 ) # 0.1318667
于 2013-03-07T17:48:09.893 に答える