просто усредненная выглядит нормально

Выглядит нормально или правильное значение? Можешь проверить?

Кол на С читать лень, но я бы проверил, везде ли я взял корень, при работе с центрами масс — где мой 0.

Я обновил слегка пост: питон выдаёт те же значения корреляции, значит мой алгоритм проблем не имеет.

Да, я вот сижу и проверяю опять, где я мог накосячить.

Ещё можно умножить значения на константу для проверки нормирования, если не ошибаюсь.

Да я делил ещё изначально на число частиц, а потом решил убрать: мы все равно нормируем. Получил такой же результат. Другое дело, что как отдельную процедуру проверки я это не рассматривал.

Ну вот, опять какая-то магия: я увеличил число итерации и получил все значения в пределах нормы. Неужели всё же numerical error?

Ну вот, опять какая-то магия: я увеличил число итерации и получил все значения в пределах нормы.

Проверь автокорреляцию через фурье-образ спектральной плотности мощности ©. Там проще обнаружить и удалить постоянную составляющую в спектре.

Поясни?

Поясни?

Дык, это просто разные способы проверить «numerical error». Если нет погрешностей вычисления АКФ через FFT, то мнимая часть FFT будет равна нулю.

P.S. Ещё можно проверить через xcorr в MATLAB.

Аа, так я пробовал np.fft использовать для подсчёта корелляции, там по умолчанию будет нулевая мнимая часть. Юмор в том, что я получаю те же ответы. Матлаба у меня нет, увы, но математика выдаёт схожие результаты (правда нормировка такая же, как и у R, что мне не подходит).

если честно, то я не очень понимаю что требуется, но

немного переписав код, явно вижу, что там ничего и не должно быть <1. Правильно ли, что vx и centre_of_mass_vx - одно и то же? Мы вектора коррелируем?

for (size_t l = 0; l < t_max; l++)
	{
		for (size_t t = 0; t < t_max - l; t++)
		{
			vx2 = vx[t] * vx[t+l];
			vy2 = vy[t] * vy[t+l]; 
			mag  = vx[t]^2   + vy[t]^2;
			magS = vx[t+l]^2 + vy[t+l]^2;

			VACF_com[l]             += vx2 + vy2;
			VACF_com_module[l]      += sqrt(magS * mag);
			VACF_com_orientation[l] += (vx2 + vy2) / sqrt(magS * mag);
		}

		VACF_com[l] /= t_max - l;
		VACF_com_module[l] /= t_max - l;
		VACF_com_orientation[l] /= t_max - l;

	}

	norm1 = VACF_com[0];
	norm2 = VACF_com_module[0];
	norm3 = VACF_com_orientation[0];

Мы вектора коррелируем?

Да, вектора

что там ничего и не должно быть <1.

Ага, не должно, но есть.

Думаю, что дело в данных (tail problem). Даже если одну эту колонку использовать, при той нормировке, что мне нужна, всё равно вылезает больше 1.

7.162431e+00  
7.281370e+00 
7.353923e+00 
7.367500e+00
7.377404e+00 

Ирония в том, что пару лет назад я начал читать про это, но не мог найти приличный код, чтобы понять, как свой писать. И тут я полез в исходники R, где увидел совершенно ужасную вещь:

if autocorrelation > 1.0
autocorrelation = 1.0

Посчитай ручками для l=0...

твои невероятные (не побоюсь этого слова) числа ~7 не содержат никакого криминала.

Мы ж не допускаем версии, что произошло восстание машин и компилятор с процессором именно в твоем случае выдают неверные результат? Верно?

Я, например, не очень представляю что значит «коррелировать вектора» :)

Посчитай ручками для l=0...

результат-то тот же выходит

твои невероятные (не побоюсь этого слова)

что же в них такое невероятное?

Я, например, не очень представляю что значит «коррелировать вектора» :)

Каждый вектор это набор скоростей в определённые моменты времени, очевидно же из определения velocity autocorrelation function.

т.е. если ты руками

вбиваешь 5 раз суммируешь 7.1^2, то получаешь нечто большее, чем 1? Значит твоя «аналитическая» формула не верна.

Вектор скорости (да хоть чего) часто имеет 2 компонетны (а то и все 3). А корреляционная функция для l=0 должна быть числом. Как ты из пары по 2 числа получаешь 1?

Hint: у меня проблем с корреляцией нет; мне так, потрепаться :)

т.е. если ты руками вбиваешь 5 раз суммируешь 7.1^2, то получаешь нечто большее, чем 1? Значит твоя «аналитическая» формула не верна.

А нормировку для кого придумали? Этот элемент и есть общая норма + количество отсчётов.

Вектор скорости

А сама скорость - одну.

А корреляционная функция для l=0 должна быть числом. Как ты из пары по 2 числа получаешь 1?

примерно так же, как из вектора - его модуль.

Суммируешь, нормируешь по l=0, количеству отсчётов (иногда ещё и по количеству частиц), получаешь всегда 1.

слушай, родной

тебе помощь нужна, или повыпендриваться хочешь? Если второе - то я сразу сдаюсь :)

Если первое, то я уже предложил следующее: взять, и сделать руками свою «программу» для l=0. Ты сначала сказал что руками тоже выходит число >1; а теперь спрашиваешь для кого нормировку придумали? Если придумали, почему ты её не применил?

После применения твоего алгоритма у меня выходит VACF_com[0] ~ 7.3^2 + delta VACF_com_module[0] ~7.3^2*1.44 VACF_com_orientation[0] ~ 1/1.44

и, если честно, я не понимаю что ты делаешь; пытаюсь выяснить чего тебе нужно.

А с чего вдруг автокорреляция скорости должна быть нормирована? Это условие очевидно неверно для NVE ансамбля в системе, где идет релаксация

Я ничего не понял, давай ещё раз, с деталями.

А с чего вдруг автокорреляция скорости должна быть нормирована?

Так делают умные дядечки, которые в моей области уже лет 20-30.

Это условие очевидно неверно для NVE ансамбля в системе, где идет релаксация

У меня неравновесные системы

Если ты меня спрашиваешь,

а не другого благородного Дона, то мой совет таков:

а) выяснить, корреляцию какой величины мы хотим посчитать. Если нас интересует модель скорости, то нужно на певром же этапе посчитать массив v=sqrt(vx^2+vy^2) и далее оперировать только с ним; б) выписать на бумажку полученные величины модуля скорости, увидеть, что это величины типа 7.1, 7.2 и т.д. (впрочем, их величина не важна). в) написать на бумажке формулу, какую хотим реализовать. Если это обычная автокорреляционная функция, то нужно сперва посчитать среднее (записать формулу), а потом формулу для искомой величины. г)Cреднее считается типа

for (int i = 0; i < t_max; i++) { av += v[i]; } 
av /= t_max 

 for (int i = 0; i < X3; i++) { num[l] += (v[i]-av)*(v[i+l]-av); } 

 for (int i = 0; i < X3; i++) { denum[l] += (v[i]-av)^2; } 

чисто на удачу:

		VACF_com[lag] /= TLIM / PCOUNT_LIMIT_VALUE - lag;
		VACF_com_module[lag] /= TLIM / PCOUNT_LIMIT_VALUE - lag;
		VACF_com_orientation[lag] /= TLIM / PCOUNT_LIMIT_VALUE - lag;

А вот здесь ты случайно не потерял скобки? Интуиция мне подсказывает что у тебя что-то плохо с нормировкой.

Лично для моих задач всегда имела смысл только ненормированная автокорреляция.

Если нас интересует модель скорости,

и модуль, и вектор

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

Мне нужно считать автокорреляцию как она есть (без вычитания среднего)

Я знаю, как работает алгоритм, потратил немало времени, чтобы разобраться. Меня интересует, почему такие странные значения получаются.

А вот здесь ты случайно не потерял скобки?

не, точно не потерял Нормировка по количеству отсчётов ведётся, просто я состояния сохраняю не каждую итерацию. Здесь

TLIM / PCOUNT_LIMIT_VALUE

lag

Видел другую нормировку, где просто делят на

TLIM / PCOUNT_LIMIT_VALUE

Лично для моих задач всегда имела смысл только ненормированная автокорреляция.

Ну у меня неравновесные системы, где законы сохранения вылетели в трубу, поэтому абсолютные значения теряют смысл. А вот как мера изменения состояния, нормированная автокорреляция очень даже неплохо работает.

Ну так значит ты не понимаешь или не знаешь как нормировать

Я вообще не вижу, с чего она должна быть нормирована, поскольку ни средняя, ни мгновенная скорость частиц в неравновесной системе не сохраняются. Сдается мне, что нормировкой ты называет банальное умножение на шаг интегрирования, которое нужно только чтобы согласовать определение через интеграл с дискретным семплированием по шагу.

Ну у меня неравновесные системы, где законы сохранения вылетели в трубу, поэтому абсолютные значения теряют смысл. А вот как мера изменения состояния, нормированная автокорреляция очень даже неплохо работает.

Классическая автокорреляция корректна только для стационарных процессов.
Для нестационарных используют ARIMA ©.

Прошу прощения, а зачем мне экономические метрики? Я использую то, что люди в моей области используют (и считают корректным).

А какая еще бывает нормировка? Определение для VACF вполне себе чёткое, ни разу не видел чего-то более экзотического.

Прошу прощения, а зачем мне экономические метрики?

Не экономические, а общестатистические. Чтобы метрика соответствовала процессу.

Я использую то, что люди в моей области используют (и считают корректным).

Тогда не удивляйся, что кривизна корректной «Плоской Земли» бывает > 0 :)

Ну если у тебя кроме этого нечего предложить, то я всё же буду следовать примерам умных дядечек, которые сделали себе степени в конкретной области soft matter, а не занимались «общей статистикой».

По-моему, ты путаешь статическое и динамическое определение. Velocity-velocity autocorrelation function можно отнормировать например, на значение при t = 0, которое для системы, состоящей из одинаковых частиц будет пропорционально средней кинетической энергии на частицу. Из физики так же понятно, что средняя кинетическая энергия сохраняется только в равновесной системе — в любом неравновесном состоянии вначале будет наблюдаться перекачка кинетической в потенциальную (или обратно) и только после характерного времени возвращения (man теорема о возвращении Пуанкаре) можно будет сказать, что мы наконец достигли стационарного (не обязательно равновесного!) состояния. Поэтому на коротких временах скорости частиц могут только расти, при этом растет и скалярное произведение v_i(0)*v_i(t) просто потому, что растет модуль v_i(t). Можно, конечно, отнормировать обе скорости на 1, получив таким образом корреляционную функцию направления, но физ. смысл её мне лично не ясен.

Для самой же автокорелляционной функции скоростей есть ясные применения: при расчете коэфф. диффузии, так же её Фурье-образ есть power spectrum системы, который можно считать пропорциональным IR спектру.

(man теорема о возвращении Пуанкаре)

Но она же только для эргодических систем.

Поэтому на коротких временах скорости частиц могут только расти, при этом растет и скалярное произведение v_i(0)*v_i(t) просто потому, что растет модуль v_i(t).

Что-то так и не увидел, почему они могут только расти.

таким образом корреляционную функцию направления, но физ. смысл её мне лично не ясен.

По факту будет показывать изменение поляризации.

Но она же только для эргодических систем.

Ты не на экзамене, студент, должен иметь какое-то ощущение физики, а не просто повторять заученные фразы. В твоей soft matter почти все системы эргодические, поскольку барьеры малы. Ладно бы кристаллы какие, так и там усть методы типа simulated annealing и генетических, которые позволяют переваливать через высокие барьеры, не дожидаясь миллиардов лет на флуктуацию нужной амплитуды.

Что-то так и не увидел, почему они могут только расти.

:facepalm: Ну представь, что у тебя в системе идет экзотермическая реакция. Что наблюдается? Рост температуры = рост средней кин. энергии.

По факту будет показывать изменение поляризации.

ты знаешь это или гадаешь? Какая поляризация в изотропной, скажем, водичке?

Ох, мы с тобой говорим про совсем разные вещи.

В твоей soft matter почти все системы эргодические, поскольку барьеры малы.

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

Ну представь, что у тебя в системе идет экзотермическая реакция. Что наблюдается? Рост температуры = рост средней кин. энергии.

Нет, не идёт.

ты знаешь это или гадаешь? Какая поляризация в изотропной, скажем, водичке?

Знаю, у меня анизотропные частицы.

Ох, мы с тобой говорим про совсем разные вещи.

Мы говорим с тобой о том, на что ты соизволил предоставить информацию — всё остальное есть лишь в твоей голове, а в моей — только выданная информация, плюс некоторая доля экстраполяций, которые возможны и которые я знаю как сделать.

Ты, судя по всему, студент или аспирант (PhD/grad student). Так вот, основная проблема людей на этом уровне профессионального развития — необоснованные предположения, что остальные достаточно квалифицированы, чтобы сходу понять твою проблему. При этом ты забываешь а) сообщить достаточно сведений о собственно проблеме б) сообщить контекст, в котором эти сведения валидны в) выяснить, правильно ли понял собеседик(и), где у тебя затык. Иногда тут добавляется пункт г) про ЧСВ спрашивающего и его фактическое нежелание понять тех, кто пытается тебе ответить, на основании понимания с предыдущих пунктов, потому что цель была не получить ответ, а показать «смотрите какой я умный».

Так вот, если ты надеешься получить конструктив в дальнейшем, либо попытайся выполнить пункты a-в, либо отказывайся от этой идеи и ищи помощи у коллег/руководителя, у которых в той или иной степени представления о этих трёх пунктах имеются.

Я честно попытался вникнуть в твою проблему, сделав той или иной степени обоснованные предположения, но после твоего заявления «здесь все совсем не так как на самом деле» — я умываю руки. Пока мне проблема казалась интересной и было ощущение прогресса — я что-то думал и советовал. Теперь могу только развести руками и надеяться, что ты изволишь дать больше информации по пунктам а-в. Если нет---например, проблема более неактуальна или ты потерял к ней интерес---ну так и суда нет, закончим и всё.

Обоже, но тут то зачем изобретать велосипед? Эти автокорреляции уже давно посчитаны, бери R или чего там ещё выдумали, и не мучайся.