Pull to refresh
2
0
Send message
На всякий случай, на известной страничке www.lfd.uci.edu/~gohlke/pythonlibs уже собранный pywin32 имеется для всех (не сильно устаревших) вариантов питона.
Если Галилей проделает очень много опытов и построит общую гистограмму (по вертикали число опытов, по горизонтали время скатывания), то получит чистый гаусс.
Нет. Гаусс получается, когда «работает» центральная предельная теорема. А вашем гипотетическом эксперименте есть некоторый неслучайный (квази)периодический эффект, который вылавливается, например, спектральным анализом временного ряда.
Если бы вместо одной гистограммы он построил бы множество, гистограммы получились бы серьёзно отличные от гаусса
Именно. И друг от друга тоже отличные. Все правильно Вы говорите — нужно доказать, что распределения в разные периоды времени — разные. А вот автор обсуждаемой публикации ищет и пытается показать наличие сходства распределений с использованием весьма неестественной метрики, да еще такой, что иногда distance(x,y) != distance(y,x). Вот взял бы он критерий Колмогорова-Смирнова (максимальное различие между кумулятивными распределениями), спорить было бы сложнее. Фокус только в том, что на маленьких выборках большие значения критерия статистически не значимы… Та же фигня творится, кстати, и с точностью оценки/сравнения параметров распределений на малых выборках. Ну, Вы же должны понимать, что это фундаментальная проблема, т.е. как данные не обрабатывай, а точности не будет? Только хуже может стать.
А потому такие «несостоятельные» гистограммы, которые приводит автор статьи, вообще-то говоря имеют полное право на жизнь и представляют отдельный интерес.
Вы можете привести хотя бы еще один пример использования такого подхода? Корректного! Я еще раз повторю, что что для дискретных случайных процессов нельзя выбирать границы и число интервалов для построения гистограмм произвольно, как это сделал автор. А последующее 5-кратное сглаживание с использованием прямоугольного окна длиной 4? Эти параметры откуда взялись?

Ну еще про использование гистограмм, построенных из последовательных выборок, Автор использовал серии по 60 измерений. Ну я я сдублировал его файл и отрезал от дубля первые 30 секунд, А потом посмотрел на кросс-корреляцию минутных гистограмм из обоих файлов. Так-вот — пропала корреляция. Как будто сравнивались два совсем разных файла… И хватит на этом.
Добрался я наконец до моей рабочей станции с двумя 5-ядерными X5650 и проверил ваши данные из обоих массивов на наличие «кросс-корреляций».

Но начну я с того, что ваш код, вычисляющий «расстояние» между сглаженными «гистограммами» (строки 182-204 из вашего кода) часто выдает разный результат в зависимости от порядка, в каком ему «выдаются» данные. Например, если рассмотреть 3-ю и 4-ю выборку (индексация, начиная с нуля) из файла 0505.txt, то:

distance between #3 and #4 = 7.98634661639
distance between #4 and #3 = 18.7975453121

Несимметричная мера, это очень серьезная ошибка, то есть теперь все ваши утверждения про наличие всяких циклов длиной в звездный день, год в процессе радиоктивного распада уже можно считать недоказанными. Исправляйте. Ну а я пока в своем анализе буду использовать полусумму (d(x,y)+d(y,x))/2.

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

  • Из первого 48 часового массива отсчетов дектектора (0505.txt) были с использованием Вашего алгоритма получены 2880 «минутных» сглаженных «гистограмм» G1[i];
  • Из второго 48 часового массива отсчетов дектектора (0506.txt) были получены 2880 «минутных» сглаженных «гистограмм» G2[i];
  • Далее для диапазона временных сдвигов dt от -2880 мин до 2880 мин были получены выборки вида distance(G1[i], G2[i+dt) для максимально возможного диапазона индекса i (чтобы не вылететь за границы массив, т.е. от -2879 до 2879).
  • Для каждой такой выборки были посчитаны ее минимум, среднее, максимум и 5%, 50% и 95% квантили;
  • Для этих данных были построены графики их зависимостей от dt


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

Графики указанных зависимостей и их обсужение
Первый график — зависимость максимума выборки расстояний между гистограммами из разных годов в зависимости от сдвига времени dt (тут на самом деле, минимум, среднее и максимум — но смотрите только на красную кривую, увы, шкалы сильно различаются). Видно, что максимум максимален там, где размер выборки максимален. Ну и больше ничего…

Если поменять шкалу по оси Y, то становится видна зависимость среднего расстояния между гистограммами из разных годов в зависимости от того же сдвига dt (синий цвет). Опять — никаких следов годовой кросс-корреляции.

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

И, наконец последний график — 5% квантиль (зеленый цвет) и медиана. И снова, кроме зависимости разброса от длины выборки, ничего нет. С моей точки зрения плоский график 5% квантиля — это самое сильное доказательство отсутствия периодичности. Т.е. если хотя бы 5% диаграмм из одного года были сильно похожи на гистограммы другого года на периоде, то график бы имел минимум. Но его нет.



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

Может быть, вас есть «правильный» код? Покажите его, пожалуйста. Или, хотя бы опубликуйте алгоритм, чтобы его можно было повторить. А иначе окажется, что Вы совсем зря здесь опубликовались…
Для меня было шоком, когда нашелся звездный год. Проверьте его обязательно
Для того, чтобы проверить, хотелось бы понять, что находится в ваших файлах 0505.txt и 0506.txt? В каждом из них 172800 строк, т.е. по 48 часов наблюдений.

  • Один файл — данные типа за первые 48 часов первого года, второй файл — 48 часов второго года?
  • Или первые сутки в каждом файле — один год, вторые сутки — другой?
  • Или еще как-нибудь иначе?

А то я хоть и сделал ваш код более чем в 10 раз быстрее, но все равно, если статистику собирать, не очень быстро получается.
но максимальная сходство функций плотности вероятности неожиданно именно на +369 минуте.
Судя по Вашему коду на гитхабе, Вы очень странным образом оцениваете плотности распределения по исходным данным:

  1. Сначала исходная последовательность количества отсчетов детектора в секунду делится на выборки по koeff=60 (целочисленных) отсчетов;
  2. Потом для каждой такой выборки строится гистограмма из ko=60 столбиков;
  3. Далее гистограмма сглаживается.

Требование, что число столбиков должно быть много меньше числа отсчетов, очевидно не выполняется. В результате, то, что у Вас получается, гистограммой назвать сложно. Ниже — результат вашей обработки первой 60-секундной выборки из файла 0505.txt:

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

Ну, да, Вы потом эти «зебры» сглаживаете, но где доказательство, что при этом функция распределения не искажается? Ага, хотя бы с использованием критерия Колмогорова-Смирнова. Да еще стоит добавить к сравнению исходное кумулятивное распределение числа отсчетов, до того, как его, возможно, «испортили» описанными выше процедурами.

На всякий случай, «правильная» максимально подробная часовая гистограмма, полученная с помощью
plt.hist(values, bins=numpy.max(values)-numpy.min(values)+1, range=(numpy.min(values)-0.5, numpy.max(values)+0.5))
выглядит так:

максимальная сходство функций плотности вероятности неожиданно именно на +369 минуте.
Про «неожиданно»… Разве это не означает, что для всех других разностей времен формы «минутных» распределений заметно непохожи друг на друга? Критерий Колмогорова-Смирнова это подтверждает? А иначе получается, что Вы берете ряд типа 1,1,1,1,1,… и утверждаете, что каждое измерение совпадает с плюс десятым в списке.

И в конце небольшой философский вопрос. Почему Вы использовали минутные выборки (по 60 значений)? Природа про минуты, впрочем как и другие шкалы измерения времени, не знает. Т.е., что случится с вашими выводами, если размер выборки изменить? Скажем, на 77? Я, после того, как окончательно пойму, что Ваш код делает, проверю.
Ага, только «умирает» он вот так:
Traceback (most recent call last):
  File "D:/My Documents/PycharmProjects/radioactivity/main.py", line 222, in <module>
    stat_distanc(-1439)
  File "D:/My Documents/PycharmProjects/radioactivity/main.py", line 181, in stat_distanc
    else: x=harra[mk+1440-counter];y=harra[mk];
IndexError: list index out of range
Кстати, у вас на github'е в репозитории https://github.com/a-nai/radioactivity только данные, а кода нет. А выдирать его из текста, собирая по кускам и плача над отсутствием форматирования по стандарту PEP8, сил нет. Так-то я в pycharm'е его переформатирую…
Какой общепринятый (видимо, непараметрический) статистический критерий Вы использовали для доказательства этого утверждения?
DTW или Dynamic time warping, это ведь, судя по википедии, один из методов сравнения временных рядов? Причем, как я понял, не очень то и статистический, а больше из области data mining, machine learning etc.?

Объясните тогда, пожалуйста, почему Вы считаете возможным использовать его для сравнения гистограмм, которые временными рядами совсем не являются?

Почему Вы сравниваете именно гистограммы (являющиеся уже сгруппированными и обработанными данными, что может внести искажения) а не сырые исходные данные?

Какие они, кстати? Временные отсчеты детектора? Какого? Как данные оцифровывались и преобразовывались в гистограммы?
По строке поиска «коды abc def» в Google требуемая страничка Россвязи — 4-я сверху.
И вам спасибо за совет не читать статьи по диагонали.
Как забавно получилось, вы молча добавили первый абзац к своей публикации после моего второго комментария, а я в результате получил кучу минусов.

И, да, вы ошиблись, утверждая, что обсуждаемый метод:

Используется по умолчанию в функции optimize из модуля scipy.optimize

Нас самом деле по умолчанию используются: one of BFGS, L-BFGS-B, SLSQP, depending if the problem has constraints or bounds.

Может быть, тогда стоило начать статью с какой-нибудь подобной фразы?
Также, генерация случайной тестовой выборки в виде (см.первая строка второй части программы):

p=np.arange(0,n,1)

это что-то с чем-то. Сгенерируйте «честно». Ну и см. предыдущий комментарий, естественно.
Не увидел в статье, а Linux умеет загружаться с ZFS, как та же самая FreeBSD?
Python 3.x — print не оператор, а функция. Ну, пробел перед скобкой… Строчки «from __future__ import print_function» нет.
Спасибо, как-то прозевал появление этой возможности. Вот еще-бы Pycharm при вызове конструктора такого класса не подсвечивал правильные аргументы как ошибочные… Впрочем, судя по PY-22102, ошибку обещают исправить в 2017.2
И тот человек, который обнаружит, что старая теория не верна, получит Нобелевскую премию, а не будет сожжен на костре.
Лично меня тут сильно раздражает необходимость дублирования имени класса (слева от знака присваивания и в качестве первого аргумента namedtuple). На самом-то деле они могут и не совпадать, например, после рефакторинга в Pycharm, если не поставить галочку насчет поиска в строках и комментариях.

Код ниже рабочий:

from collections import namedtuple
Machine = namedtuple('Car', ['a', 'b'])
machine = Machine(1, 2)
print(machine)


Только выводит:

Car(a=1, b=2)


Т.е. хотелось бы интеграции этого механизма прямо в интерпретаторе, а не «сбоку», как сейчас.

Information

Rating
5,563-rd
Registered
Activity