LINUX.ORG.RU

Автокорреляция > 1

 , , ,


1

3

Решил я тут кое-какие расчёты сделать и внезапно начал получать значения нормированной автокорреляции > 1, чего быть не должно. Код выглядит примерно так, проблемы возникают, когда я начинаю считать корреляцию для центра масс, просто усредненная выглядит нормально (алгоритм и циклы те же):

for (size_t lag = 0; lag < TLIM / PCOUNT_LIMIT_VALUE; lag++)
	{
		for (size_t time = 0; time < TLIM / PCOUNT_LIMIT_VALUE - lag; time++)
		{
			vx_squared = centre_of_mass_vx[time] * centre_of_mass_vx[time + lag];
			vy_squared = centre_of_mass_vy[time] * centre_of_mass_vy[time + lag]; 

			magnitude = centre_of_mass_vx[time] * centre_of_mass_vx[time] +
			            centre_of_mass_vy[time] * centre_of_mass_vy[time];

			magnitudeShifted = centre_of_mass_vx[time + lag] * centre_of_mass_vx[time + lag] +
			                   centre_of_mass_vy[time + lag] * centre_of_mass_vy[time + lag];

			VACF_com[lag] += vx_squared + vy_squared;

			VACF_com_module[lag] += sqrt(magnitudeShifted * magnitude);

			VACF_com_orientation[lag] += (vx_squared + vy_squared) / sqrt(magnitudeShifted * magnitude);
		}

		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;

	}

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

Если честно, ощущения, что слегка едет кукуха: отклонения невелики, но на ошибку округления не похоже. Плохих значений вроде нет (хотя нужно проверить). Я что-то упускаю, но не могу понять, что именно.

Upd. Нет, проблема не в алгоритме: я проверил с помощью np.correlate, получил тот же результат. Ирония в том, что начальная точка всегда максимум, такого быть не может. Upd2. Минимальный набор vx, vy для проверки

7.162431e+00  1.572510e-03 
7.281370e+00 -1.435231e-01 
7.353923e+00 -2.122551e-02 
7.367500e+00  1.025319e-01 
7.377404e+00  9.108236e-02 

★★★★★

Последнее исправление: ZERG (всего исправлений: 2)

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

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

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

tyakos ★★★
()
Ответ на: комментарий от tyakos

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

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

ZERG ★★★★★
() автор топика
Ответ на: комментарий от tyakos

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

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

ZERG ★★★★★
() автор топика
Ответ на: комментарий от ZERG

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

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

quickquest ★★★★★
()
Ответ на: комментарий от ZERG

Поясни?

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

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

quickquest ★★★★★
()
Последнее исправление: quickquest (всего исправлений: 1)
Ответ на: комментарий от quickquest

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

ZERG ★★★★★
() автор топика

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

немного переписав код, явно вижу, что там ничего и не должно быть <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];
sshestov
()

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

Да, вектора

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

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

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

7.162431e+00  
7.281370e+00 
7.353923e+00 
7.367500e+00
7.377404e+00 
Я уже обратился к первоисточникам (allen tildesley computer simulation of liquids), результат тот же.

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

if autocorrelation > 1.0
autocorrelation = 1.0
ZERG ★★★★★
() автор топика
Ответ на: комментарий от ZERG

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

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

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

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

sshestov
()
Ответ на: Посчитай ручками для l=0... от sshestov

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

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

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

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

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

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

ZERG ★★★★★
() автор топика
Ответ на: комментарий от ZERG

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

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

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

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

sshestov
()
Ответ на: т.е. если ты руками от sshestov

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

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

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

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

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

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

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

ZERG ★★★★★
() автор топика
Ответ на: комментарий от ZERG

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

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

Если первое, то я уже предложил следующее: взять, и сделать руками свою «программу» для 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

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

sshestov
()
Ответ на: комментарий от ZERG

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

unanimous ★★★★★
()
Ответ на: комментарий от unanimous

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

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

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

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

ZERG ★★★★★
() автор топика
Ответ на: комментарий от ZERG

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

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

а) выяснить, корреляцию какой величины мы хотим посчитать. Если нас интересует модель скорости, то нужно на певром же этапе посчитать массив 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 
Потом сделать массив num (числитель), такой чтобы
 for (int i = 0; i < X3; i++) { num[l] += (v[i]-av)*(v[i+l]-av); } 
и похожий массив denum:
 for (int i = 0; i < X3; i++) { denum[l] += (v[i]-av)^2; } 
е) посчитать num/denum и профигарить цикл по l

sshestov
()

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

		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;

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

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

cvv ★★★★★
()
Последнее исправление: cvv (всего исправлений: 1)
Ответ на: Если ты меня спрашиваешь, от sshestov

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

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

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

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

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

ZERG ★★★★★
() автор топика
Ответ на: комментарий от cvv

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

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

TLIM / PCOUNT_LIMIT_VALUE
это количество сохраненных состояний системы, читай - количество итераций, ну а
lag
он и в африке лаг.

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

TLIM / PCOUNT_LIMIT_VALUE
, но такое только в статистике встречал, во всех областях soft matter/active matter нормируют именно так.

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

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

ZERG ★★★★★
() автор топика
Ответ на: комментарий от unanimous

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

unanimous ★★★★★
()
Ответ на: комментарий от ZERG

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

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

quickquest ★★★★★
()
Ответ на: комментарий от quickquest

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

ZERG ★★★★★
() автор топика
Ответ на: комментарий от unanimous

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

ZERG ★★★★★
() автор топика
Ответ на: комментарий от ZERG

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

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

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

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

quickquest ★★★★★
()
Ответ на: комментарий от quickquest

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

ZERG ★★★★★
() автор топика
Ответ на: комментарий от ZERG

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

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

unanimous ★★★★★
()
Ответ на: комментарий от unanimous

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

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

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

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

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

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

ZERG ★★★★★
() автор топика
Ответ на: комментарий от ZERG

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

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

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

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

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

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

unanimous ★★★★★
()
Последнее исправление: unanimous (всего исправлений: 1)
Ответ на: комментарий от unanimous

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

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

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

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

Нет, не идёт.

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

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

ZERG ★★★★★
() автор топика
Ответ на: комментарий от ZERG

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

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

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

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

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

unanimous ★★★★★
()

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

yvv ★★☆
()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.