?

Log in

No account? Create an account
Previous Entry Share Next Entry
(no subject)
sky
northern_wind
Мне подбросили задачку, я от нее в некой растерянности.
Есть цель сделать красивую рентгенограмму из плохой, негодной, шумной рентгенограммы. В идеале стоило бы усреднить несколько измерений, но такой возможности по какой-то причине нет (например, измеряемые процессы настолько быстры, что не получается накопить измерений на хороший сигнал).

Человек предлагает сделать следующую имитацию нескольких измерений:
Считается, что интенсивность Y в точке X - это истинная величина плюс-минус некая ошибка, зависящая от значения Y. Ошибка принимает значение (-y0,+y0). Из этого интервала выбирается некое случайное значение (распределение ошибки считается гауссовым), делается поправка для интенсивности (получается "как бы" еще одно измерение), для новой интенсивности процедура повторяется, и так несколько раз. Random walk, короче.
Затем для полученных значений интенсивности считается среднее арифметическое. С ним дальше и работаем.

Мне кажется, что так нельзя. Может появиться ошибка, которая будет накапливаться.
Но насколько нельзя? Как будет выглядеть эта ошибка?
Если чуть-чуть нельзя, то, может быть, стоит пожертвовать ошибкой в пользу такого метода фильтрации шума?
Или в силу случайности процесса ошибка накапливаться не будет? Все же это не последовательные приближения.

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

  • 1
Если что, то могу дать дополнительные комментарии по задаче.
У нас многое нельзя в рентгене, но некоторое из того, что нельзя, иногда с оговорками можно, если очень хочется.

Знаешь в чем я подозреваю еще одну проблему? При достаточном количестве измерений значение среднего арифметического все равно, скорее всего, сведется к исходному значению интенсивности. А недостаточное количество измерений типа недостоверно и приведет к еще большему шуму.
Это нужно сесть и доказывать с бумажкой, конечно.

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

В сторону bootstrapping не копала еще?

Неа, это не моя задача и не моя идея.
Но я, честно говоря, не знаю как применить bootstrapping, если у нас (если я правильно понимаю устройство рентгена) не набор независимых измерений одной величины, а набор измерений разных величин вообще (сейчас призову автора алгоритма в тред, пусть он меня поправит).

Именно так, одна рентгенограмма - набор независимых измерений в разных точках Х разных величин Y. Ошибка измерения в одной точке (я о статистической ошибке) не зависит от ошибки в других точках. И не влияет на них, соответственно.

да было какое-то забавное использование bs для noise reduction в сигналах, если вспомню, как, скажу:) или поищу...

Предмет обсуждения))

Image and video hosting by TinyPic

Сигнал одномерный или двухмерный? Это все решается методами обработки сигналов/изображений, которым уже лет как 60. В данном случаи видно что шум локализован в спектральной области, соответственно вполне может хватить тупого линейного стационарного фильтра. Можно попробовать медианный фильтр, он даже проще в реализации и дает хорошие результаты, особенно в присутствии выбросов, берется окно из ближайших значений (x - window_size/2, x + window_size/2) , в нем находится медиана, которая и берется за значение в точке.

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

Сигнал одномерный.

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

Я сейчас сходу не соображу, возможно ли вообще придумать такую зависимость y0 (то есть, фактически, дисперсии ошибки) от Y, при которой этот процесс бы не вредил, но таки вряд ли. Если же взять идеальный случай, где на каждом шаге мы добавляем гаусса с нулевым матожиданием и фиксированной дисперсией sigma^2, совпадающей с дисперсией изначальных измерений, то получаем:

Y[0] = Y[true] + f[0]
Y[1] = Y[0] + f[1]
Y[2] = Y[1] + f[2] = Y[0] + f[1] + f[2]
...
Y[n] = Y[n-1] + f[n] = Y[0] + f[1] + f[2] + ... + f[n]

где Y[true] - это истина, которая где-то рядом, Y[0] - результат собственно измерения, f[k] - набор независимых гауссовых сл.в. с дисперсией sigma^2 и нулевым матожиданием.

Результат таких вычислений:
Y[calculated] = (Y[0] + Y[1] + ... + Y[n]) / (n+1) = Y[0] + (n*f[1] + (n-1)*f[2] + ... + 2*f[n-1] + 1*f[n]) / (n+1)

Y[calculated] - Y[true] = f[0] + (n*f[1] + (n-1)*f[2] + ... + 2*f[n-1] + 1*f[n]) / (n+1) ---- выражение в правой части является гауссианой с нулевым матожиданием и дисперсией (1 + (n+1)*n/2/(n+1)) * sigma^2 = (1 + n/2) * sigma^2

То есть, чем дальше в лес, тем толще Гаусс: после n итераций дисперсия ошибки возрастёт в (1 + n/2) раз

Ну разве что "делается поправка для интенсивности" каким-то образом всех спасёт, если там будет совсем магическая функция, но я в ней сомневаюсь.

Edited at 2014-05-17 10:15 pm (UTC)

Ага, спасибо!
Я как раз этого и боялась.


Спасибо за расчёты!
Да, по сути именно так.
Когда я делаю эксперимент, то при физическом последовательном усреднении, когда я измеряю величину несколько раз в одной точке, величина в итоге стремится к истинной, поскольку ошибка уменьшается.
Математически, видимо, повторить такое не удастся, поскольку истинная величина неизвестна.
Жаль...

Только я там при рассчёте дисперсии наглюкал.
На самом деле у вот этой сл.в.
f[0] + (n*f[1] + (n-1)*f[2] + ... + 2*f[n-1] + 1*f[n]) / (n+1)
дисперсия будет
(1 + (n/(n+1))2 + ((n-1)/(n+1))2 + ... + (2/(n+1))2 + (1/(n+1))2) * sigma2

но суть от этого не меняется: чем дальше в лес, тем толще Гаусс, потому что ошибка на каждом предыдущем шаге Y[k] добавляется и ко всем последующим шагам.

Спасибо! То есть лучше использовать имеющиеся алгоритмы фильтрации...

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

Да, это может сработать.
Так визуально можно вычленить выбросы, например.

А в мой моск вчера врезалась вот эта задачка, и теперь его потихоньку ест: http://ru-xkcd.livejournal.com/57371.html

Впусти Иисуса в своё сердце задачку по резисторы в свой моск!

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

> а выбираем из них самые хорошие, в смысле гладкости например
Вот это довольно интересно.

Еще можно сделать так: предположить закон распределения для ошибки в каждой точке,
Y[measured] = Y[true] + err
например пусть err ~ N(0, s^2) и ошибки в разных точках независимы,
синтезировать новые рентгенограммы по методу Монте-Карло (НЯП, при очередном испытании для каждой точки получаем некую реализацию err и вычитаем ее из Y[measured]) и выбирать те из них, которые более "хорошие": гладкие и т.п.

А еще ошибки измерения в соседних точках могут быть зависимы :) Как зерно на фотопленке, например.

Edited at 2014-05-18 02:12 pm (UTC)

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

Выше писали про фильтры для обработки сигналов: медианный, гаусса и т.д. Они не дали нужный результат?

Фильтрация такого рода экспериментальных данных существует давно. Меня не устраивают некоторые артефакты, которые дают известные методики.

Как задача целиком звучит? Зачем нужна "красивая рентгенограмма"? Может вам надо отличить здорового от не здорового, или больного заболеванием А от не больного этим заболеванием.
Тогда, например, вот вам реализация многомерной с.в.: каждое измерение Y(t_1), Y(t_2)- компонент вектора. У образца (нулевой гипотезы) - вот такое матожидание этой с.в. и такая дисперсия, и можно пробовать проверять гипотезу, вычислять p-value, или может получится принцип максимума правдоподобия (max likelihood) применить.

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

Красивая нужна затем, что из неё более адекватно данные извлекаются. Форма пиков точнее, есть уверенность, что за шумами не скрыт какой-то слабый сигнал.
Это на самом деле не медицина, а физика и химия.

И всё-таки меня одолевает любопытство, что за задача?

Физическая?
Регистрация in situ дифракционных сигналов быстрых процессов. В этой области хороший (в смысле, с достаточно высокой интенсивностью = низкими шумами) эксперимент часто роскошь.

О, точняк! Броуновское движение!

Чего-то я не понимаю. Во-первых, зачем усреднять измерения (при их наличии), если лучше брать median?

Во-вторых, усредненный random walk - это всего лишь аддитивный гауссов шум. Откуда здесь фильтрация?

В-третьих, с точки зрения теории информации, вынимать хоть что_то можно сугубо из spatial information, добавить туда temporal "от балды" невозможно по-моему ну никак.

Короче, обьясните плз, бо пахнет карго-культом симуляции эксперимента.

  • 1