Как настроить метод трапециевидной интеграции на произвольную нулевую точку? - PullRequest
0 голосов
/ 11 января 2019

Я хочу вычислить интеграл от вектора последовательности. Поскольку функции недоступны, я использую трапециевидный метод 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() throws, потому что 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 должен проиллюстрировать принцип немного больше:

enter image description here


данные

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)
})

1 Ответ

0 голосов
/ 11 января 2019

Во-первых, есть еще одна небольшая, но более важная проблема, хотя после ее устранения ваш вопрос все равно останется в силе. Я имею в виду, что порядок x и y в качестве аргументов вашей функции должен быть обратным из-за того, как вы используете функцию в 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 

Бывает, что только два абсолютных значения отличаются от версии, в которой x и y обращены. Давайте посмотрим на первую функцию: она должна быть 9 или 15? У нас есть это 2 * 2/2 + 1 * 2 + 1 * 2/2 + 2 * 4/2 = 9, поэтому мы действительно хотим обратить x и y.

Другой способ написать функцию - это

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 (значения сетки) в порядке.

...