Comments 23
jax.grad с allow_int = True не пробовали?
Спасибо за комментарий. Обязательно попробую. Я, собственно, собирался про МКЭ сразу писать с переходом на 2- и 3-мерные задачи. Там с визуализации и градиентные методы и т.д. Только неожиданно застрял в Numpy. Оказывается, к opensource радо аккуратно относиться. Причем всегда :))
Попробовал. Если задать функцию аналитически, выдает производную в точке или на массиве точек. Это не то, что нужно. Как ни пытался "запихнуть" массив значений в аргументы jax.grad ( в том числе как return функции целочисленного аргумента), не получилось. Возможно, не догоняю.
В NumPy (насколько я понимаю) для этого есть функция vectorize. В Jax, соответственно jax.vmap и jax.numpy.vectorize.
Проблема, однако, в том, что для получения такой "векторизированной" функции, принимающей на вход дискретный массив значений, надо этой самой vectorize скормить обычную "формульную" функцию. Получается замкнутый круг.
Спасибо за комментарий. Надо сказать, я сильно удивлен интересом к моей маленькой заметке. Сначала я думал, что ограничусь одной публикацией на тему дифференцирования, и то только потому, что не все оказалось в порядке в библиотеках Питона. Причем именно с дискретными функциями, например, полученными в эксперименте (физическом или математическом). Потом планировал сразу перейти к диффурам, потому что с ОДУ 2-го порядка также все непросто в NumPy и SciPy. От них уже можно перейти ближе к ML, в частности, к обсуждению методов выбора функций ошибок в плане сходимости градиентых методов. На тот счет также есть идеи. Однако все пошло по другому пути, особенно, когда я прочитал публикацию https://habr.com/ru/articles/749288/. В результате МКЭ переместился в планируемую вторую часть публикации и туда же добавилась попытка приблизиться к реальным задачам. Как всегда бывает, теория не всегда удовлетворительно работает на практике. Приходится додумывать.
Недельку посачковал, думаю на следующей неделе выгрузить продолжение (боюсь, не последнее). А так просто дифференцировать функции, конечно, скучно.
ОДУ 2-го порядка - это система ОДУ 1-го порядка, где, скажем, dy/dx = u, а d^2y/dx^2 = du/dx. Ну и так далее. Это если мы именно про ОДУ, а не про эллиптические/гиперболические/параболические уравнения. Разве нет?
я сильно удивлен интересом к моей маленькой заметке
На фоне массы постов про то как "тестировать фронтенд" и "как джуну найти первую работу" просто как глоток свежего воздуха :)
Пока ОДУ. Для ОДУ 2-го порядка (и далее) в Питоне почти все про задачу Коши. А вот с краевыми задачами, хотя бы с условиями Дирихле что-то непонятное и очень накрученное. Я пока код гоняю, смотрю получится ли что-то содержательное.
Не уверен, что Хабр - подходящая площадка для УЧП. Будем двигаться постепенно, у меня есть старые наработки.
А по-поводу тематики понятно, всем веб-программирование и ИИ нужно. Кстати, очень много просмотров получает науч-поп.
Мы решали похожую задачу методом экстраполяции тренда. На границах ряда это позволяет избежать укорочения сигнала, а внутри сигнала - единообразно с границами обрабатывать перерывы и пропуски. Пример того, как это работает для одного частного случая, можно найти вот
в этой статье
Ссылка на статью, а ее текст можно взять вот здесь или вот здесь.
Но, это все сделано не на питоне, а в собственном
пакете для анализа временных рядов
Описание и программы можно скачать вот отсюда, а исходные тексты лежат вот здесь. Но у нас достаточно
специфические данные и задачи
ряды только эквидистантные, но с кучей пропущенных наблюдений, а обработка в скользящем окне обычно включает несколько итераций разными методами. Поэтому длина сигнала при использовании алгоритмов скользящего окна не должна уменьшаться
и пакет заточен строго под их решение, поэтому его сложно назвать универсальным.
К тому же это все написано
на кондовом фортране
с кучей зависимостей между разными функциями, так что вытащить из него какие-то отдельные процедуры проблематично. А пытаться прицепить весь пакет -
безнадежно
как сказал один наш коллега, заглянув внутрь - "у них там свой Excell"
.
Да и идеология у нас немного другая: наш пакет - это фактически не готовые решения, а скорее конструктор, из которого надо самому собирать алгоритм, подходящий к той или иной ситуации.
только вот
сборка этого алгоритма происходит не методом написания кода, а путем выполнения нужных процедур
в интерфейсе
по принципу "график" - "вычисления" - "график"
в результате чего получается не программа, а в лучшем случае макрокоманда на очень примитивном собственном языке. Сейчас мы бы уже сделали по-другому, но ноги всей этой системы растут из 1980-х годов, когда даже такая архитектура была откровением. Так с тех пор и прикручиваем к старому ядру всякие финифлюшки по принципу "работает-не трогай!"
В общем, конкретные функции оттуда вытащить сложно, а вот идеи - пожалуйста.
А почему у вас на последнем графике линия Test с изломом?
ЗЫ: и чтоб два раза не бегать – s/окресности/окрестности/
Как только из аналитической функции получен массив значений, о ней (этой функции) забываем. Работаем только с дискретными значениями. В последнем графике помещаются только три узла, поэтому виден излом. Вообще-то, лучше точки поставить, но тогда не очень красиво будет. Интерполяции между точками пока не касаемся. Без пропавшей согласной писал интуитивно. Наверное, когда в школе учился, можно было и так и так. А сейчас, наверное, лучше попрапавить. Спасибо.
Не надо забывать о функции: через те же точки проходит бесконечное множество кривых, с производными по вкусу (см. интерполяционные полиномы – правда, я не помню, в какой их разновидности можно производные задавать, но она есть... Ну или обычный b-сплайн).
Т.е. для вашего набора точек любое из решений может оказаться верным, поэтому есть резон указать ограничения на функцию (для стандартных методов – обычно это ограничение на модуль второй производной)
Ну и задача с шумом, которую тут упомянули, ещё интереснее.
Согласен. Верных решений много, выбор зависит от разных факторов, например от вычислительных затрат, или от степени гладкости и т.д. В следующей публикации напишу еще об одной интересной задаче, недавно опубликованной на Хабре (если я ее, конечно, правильно понял). На самом деле каждая конкретная задача требует своего подхода. Может быть неправильным ограничивать модуль второй производной.
С модулем второй производной – просто очень легко доказывается, что погрешность вычисления производной будет не больше максимума этого самого модуля на "рабочем отрезке", умноженного на длину того самого отрезка и какую-то константу. Вот этим и пользуются. Это по трём точкам. По пяти – вылезает, соответственно, четвёртая производная (так что считай вы по пяти точкам – последний график бы чётко лёг на эталон) и куб длины отрезка.
Но при наличии шума всё опять хуже, интерполяция высоких порядков начинает проигрывать. Вроде тут в идеале нужен плюс-минус тот же подход,что в фильтре Калмана, когда изначально у нас сигнал и шум, шум оцениваем и аппроксимируем с его учётом... Ну а по-простому – интерполяция низких порядков, да производные по отрезку подлиннее, или заранее сгладить исходные данные. Интересно, что у вас в статье будет.
Наилучшая аппроксимация в окрестности точки (тут, конечно, по пяти точкам точнее будет) - не всегда основное требование. В методе конечных элементов, о котором напишу в продолжении этой публикации, основная идея - наилучшая аппроксимация "в среднем" по всему отрезку. Кроме этого будет рассмотрена одна случайно образовавшаяся задача, где как раз надо принимать решение - какое ограничение накладывать на вторую производную по всему отрезку. Надеюсь, Вас также заинтересует.
Конкретно в этой публикации можно было сразу начать с 5-точечной схемы, также придумав какой-нибудь костыль по краям, там уже по 2 точки "отвалятся". Но это все-таки именно костыли, что неправильно идеологически.
О, метод конечных элементов – прям как маслом по сердцу, студенческие годы вспомнил (кафедра гидрогазодинамики, там числяка много было).
Надеюсь, вы его в приложении к условно-реальным задачам опишете – когда на элемент сразу накладываются законы сохранения.
Понимаете, какая штука. Дифференцирование функции одной производной хоть МКР хоть МКЭ - не слишком увлекательное занятие. Я решил на эту тему написать публикацию только потому, что обнаружил, что в библиотеках Python нет соответствующей адекватной процедуры. Почти то же самое с ОДУ. По задаче Коши много чего, а по краевым задачам немного и невнятно. Но основной интерес, конечно в 2- и 3-мерных задачах. Я в свое время много чего "в стол" написал, кстати в МДТТ. Вот тут самое интересеное должно быть.
А, и ещё: будет не лень – проверяйте методы на функции exp(-1/x²). Очень интересная функция, непрерывная, гладкая, но не раскладывается в ряд Фурье (точнее, Мак-Лорана – в окрестности нуля).
слышали про теорему Такенса
Нет, на практике я с такими задачами не встречался. Как я понимаю, идея там простая и сводится к анализу траекторий в фазовом пространстве с учетом их возможного искажения шумом. Но на практике это будет работать только для стационарных сигналов. А почти все без исключения экспериментальные временные ряды - нестационарны. Начиная от геофизики и кончая макроэкономикой и социологией. Если об этой нестационарности "подзабыть", можно сделать массу интересных открытий ;-) И это относится не только к корреляциям, но и к любым другим методам.
Да, мы всегда пытаемся повысить отношение сигнал/шум, накапливая данные (и неявно при этом предполагая, что полезный сигнал повторяется, а шумы случайны). Но делать это бессмысленно, если измеряемый сигнал каждый раз относится к разным явлениям, а не к одному и тому же. Никто ведь не будет сегодня измерять суточный ход температуры у гриппозного больного, завтра у обезьянки, послезавтра у крокодила, потом в комнате, а потом это все осреднять. Все понимают, что эти явления разные. Однако почему-то считается совершенно нормальным измерять один и тот же показатель на протяжении многих лет, и без каких-либо проверок постулировать априори, что раз мы ежедневно измеряли одно и то же, то значит вправе наши измерения осреднять, вычислять дисперсии-корреляции и т.д. Хотя любая нестационарность изучаемого процесса совершенно четко и однозначно указывает:
нет, это не одно и то же!
На самом деле при повторяющихся измерениях очень важно учитывать: изменилось ли то окружение, в котором протекают наши процессы. Если оно стало другим - то мы измеряем уже какую-то иную величину, хотя и датчики, и вольтметры по-прежнему стоят под той же самой горой. То есть, имеем дело со случайным процессом, а значит, уже не вправе применять для анализа этих рядов аппарат, разработанный для случайных величин.
На первый взгляд, это может показаться буквоедством и формалистикой, но вы удивитесь, к каким катастрофическим последствиям может приводить игнорирование этого факта. За примерами снова отсылаю вот к этой статье.
Но если у нас не случайная величина, а что-то другое, то вопрос по возможности применения большинства типовых методов статистического анализа для таких данных требует особого изучения. Так как говоря строго, условия их применимости уже однозначно нарушены. И пока мы не выясним, какие последствия из этого могут быть, любые попытки содержательного анализа результатов такой обработки фактически беспредметны.
P.S. Точнее, в лаборатории стационарные временные ряды, наверно, бывают. Но вот у нас в геофизике я с ними за много лет еще ни разу не сталкивался.
На практике приходится работать с зашумленными данными. Подсчет производной вышеприведенными способами приводит к усилению шума, причем довольно существенно. Я использую следующий способ - считаю линию в окрестности текущего значения по методу наименьших квадратов и в качестве производной беру тангенс угла наклона этой линии (k из уравнения y=kx+b). Окрестность подбирается под наибольшую частоту спектра исходной функции без учета шумов.
Спасибо за важный комментарий. Обработка зашумленных данных - отдельная важная задача. Часто бывает не до производных, а просто надо "сгладить" экспериментальные данные. Я бы к этому вопросу подходил системно: определил базис функций, на котором будет проводиться аппроксимация, определил меру (например, 
), а после этого Вы автоматически попадете в сопряженные аппроксимации, как наилучшие. Тут есть опасность "переборщить" со сглаживанием и потерять физику процесса. Я скоро выпущу часть вторую, коснусь темы МКЭ. Возможно, пройдусь по теории. Метод наименьших квадратов, конечно, имеет право на существование, но есть много альтернатив.
Как правильно дифференцировать дискретные функции (Часть 1. Тестируем и улучшаем Numpy)