【问题标题】:How to adjust trapezoidal integration method to a custom zero point?如何将梯形积分方法调整为自定义零点?
【发布时间】:2019-06-06 14:27:19
【问题描述】:

我想计算一个序列向量的积分。由于没有可用的函数,我使用trapezoidal方法1

iglTzm <- function(x, y) sum(diff(x) * (head(y, -1) + tail(y, -1))) / 2

序列的第一个元素应该是零点,所以原则是:如果序列的值主要低于第一个值,则积分应该为负,否则为正,或为0。

考虑矩阵m1

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    6    7    8    8    6    8   10
[2,]    9    9    8    9    9    8    9
[3,]    9   10   10    9    9    9    9
[4,]    9    8    8    8    6    8    9
[5,]   10   10   10    9   10    8    0
[6,]    9    8    9   10    9    9    9

与这些原始值集成很可能会导致值不一致:

> setNames(apply(m1, 1, iglTzm, 0:6), 1:6)
  1   2   3   4   5   6 
 15   2  -2   7 -52   0 

所以我根据第一个值(第 1 列)调整序列(行),以便设置正确的符号,并获得矩阵 m2

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    0    1    2    2    0    2    4
[2,]    0    0   -1    0    0   -1    0
[3,]    0    1    1    0    0    0    0
[4,]    0   -1   -1   -1   -3   -1    0
[5,]    0    0    0   -1    0   -2  -10
[6,]    0   -1    0    1    0    0    0

从逻辑上讲,这不会改变 iglTzm() 抛出的任何值,因为 diff() 是相同的:

> setNames(apply(m2, 1, iglTzm, 0:6), 1:6)
  1   2   3   4   5   6 
 15   2  -2   7 -52   0 

无论如何,因为我不能简单地缩放或反转它,所以我还没有一个绝妙的主意如何调整函数以获得正确的符号,假设是:

#  1   2   3   4   5   6 
# 15  -2   2  -7 -52   0

有谁知道如何调整iglTzm() 以获得正确符号的积分?

m2 的情节应该更能说明原理:


数据

m1 <- matrix(c(6, 7, 8, 8, 6, 8, 10,
                9, 9, 8, 9, 9, 8, 9,
                9, 10, 10, 9, 9, 9, 9,
                9, 8, 8, 8, 6, 8, 9, 
                10, 10, 10, 9, 10, 8, 0, 
                9, 8, 9, 10, 9, 9, 9), 6, byrow=TRUE)

m2 <- t(apply(m1, 1, function(x) scale(x, center=x[1], scale=FALSE)))

# plot
par(mfrow=c(2, 3))
lapply(1:nrow(m2), function(x) {
  plot(m2[x, ], type="l", main=x)
  abline(h=m2[x, 1], col="red", lty=2)
})

【问题讨论】:

    标签: r numerical-integration


    【解决方案1】:

    首先,还有一个小但更重要的问题,尽管在修复它之后你的问题仍然有效。我的意思是,xy 作为函数参数的顺序应该颠倒,因为你在 apply 中使用函数的方式。

    但这还不够,现在我们回到您的问题。为此,让我们回顾一下通常的积分:ʃf(x)dx(限制从 a 到 b)将积分 f 下方的区域,这是您的函数已经成功完成的。现在你想要的是调整它的水平。但是如果我们从a积分到b,则等于ʃ(f(x)-f(a))dx = ʃf(x)dx - (b-a)f(a),从而导致

    iglTzm <- function(y, x) sum(diff(x) * (head(y, -1) + tail(y, -1))) / 2 - y[1] * (max(x) - min(x))
    setNames(apply(m1, 1, iglTzm, 0:6), 1:6)
    #  1  2  3  4  5  6 
    #  9 -2  2 -7 -8  0 
    

    碰巧只有两个绝对值与 xy 颠倒的版本不同。让我们看看第一个函数:它应该是 9 还是 15?我们有 2*2/2 + 1*2 + 1*2/2 + 2*4/2 = 9,所以我们确实想要反转 xy

    另一种编写函数的方法是

    iglTzm <- function(y, x) sum(diff(x) * (head(y - y[1], -1) + tail(y - y[1], -1))) / 2
    setNames(apply(m1, 1, iglTzm, 0:6), 1:6)
    #  1  2  3  4  5  6 
    #  9 -2  2 -7 -8  0 
    

    编辑:我说的颠倒只是指函数定义中的顺序,或者你在apply 中如何使用它;就 y(函数值)和 x(网格值)而言,函数本身很好。

    【讨论】:

    • 你好!我注意到了奇怪的值,但我还没有真正意识到错误。但我相信现在的顺序更好,可以直观地与apply或其他类似架构的函数一起使用。我确定发生混淆是因为我之前使用过integrate()。现在 - 不管你已经很好的答案 - 在这一点上我有一个直觉问题:integrate(function(x) x, 0, 1)iglTzm(0:1, 0:1)yield 相同,而 iglTzm(0:10, 0:1) 不同,但不应该相等吗?我们在所有情况下从 0 到 1 进行积分,并且所有斜率都相等。希望这是有道理的。
    • @jay.sf, iglTzm(0:10, 0:1) 是错误的(由于回收,您不会收到警告/错误)。你有 11 个函数值 0、1、...、10,而你只给出了两个网格点 0、1。那么函数在哪里取值 3、6、10?现在iglTzm(0:10, seq(0, 1, length = 11)) 将是从 0 到 1 的正确积分,但斜率不一样。 0:1的斜率是1,而0:10的斜率是10。(由于回收,iglTzm(0:10, 0:1)iglTzm(0:10, 0:10)相同。)
    • 啊回收!好的,但是当我 plot(0:1)plot(0:10) 时,两个斜率应该相同,并且从 0 到 1 的两个积分也相同,不是吗?也许我缺少什么。
    • @jay.sf, plot(0:1) 对应于 f(x)=x 对于x 的范围为[0,1],而plot(0:10) 与范围[0] 相同的函数,10]。那么第一个积分是 1*1/2=1,而第二个积分是 10*10/2=50(它们只是两个三角形)。 xy 之间似乎有些混淆,你传递给 iglTzmxy 就 f(x)=y 作为实际函数而言,并且可能是什么斜率(例如,通过重新缩放绘图,斜率可能看起来不同,但事实并非如此)。
    • 是的,我认为你的话在我的脑海中打了个结。非常感谢您的帮助,谢谢!
    猜你喜欢
    • 1970-01-01
    • 2020-05-23
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2019-07-31
    • 1970-01-01
    • 2021-06-02
    • 1970-01-01
    相关资源
    最近更新 更多