Для шумных или тонко структурированных данных существуют ли лучшие квадратуры, чем правило средней точки?

Только первые два раздела этого длинного вопроса важны. Остальные только для иллюстрации.

Фон

Продвинутые квадратуры, такие как составные выражения Ньютона – Котса, Гауза – Лежандра и Ромберга более высокой степени, по-видимому, в основном предназначены для случаев, когда можно тонко выбрать функцию, но не интегрировать аналитически. Однако для функций со структурами, меньшими, чем интервал дискретизации (пример см. В Приложении A) или шум измерения, они не могут конкурировать с простыми подходами, такими как правило средней точки или трапеции (для демонстрации см. Приложение B).

Это несколько интуитивно понятно, поскольку, например, составное правило Симпсона по существу «отбрасывает» одну четверть информации, присваивая ей меньший вес. Единственная причина, по которой такие квадратуры лучше подходят для достаточно скучных функций, заключается в том, что правильная обработка граничных эффектов перевешивает эффект отбрасываемой информации. С другой точки зрения, мне интуитивно понятно, что для функций с тонкой структурой или шумом выборки, удаленные от границ области интегрирования, должны быть почти равноудаленными и иметь практически одинаковый вес (для большого числа выборок ). С другой стороны, квадратура таких функций может выиграть от лучшей обработки граничных эффектов (чем для метода средней точки).

Вопрос

Предположим, что я хочу численно интегрировать шумные или тонко структурированные одномерные данные.

Количество точек выборки является фиксированным (из-за того, что оценка функции обходится дорого), но я могу свободно разместить их. Однако я (или метод) не могу размещать точки выборки в интерактивном режиме, т. Е. На основе результатов из других точек выборки. Я также не знаю потенциальных проблемных регионов заранее. Итак, что-то вроде Гаусса-Лежандра (неэквидистантные точки отбора проб) в порядке; адаптивной квадратуры нет, так как она требует интерактивных точек выборки.

• Были ли предложены какие-либо методы, выходящие за пределы метода средней точки, для такого случая?

• Или: есть ли доказательства того, что метод средней точки лучше в таких условиях?

• В целом: есть ли какие-либо работы по этой проблеме?

Приложение A: Конкретный пример тонко структурированной функции

Я хотел бы оценить для: си. Типичная функция выглядит так:$\int_0^1f(t)\, \mathrm{d}t$

$$f(t) = \sum_{i=1}^{k} \frac{\sin(ω_i t-φ_i)}{ω_i},$$ log ω i ∈[1,1000$φ_i∈ [0,2π]$$\log{ω_i} ∈ [1,1000]$

Я выбрал эту функцию для следующих свойств:

• Он может быть интегрирован аналитически для контроля результата.
• Он имеет тонкую структуру на уровне, который делает невозможным захват всего этого с количеством образцов, которые я использую ( ).$<10^2$
• В нем не доминирует его тонкая структура.

Приложение B: контрольный показатель

Для полноты, вот пример в Python:

import numpy as np
from numpy.random import uniform
from scipy.integrate import simps, trapz, romb, fixed_quad

begin = 0
end   = 1

def generate_f(k,low_freq,high_freq):
    ω = 2**uniform(np.log2(low_freq),np.log2(high_freq),k)
    φ = uniform(0,2*np.pi,k)
    g = lambda t,ω,φ: np.sin(ω*t-φ)/ω
    G = lambda t,ω,φ: np.cos(ω*t-φ)/ω**2
    f = lambda t: sum( g(t,ω[i],φ[i]) for i in range(k) )
    control = sum( G(begin,ω[i],φ[i])-G(end,ω[i],φ[i]) for i in range(k) )
    return control,f

def midpoint(f,n):
    midpoints = np.linspace(begin,end,2*n+1)[1::2]
    assert len(midpoints)==n
    return np.mean(f(midpoints))*(n-1)

def evaluate(n,control,f):
    """
    returns the relative errors when integrating f with n evaluations
    for several numerical integration methods.
    """
    times = np.linspace(begin,end,n)
    values = f(times)
    results = [
            midpoint(f,n),
            trapz(values),
            simps(values),
            romb (values),
            fixed_quad(f,begin,end,n=n)[0]*(n-1),
        ]

    return [
            abs((result/(n-1)-control)/control)
            for result in results
        ]

method_names = ["midpoint","trapezoid","Simpson","Romberg","Gauß–Legendre"]

def med(data):
    medians = np.median(np.vstack(data),axis=0)
    for median,name in zip(medians,method_names):
        print(f"{median:.3e}   {name}")

print("superimposed sines")
med(evaluate(33,*generate_f(10,1,1000)) for _ in range(100000))

print("superimposed low-frequency sines (control)")
med(evaluate(33,*generate_f(10,0.5,1.5)) for _ in range(100000))

(Я здесь использую медиану, чтобы уменьшить влияние выбросов из-за функций, которые имеют только высокочастотное содержание. В среднем, результаты похожи.)

Медианы относительных ошибок интеграции:

superimposed sines
6.301e-04   midpoint
8.984e-04   trapezoid
1.158e-03   Simpson
1.537e-03   Romberg
1.862e-03   Gauß–Legendre

superimposed low-frequency sines (control)
2.790e-05   midpoint
5.933e-05   trapezoid
5.107e-09   Simpson
3.573e-16   Romberg
3.659e-16   Gauß–Legendre

Примечание: после двух месяцев и одного вознаграждения без результата я разместил это на MathOverflow .

Прежде всего, я думаю, что вы неправильно понимаете концепцию адаптивной квадратуры. Адаптивная квадратура не подразумевает «интерактивное размещение точек выборки». Вся идея адаптивной квадратуры состоит в том, чтобы разработать схему, которая будет интегрировать определенную функцию с определенной (оценочной) абсолютной или относительной погрешностью с как можно меньшими оценками функции.

Второе замечание: вы пишете: «Количество точек выборки фиксировано (из-за того, что оценка функции обходится дорого), но я могу свободно разместить их». Я думаю, что идея должна состоять в том, чтобы число точек отбора проб (или оценки функций в квадратурной терминологии) было как можно меньше (т.е. не фиксировано).

Так в чем же идея адаптивной квадратуры, реализованной, например, в QUADPACK ?

1. Основным ингредиентом является «вложенное» квадратурное правило: это сочетание двух квадратурных правил, где одно имеет более высокий порядок (или точность), чем другое. Почему? Основываясь на разнице между этими правилами, алгоритм может оценить квадратурную ошибку (конечно, алгоритм будет использовать наиболее точную в качестве контрольного результата). Примерами могут служить правило трапеции с узлами и узлами. В случае с QUADPACK правила являются правилами Гаусса-Кронрода. Это интерполяционные квадратурные правила, которые используют квадратурное правило Гаусса-Лежандра определенного порядка 2 n + 1 N N 2 N - $2^{n}$$2^{n+1}$$N$и оптимальное расширение этого правила. Это означает, что более высокий квадратурный порядок может быть получен повторным использованием узлов Гаусса-Лежандра (т. Е. Оценки дорогостоящих функций) с различными весами и добавлением ряда дополнительных узлов. Другими словами, оригинальное правило Гаусса-Лежандра порядка будет точно интегрировать все многочлены степени , тогда как расширенное правило Гаусса-Кронрода точно интегрирует некоторый многочлен более высокого порядка. Классическим правилом является G7K15 (Гаусс-Лежандр 7-го порядка с Гаусс-Кронродом 15-го порядка). Волшебство заключается в том, что 7 узлов Гаусса-Лежандра являются подмножеством 15 узлов Гаусса-Кронрода, поэтому с 15 оценками функций у меня есть квадратурная оценка вместе с оценкой ошибки!$N$$2N-1$

2. Следующий компонент - стратегия «разделяй и властвуй». Предположим, что вы выпустили этот G7K15 на свою интегральную функцию и наблюдаете квадратурную ошибку, которая по вашему вкусу слишком велика. Затем QUADPACK подразделяет исходный интервал на два равных интервала. И затем он пересмотрит эти два субинтеграла, используя основное правило, G7K15. Теперь алгоритм имеет глобальную оценку ошибки (которая, как мы надеемся, должна быть ниже первой), а также две локальные оценки ошибки. Он выбирает интервал с наибольшей ошибкой и делит его на два. Оцениваются два новых интеграла и обновляется глобальная ошибка. И так до тех пор, пока глобальная ошибка не окажется ниже запрошенной цели или не будет превышено максимальное количество подразделений.

Поэтому я призываю вас обновить код выше, используя scipy.quadметод. Возможно, в случае подынтегральной функции с большой «тонкой структурой» вам может потребоваться увеличить максимальное количество подразделений ( limitопция). Вы также можете поиграть с параметрами epsabsи / или epsrel.

Однако, если у вас есть только экспериментальные данные, я вижу две возможности.

1. Если у вас есть возможность выбрать точки измерения, то есть значения , я бы выбрал их на равном расстоянии, предпочтительно в виде степени чтобы вы могли применить вложенное правило трапеции (и получить прибыль от экстраполяции Ромберга).$t$$2$
2. Если у вас нет средств для выбора узлов, то есть измерения производятся в случайное время, лучшим вариантом, на мой взгляд, по-прежнему является правило трапеции.

Я не уверен, что ваш код демонстрирует что-то фундаментальное о различных квадратурных правилах и о том, насколько хорошо они справляются с шумом и тонкой структурой, и полагаю, что если вы выберете различную структуру штрафов, вы найдете что-то другое. Вот теорема:

Ни один квадратурный метод не может дать низкую абсолютную или относительную погрешность для функции с неограниченным полным изменением. В системе с плавающей запятой с единичным округлением мы имеем оценку$\mu$ где - квадратурная сумма, влияющая на численную реализацию функции .

$$\left| \int_{a}^{b} f \, \mathrm{d}x - \hat{Q}[\hat{f}] \right| \le \left| \int_{a}^{b} f \, \mathrm{d}x - Q[f] \right| + \mu\left[ 4\int_{a}^{b} |f| \, \mathrm{d}x + \int_{a}^{b} |xf'| \, \mathrm{d}x \right]$$ $\hat{Q}$$\hat{f}$$f$

Доказательство. Пусть квадратурные узлы будут а (неотрицательные) квадратурные веса - и обозначим их приближения с плавающей запятой через и . Предположим, что удовлетворяет где где - округление единицы. потом $\{x_i\}_{i=0}^{n-1}$$\{w_i\}_{i=0}^{n-1}$$\hat{w}_{i}$$\hat{x}_i$$\hat{f}$$\hat{f}(x) = f(x)(1+2\delta)$$|\delta| \le \mu$$\mu$

$$\begin{align*} \hat{Q}[\hat{f}] &= \sum_{i=0}^{n-1} \hat{w}_i \otimes \hat{f}(\hat{x}_i) \\ &= \sum_{i=0}^{n-1} w_i (1+\delta^w_i)f(x_i + \delta_i^x x_i)(1+2\delta_i^{f})(1+\delta_i^{*}) \\ &\approx \sum_{i=0}^{n-1} w_i \left[f(x_i)+ \delta_i^x x_i f'(x_i) \right] (1+\delta^w_i + 2\delta_i^{f} + \delta_i^*) \\ &\approx \sum_{i=0}^{n-1} w_i f(x_i) + \sum_{i=0}^{n-1}\delta_i^x w_i x_i f'(x_i) + w_i f(x_i) (\delta^w_i + 2\delta_i^{f} + \delta_i^*) \\ \end{align*}$$ так, чтобы Это предполагает, что сумма вычисляется без ошибок; умножить на чтобы отбросить это предположение. $$\begin{align*} |\hat{Q}[\hat{f}] - Q[f]| &\le \mu \sum_{i=0}^{n-1}w_i(|x_i f'(x_i)| + 4|f(x_i)|) \\ &\approx 4\mu \int |f| \, \mathrm{d}x + \mu \int |xf'| \, \mathrm{d}x \end{align*}$$ $n$

Mutatis mutandis вы также можете показать, что результат выполняется в арифметике с фиксированной запятой.