r - R 中的 Lebesgue-Stieltjes 集成

Question

我用数值方法计算了分布函数。 首先我绘制函数。看起来错误在 0.05 左右。请问这是由于舍入错误吗？

第二，我需要找到对应的第一和第二非中心矩。也就是说，

EX = int x dF(x)

EX^2 = int x^2 dF(x)

我可以在 R 中进行这种类型的 Lebesgue-Stieltjes 集成吗？请问有内置方法吗？如果R中没有，请问什么包提供了这样的计算？

我想或者，我可以找到 F(x) 的数值微分 f(x)，然后进行通常的积分，如

EX = int x f(x) dx

但我记得在某个地方数值微分不太稳定。请问哪种方法是正确的？

仅供引用，我的功能附在下面。

library(mvtnorm)
library(matrixcalc) 

VAR <- matrix(c(1.043856e-03, 5.044899e-04, 3.239951e-04, 2.330992e-04, 0.0001779055, 0.0001403866, 0.0001127118, 9.074962e-05, 7.157144e-05,
                5.044899e-04, 5.485889e-04, 3.523165e-04, 2.534751e-04, 0.0001934568, 0.0001526582, 0.0001225642, 9.868232e-05, 7.782773e-05,
                3.239951e-04, 3.523165e-04, 3.878844e-04, 2.790645e-04, 0.0002129870, 0.0001680697, 0.0001349376, 1.086447e-04, 8.568475e-05,
                2.330992e-04, 2.534751e-04, 2.790645e-04, 3.123147e-04, 0.0002383642, 0.0001880950, 0.0001510153, 1.215896e-04, 9.589399e-05,
                1.779055e-04, 1.934568e-04, 2.129870e-04, 2.383642e-04, 0.0002728857, 0.0002153361, 0.0001728863, 1.391990e-04, 1.097820e-04,
                1.403866e-04, 1.526582e-04, 1.680697e-04, 1.880950e-04, 0.0002153361, 0.0002548851, 0.0002046389, 1.647645e-04, 1.299447e-04,
                1.127118e-04, 1.225642e-04, 1.349376e-04, 1.510153e-04, 0.0001728863, 0.0002046389, 0.0002555744, 2.057751e-04, 1.622886e-04,
                9.074962e-05, 9.868232e-05, 1.086447e-04, 1.215896e-04, 0.0001391990, 0.0001647645, 0.0002057751, 2.840218e-04, 2.239993e-04,
                7.157144e-05, 7.782773e-05, 8.568475e-05, 9.589399e-05, 0.0001097820, 0.0001299447, 0.0001622886, 2.239993e-04, 3.974881e-04),
              nrow=9, ncol=9, byrow=TRUE)

is.symmetric.matrix(VAR)
is.positive.definite(VAR)
kappa(VAR)

CDF <- function(x){
  summand <- rep(0, 5)

  for(j in 5:9){
    choice <- combn(9, j)
    for(i in 1:ncol(choice)){
      ub <- rep(Inf, 9)
      ub[choice[, i]] <- x
      summand[j-4] <- summand[j-4] + as.numeric(pmvnorm(lower=rep(-Inf, 9), upper=ub, sigma=VAR))
    }
  }

  l <- c(1, -5, 15, -35, 70)
  as.numeric(t(l)%*%summand)
}

CDF <- Vectorize(CDF)

x <- seq(-0.1, 0.1, by=0.01)
y <- CDF(x)

plot(x, y, type="l", lwd=2)

Answer 1

我最初绘制了从 numCDF <- CDF( seq(-10, 10, length=100) ) 中获取一阶差异得到的结果，但这相当令人失望，因为只有一个值不同于 0。所以我将焦点限制在:

 numCDF <- CDF( seq(-.10, .10, length=100) )
 plot( diff(numCDF) )

简单地绘制 numCDF 的值会在您表示关注的区域中产生类似的困惑结果。

所以我认为你的函数可能表现得不够好，无法产生良好的结果。