如何在 R 中进行向量积分？ [复制]答案

【问题标题】：how can I do vector integral in R? [duplicate]如何在 R 中进行向量积分？ [复制]
【发布时间】：2015-09-07 06:14:41
【问题描述】：

我想在 R 中集成一个一维向量，我应该怎么做？假设我有：

d=hist(p, breaks=100, plot=FALSE)$density

其中 p 是一个样本，例如：

p=rnorm(1e5)

如何计算 d 上的积分？

【问题讨论】：

标签： r vector integral

【解决方案1】：

如果我们假设d 中的值对应于函数的 y 值，那么我们可以使用离散近似来计算积分。例如，我们可以为此使用梯形规则或辛普森规则。然后，我们还需要输入与 x 轴上的离散区间相对应的步长，以便“近似曲线下的面积”。

离散积分函数定义如下：

p=rnorm(1e5)
d=hist(p,breaks=100,plot=FALSE)$density

discreteIntegrationTrapeziumRule <- function(v,lower=1,upper=length(v),stepsize=1)
{
  if(upper > length(v))
    upper=length(v)
  if(lower < 1)
    lower=1

  integrand <- v[lower:upper]
  l <- length(integrand)
  stepsize*(0.5*integrand[1]+sum(integrand[2:(l-1)])+0.5*v[l])
}

discreteIntegrationSimpsonRule <- function(v,lower=1,upper=length(v),stepsize=1)
{
  if(upper > length(v))
    upper=length(v)
  if(lower < 1)
    lower=1

  integrand <- v[lower:upper]
  l <- length(integrand)
  a = seq(from=2,to=l-1,by=2);
  b = seq(from=3,to=l-1,by=2)
  (stepsize/3)*(integrand[1]+4*sum(integrand[a])+2*sum(integrand[b])+integrand[l])
}

例如，让我们在假设离散 x 步长为 1 的情况下近似曲线下的整个区域，然后对 d 的后半部分做同样的事情，同时我们假设 x 步长为 0.2。

> plot(1:length(d),d) # stepsize one on x-axis
> resultTrapeziumRule <- discreteIntegrationTrapeziumRule(d) # integrate over complete interval, assume x-stepsize = 1
> resultSimpsonRule <- discreteIntegrationSimpsonRule(d) # integrate over complete interval, assume x-stepsize = 1
> resultTrapeziumRule
[1] 9.9999
> resultSimpsonRule
[1] 10.00247

> plot(seq(from=-10,to=(-10+(length(d)*0.2)-0.2),by=0.2),d) # stepsize 0.2 on x-axis
> resultTrapziumRule <- discreteIntegrationTrapeziumRule(d,ceiling(length(d)/2),length(d),0.2) # integrate over second part of vector, x-stepsize=0.2
> resultSimpsonRule <- discreteIntegrationSimpsonRule(d,ceiling(length(d)/2),length(d),0.2) # integrate over second part of vector, x-stepsize=0.2
> resultTrapziumRule
[1] 1.15478
> resultSimpsonRule
[1] 1.11678

一般来说，辛普森规则提供了更好的积分近似值。您拥有的 y 值越多（x 轴步长越小），您的近似值就会越好。

小EDIT为清楚起见：在这种特殊情况下，步长显然应该是 0.1。正如预期的那样，密度曲线下的完整面积（大约）等于 1。

> d=hist(p,breaks=100,plot=FALSE)$density
> hist(p,breaks=100,plot=FALSE)$mids # stepsize = 0.1
 [1] -4.75 -4.65 -4.55 -4.45 -4.35 -4.25 -4.15 -4.05 -3.95 -3.85 -3.75 -3.65 -3.55 -3.45 -3.35 -3.25 -3.15 -3.05 -2.95 -2.85 -2.75 -2.65 -2.55
[24] -2.45 -2.35 -2.25 -2.15 -2.05 -1.95 -1.85 -1.75 -1.65 -1.55 -1.45 -1.35 -1.25 -1.15 -1.05 -0.95 -0.85 -0.75 -0.65 -0.55 -0.45 -0.35 -0.25
[47] -0.15 -0.05  0.05  0.15  0.25  0.35  0.45  0.55  0.65  0.75  0.85  0.95  1.05  1.15  1.25  1.35  1.45  1.55  1.65  1.75  1.85  1.95  2.05
[70]  2.15  2.25  2.35  2.45  2.55  2.65  2.75  2.85  2.95  3.05  3.15  3.25  3.35  3.45  3.55  3.65  3.75  3.85  3.95  4.05  4.15

> resultTrapeziumRule <- discreteIntegrationTrapeziumRule(d,stepsize=0.1) 
> resultTrapeziumRule
[1] 0.999985

【讨论】：