R中的Lebesgue-Stieltjes积分

Question

我用数值计算了分布函数。First I绘制函数。似乎在0.05左右就错了。请问这是因为四舍五入的错误吗？

​

第二，我需要找到相应的第一和第二非中心时刻。那是,

X= int (X)

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

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

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

X= 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的值会在您表示关注的区域产生类似的混乱结果。

​

​

所以我想也许你的功能表现不够好，不能产生好的效果。