Помощь - Поиск - Пользователи - Календарь
Полная версия этой страницы: Абсцисса пересечения двух гауссиан
Форум врачей-аспирантов > Разделы форума > Медицинская статистика
nokh
Для каждого распределения известны параметры мю и сигма, а также объём выборки (в долях единицы).
Поиск по теме дал несколько аналогичных результатов.
Например, здесь дан вывод уравнения для нахождения абсциссы точки пересечения через решение квадратного уравнения:
https://stats.stackexchange.com/questions/3...asiest,2(x)%3D0.
А здесь те же формулы даны для matlab и подходят для R: https://stackoverflow.com/questions/5202142...n-distributions
Здесь на пайтоне: https://stackoverflow.com/questions/4136865...etween-gaussian

Я завёл всё это в Excel - работает (приложил). Но этот подход предполагает равенство объёмов выборок. На практике же они обычно разные и если использовать разделение смеси распределений, то тут эта формула не работает. Я приложил картинку, где реальные данные приближаются тремя распределениями. Пакет mixdist выдал:
Parameters:
pi mu sigma
1 0.09875 1.417 0.9399
2 0.84174 5.608 1.4961
3 0.05951 10.260 1.2689

Используя эти параметры я не могу найти абсциссы пересечения кривых. Например, подстановка мю и сигм в формулу выше даёт для двух первых распределений значение 3,1011, тогда как при имеющемся соотношении плотностей распределений визуально должно быть около 2,4. Ясно, что по мере уменьшения доли первой группы в выборке эта точка будет всё сильнее сдвигаться влево, пока не скатится по левой горке распределения второй группы к нулю (визуально).

Прошу помочь идеями или кодом, как найти искомое. На худой конец наверное можно как-то "выпотрошить" функцию plot, чтобы найти точку двух кривых с одинаковой ординатой и выбрать её абсциссу (хотя не хотелось бы привязываться к конкретному софту, т.к. пакет PAST выдаёт немного отличные параметры).
Диагностик
Цитата(nokh @ 18.02.2022 - 03:33) *
подстановка мю и сигм в формулу выше даёт для двух первых распределений значение 3,1011, тогда как при имеющемся соотношении плотностей распределений визуально должно быть около 2,4.

А доли подставляли?
nokh
Цитата(Диагностик @ 18.02.2022 - 03:35) *
А доли подставляли?

Нет. Я не нашёл куда)))
Игорь
Цитата(nokh @ 17.02.2022 - 22:33) *
Для каждого распределения известны параметры мю и сигма, а также объём выборки (в долях единицы).
... подход предполагает равенство объёмов выборок

1. Численности в формулы не входят - равенство объемов не требуется.
2. Формулы для нормальных кривых. Если распределение отличается от нормального (в том числе двухвершинное, возможно, что смесь), нужно выводить новые формулы.
3. Формулы основаны на логарифмировании - понятно, что нужно развалить функцию на простые составляющие. Но при этом не говорится, что применяя операцию взятия логарифма, мы накладываем ограничение на область определения функции.
4. Методика не годится для эмпирических функций распределения.

Теоретические функции распределения - не такие уж сложные математические объекты. Во всяком случае, их формулы и свойства известны. Можно попробовать сразу численно решать нелинейное уравнение f1 - f2 =0.
Диагностик
Цитата(nokh @ 18.02.2022 - 15:54) *
Я не нашёл куда)))
Домножить плотность функции распределения каждую на свою долю в смеси.
100$
Цитата(nokh @ 17.02.2022 - 22:33) *
Используя эти параметры я не могу найти абсциссы пересечения кривых.


Ну, это вы зря )

Мое решение x*=2.4289
Решение прикрепил.
В вашем исходном файле были кое-какие ошибочки (см. примечания в соответствующих ячейках)

Для корректного определения x* нельзя забывать о весах распределений в смеси. Диагностик дело говорит.

nokh
Спасибо всем огромное! Сегодня утром с подачи Диагностика расширил член С в уравнении логарифмом отношения выдаваемых пи, но искомого значения не получил - видимо из за ошибок которые нашёл 100$. Буду разбираться...

Это мы со студенткой ищем естественные границы для выделения групп в очень асимметричном распределении.Преобразовали по Боксу - Коксу, разделили смесь, а с границами я споткнулся... Далее ретрансформируем границы в исходную шкалу показателя (наивным обратным преобразованием Б-К по выписанной лямбде) и будет полезный на практике результат. Кстати, уже не первый раз сталкиваюсь с такой задачей, просто раньше точного решения не нужно было, можно было просто по ядерной плотности. Удивительно, что в сети нет готовых решений, по крайней мере быстро не находятся...

Ну и бонусом медицинскому форуму - байка из жизни. Давным давно, когда я ещё не знал про разделение смеси распределений, мы с гинекологом искали критерии для отнесения пациенток с гарднеререллой к группе больных гарднереллёзом, поскольку собственно заболевание не у всех носителей, у части - как условно-патогенный м/о. Там у врачей свои тонкости диагностики, но они не всегда срабатывали. По результатам факторного анализа в первый фактор вошли как раз нужные показатели (гарднерелла, ключевая клетка и т.п), а ненужные не вошли. Факторные метки первого фактора дали красивое бимодальное распределение (там около 250 наблюдений было) с тонюсенькой зоной трансгрессии, т.е. собственно больных Г и не больных. Я ещё тогда с формулами игрался, но до ума не довёл, а гинеколог кандидатскую бросила, так что даже публикации не осталось.
Диагностик
nokh, дайте данные по гистограмме, попробую вашу смесь расщепить.
nokh
Цитата(Диагностик @ 18.02.2022 - 16:38) *
nokh, дайте данные по гистограмме, попробую вашу смесь расщепить.

прикрепил
100$
Цитата(nokh @ 18.02.2022 - 18:39) *
прикрепил


А это уже Боксо-Коксовая цифирь? Похоже, что так.
И зело странно, что на 57 с.в. на гистограмме аж 14 разрядов.
nokh
Цитата(100$ @ 18.02.2022 - 20:45) *
А это уже Боксо-Коксовая цифирь? Похоже, что так.
И зело странно, что на 57 с.в. на гистограмме аж 14 разрядов.

Да, уже преобразованные. Ну да, 14 разрядов не по Стургесу))) Зато неоднородность видна хорошо, а то если всё будет одной группой не понятно почему кривая плотности изгибы даёт...
Ваше решение не понял, похоже там хитрость) Похоже на то, что параметр X находился функцией подбора значений. Ну или офис у меня дома старый (2003), т.к. подстановка моего значения в ячейку для Х вызвала сбой в других.
Но это уже не важно, я добавил веса в приравниваемые уравнения f(x) и получил искомое X аналитически. С учётом исправлений в первом файле LOG на LN и добавлением в параметр С квадратного уравнения логарифма отношения весов всё заработало!

Поэтому экселевский файл в первом сообщении убираю и заменяю на правильный. Ещё раз спасибо за наводки!

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

100$
Цитата(nokh @ 18.02.2022 - 19:15) *
Да, уже преобразованные. Ну да, 14 разрядов не по Стургесу)))


Мне гистограмма по Симадзаки-Синомото сразу выдала три больших столбца.
Файл прикрепил. Гляньте, когда не лень)

Цитата
Похоже на то, что параметр X находился функцией подбора значений.


Конечно. Точно так.



Диагностик
Цитата(nokh @ 19.02.2022 - 00:15) *
Да, уже преобразованные.

А можно исходные данные?
Диагностик
Цитата(nokh @ 19.02.2022 - 00:15) *
Да, уже преобразованные. Ну да, 14 разрядов не по Стургесу))) Зато неоднородность видна хорошо, а то если всё будет одной группой не понятно почему кривая плотности изгибы даёт...

Неоднородность есть следствие случайности выборки. Попробовал 19 разрядов. По критерию Пирсона не отвергается гипотеза о нормальности статистического распределения на уровне 0,25.
100$
Цитата(Диагностик @ 19.02.2022 - 03:36) *
Попробовал 19 разрядов.


Кто больше?

Цитата
По критерию Пирсона не отвергается гипотеза о нормальности статистического распределения на уровне 0,25.


И на это ушло 2 ч. 15 мин.?

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

И вообще - по ночам спать надо...
Диагностик
И на это ушло 2 ч. 15 мин.?

Это не важно. Важно то, что это не смесь.
100$
Цитата(Диагностик @ 19.02.2022 - 14:28) *
Это не важно. Важно то, что это не смесь.


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

Однако.
Нормальное распределение устойчиво по суммированию: сумма нормальных величин - нормальна.

Поэтому, если у исследователя есть причины предполагать, что здесь смешаны 3 распределения - пусть моделирует.
Другое дело, что нельзя складывать до кучи 3 распределения, не убедившись предварительно с помощью теста отношения правдоподобия, что это не противоречит имеющимся данным.
nokh
>Диагностик

Конкретно здесь неоднородность не является следствием случайности выборки, т.к. это не совсем исходные данные, а это - важно. Данные получены преобразованием исходных с помощью адаптивного к данным преобразования Бокса - Кокса, которое и предназначено для того, чтобы делать исходные распределения нормальными настолько, насколько это только возможно. Поэтому то, что тестами не обнаруживается отличие от нормальности указывает лишь на то, преобразование справилось со своей задачей. Но вот то, что даже после такой замечательной штуки плотность распределения указывает на полимодальность и является основанием предполагать смесь распределений. Теоретически это тоже оправдано. Я много работаю с преобразованием БК, т.к. в тех областях где приходится считать чаще всего нормального распределения почти не бывает, а уж так сложилось, что я люблю среднее и ДИ больше медианы с квартилями, и это любовь не иррациональная, а обоснованная практикой. Поэтому часто приходится видеть и унимодальные и полимодальные распределения после БК. Суть проводимой работы - отыскание естественных границ для разных классов объектов в том случае, когда исходные распределения настолько асимметричны, что не позволяют даже предположить неоднородность данных.

Для любителей поюзать реальные данные выложил файл целиком: донные отложения только озёр (есть загрязнённые), валовое содержание 4-х элементов, исходные значения (мг/кг сухого вещества) и преобразованные. Могу какие-нибудь цитокины поискать, там тоже всё сильно асимметрично.

Судя по всему, в этой ситуации у нас нет критериев, чтобы оценить статистическую значимость решения. Приходится полагаться на теоретическую возможность, глазомер и вспомогательные процедуры типа плотности распределения, всяких BIC и AIK, японской диаграммы, по которую узнал от 100$ (о ней ниже)
nokh
> 100$

Японская диаграмма прикольная. Я люблю такие автоматизированные техники. Оптимизированный сплайн, аддитивные модели регрессии, да то же преобразование БК. Нашёл код для R здесь:
https://web.archive.org/web/20210909021638/https://www.neuralengine.org/res/histogram.html
Пока про принципы не читал, не очень понравилось, что на данных первого примера несколько сузила средний класс откинув его крайние варианты в крайние классы. Рис. прикрепил. В структуре результата R есть границы классов - полезно. Почему-то не справляется с Zn после БК: выдаёт обычную гистограмму. Попробуйте свой экселевский код, может получится?

md<-read.table("clipboard", dec=",")
str(md)
'data.frame': 57 obs. of 1 variable:
$ V1: num 10.27 4.49 7.77 12.08 7.95 ...

sshist <- function(x){
N <- 2: 100
C <- numeric(length(N))
D <- C
for (i in 1:length(N)) {
D[i] <- diff(range(x))/N[i]
edges = seq(min(x),max(x),length=N[i])
hp <- hist(x, breaks = edges, plot=FALSE )
ki <- hp$counts
k <- mean(ki)
v <- sum((ki-k)^2)/N[i]
C[i] <- (2*k-v)/D[i]^2 #Cost Function
}
idx <- which.min©
optD <- D[idx]
edges <- seq(min(x),max(x),length=N[idx])
h = hist(x, breaks = edges )
rug(x)
return(h)
}

res<-sshist(md$V1)
str(res)
List of 6
$ breaks : num [1:4] 0.936 4.651 8.365 12.08
$ counts : int [1:3] 13 38 6
$ density : num [1:3] 0.0614 0.1795 0.0283
$ mids : num [1:3] 2.79 6.51 10.22
$ xname :8322456 "x"
$ equidist: logi TRUE
- attr(*, "class")=8322456 "histogram"
100$
Цитата(nokh @ 19.02.2022 - 21:21) *
Почему-то не справляется с Zn после БК: выдаёт обычную гистограмму. Попробуйте свой экселевский код, может получится?


Обожаю все эти контаминанты: такая гадость )
Результаты в прикрепленном файле: гистограммная оценка + непараметрическая (вся информация - в заголовках диаграмм)
- сырых данных (Zn);
- трансформированных (Zn_tr)

Интересно было бы подогнать модельное распределение именно к сырым данным: оно похоже на распределение Рэлея.
И честное слово, глядя на ядерные оценки плотности распределения сырых данных, я не вижу в них никакой "неоднородности" (ширина окна для оценивания плотности оптимизирована методом максимума правдоподобия).
Диагностик
Цитата(nokh @ 20.02.2022 - 01:57) *
Для любителей поюзать реальные данные выложил файл целиком: донные отложения только озёр (есть загрязнённые), валовое содержание 4-х элементов, исходные значения (мг/кг сухого вещества) и преобразованные. Могу какие-нибудь цитокины поискать, там тоже всё сильно асимметрично.

Каждое значение св это концентрация элемента для отдельного озера из 59?
nokh
Цитата(100$ @ 20.02.2022 - 01:49) *
...И честное слово, глядя на ядерные оценки плотности распределения сырых данных, я не вижу в них никакой "неоднородности" (ширина окна для оценивания плотности оптимизирована методом максимума правдоподобия).

Да, поэтому и "Суть проводимой работы - отыскание естественных границ для разных классов объектов в том случае, когда исходные распределения настолько асимметричны, что не позволяют даже предположить неоднородность данных" )))
Мы же не ищем лёгких путей...
Вы приложили старый файл, там где Pb

Цитата(Диагностик @ 20.02.2022 - 03:32) *
Каждое значение св это концентрация элемента для отдельного озера из 59?

Да. 59 озёр.
Диагностик
nokh, нужно найти аномальные значения концентрации?
nokh
Цитата(Диагностик @ 20.02.2022 - 11:27) *
nokh, нужно найти аномальные значения концентрации?

Вопрос можно понять двояко, так и отвечу.
1) нужно найти аномальные значения концентрации (в дополнение к задаче)?
Ну, не то что нужно, но они возможны (сильные загрязнения) и если такие аномальные значения - явные выбросы, то они ухудшают общую картину и подход в целом, т.к. преобразование пытается "поджать" и их. Поэтому скорее да, не помешает.
2) заключается ли задача в том, чтобы найти аномальные значения концентрации?
Нет, задача заключается как раз в том, чтобы найти границу, отделяющую условно фоновые значения концентраций от всех остальных. Т.е. это типа ПДК, критерия для нормирования, чтобы можно было сказать "раз значение больше ..., значит есть основания подозревать загрязнение". В отличие от воды, ПДК для донных отложений не разработаны, но это отдельная тема не для этого форума. Условность понятия "фоновые" связана с тем, что все водоёмы в той или иной степени техногенно загрязнены. Но даже если рядом нет заводов, просто живописное место, то там очень высокая автотранспортная нагрузка от отдыхающих, т.е. паттерн Zn+Cd+Pb, есть основания подозревать загрязнение карповых водоёмов свинцом от грузил рыбаков и т.п.
Диагностик
nokh, аномальные значения левого "хвоста" гистограммы о чём-то говорят?
Диагностик
nokh,поработал со свинцом. Оказалось что концентрация у него распределена логнормально.
Нажмите для просмотра прикрепленного файла
Нашёл параметры ограниченного нормального распределения (отрубил предполагаемые аномальные значения). Потом, добавляя по одному интервалу, каждый раз проверял гипотезу о нормальности и пересчитывал параметры нового распределения. Прошёл таким образом весь правый хвост и в результате аномалий не обнаружил. На левом хвосту аномалии есть (5 озер), но я их не учитывал при подборе распределения. Затем нашёл допустимое значение по распределению крайней статистики для 59-го озера на уровне 0,95. Получилось 111,7. Далее представил своё видение оценок загрязнения в зависимости от концентрации.
Нажмите для просмотра прикрепленного файла
Проверил таким-же образом трансформированные данные по Б-К, справа выбросов также нет, слева 7 озер.
Нажмите для просмотра прикрепленного файла
Диагностик
Цитата(Диагностик @ 21.02.2022 - 08:37) *
nokh,поработал со свинцом. Оказалось что концентрация у него распределена логнормально.
Принципиальная ошибка. Нельзя с исходной величиной совершать преобразования, приводящие к снижению неоднородности исходного распределения. Подобные логарифмированию и Б-К. При этой процедуре выбросы маскируются, а мы стремимся их наоборот, выявить. В связи с вышесказанным провел анализ чистых исходных данных для каждого элемента. Вот результаты:
Нажмите для просмотра прикрепленного файла
Нажмите для просмотра прикрепленного файла
Нажмите для просмотра прикрепленного файла
Нажмите для просмотра прикрепленного файла
Нажмите для просмотра прикрепленного файла
Нажмите для просмотра прикрепленного файла
Нажмите для просмотра прикрепленного файла
Нажмите для просмотра прикрепленного файла

...
Диагностик
Цитата(Диагностик @ 24.02.2022 - 11:18) *
Принципиальная ошибка. Нельзя с исходной величиной совершать преобразования, приводящие к снижению неоднородности исходного распределения. Подобные логарифмированию и Б-К. При этой процедуре выбросы маскируются, а мы стремимся их наоборот, выявить.
Однако случай со свинцом и остальными элементами показал, что можно и даже нужно. Займусь на досуге.
nokh
Цитата(Диагностик @ 24.02.2022 - 08:18) *
Принципиальная ошибка. Нельзя с исходной величиной совершать преобразования, приводящие к снижению неоднородности исходного распределения. Подобные логарифмированию и Б-К. При этой процедуре выбросы маскируются, а мы стремимся их наоборот, выявить. В связи с вышесказанным провел анализ чистых исходных данных для каждого элемента. Вот результаты:

Благодарю за интерес к проблеме и труд. Разгребаю другие дела, поэтому пока посмотрел не вникая, но потом погляжу повнимательней. В принципе, то что я делаю - я уверен в работоспособности такого подхода. Но к критике нужно быть готовым - я планирую это публиковать (с добавлением биологических примеров). Поэтому буду признателен за ссылки, где такие мнения (нельзя ... и т.д.) прописано, чтобы вступить в виртуальную полемику))) В принципе никто не подвергает сомнению то, что полимодальность указывает на внутреннюю неоднородность данных. А вот то, как с этой неоднородностью работать и как на неё выходить - нет готовых рецептов и то, что делаете с данными вы является одним из возможных подходов, уже вашей наработкой.
Диагностик
Цитата(nokh @ 25.02.2022 - 18:20) *
как с этой неоднородностью работать

Работал с исходными данными замеров (непреобразованными). Отсекал подозрительные на выброс крайние элементы выборки. По ММП находил параметры исходного распределения (ограниченный нормальный закон). Проверял эту гипотезу. ПДК находил через предикционные интервалы по ГОСТ Р ИСО 16269-8-2005. Получил следущие значения:

никель - 94,5; выбросов 21 шт.
медь - 49,0; выбросов 21 шт.
цинк - 235,2; выбросов 5 шт.
свинец - 60,7; выбросов нет.
100$
Цитата(Диагностик @ 4.03.2022 - 04:28) *
Работал с исходными данными замеров (непреобразованными). Отсекал подозрительные на выброс крайние элементы выборки. По ММП находил параметры исходного распределения (ограниченный нормальный закон). Проверял эту гипотезу. ПДК находил через предикционные интервалы по ГОСТ Р ИСО 16269-8-2005. Получил следущие значения:

никель - 94,5; выбросов 21 шт.
медь - 49,0; выбросов 21 шт.
цинк - 235,2; выбросов 5 шт.
свинец - 60,7; выбросов нет.


Подведем некоторые промежуточные итоги.

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

Для того, чтобы прервать размещение юзером "Диагностик" зубодробительной чуши, необходимо произнести заклинание "Раз,два,три - горшочек, больше не вари".

Позвольте пояснить.

1. Во-первых, понятие "выброс" не существует само по себе вне привязки к соответствующему гипотетическому распределению.

Поэтому существуют методы для (априорно) нормального распределения, экспоненциального, Вейбулла - Гнеденко.

В качестве исключения могу вспомнить разве что критерий Дарлинга, детектирующий по одному выпадающему наблюдению (max или min) в случае, если про распределение известно лишь то, что оно непрерывно.

2. Если есть подозрение, что выбросов несколько, то их число должно быть заранее известно. Тогда можно воспользоваться критерием Титьена - Мура (Tietjen G., Moore H.,1972).

Если же число выбросов заранее неизвестно, и применяется некоторая последовательная процедура детектирования выбросов, то необходимо заранее знать ее статистические свойства: прежде всего, способность удерживать фактический уровень значимости вблизи заявленного номинального. Единственная известная мне процедура, удовлетворяющая данному требованию, - критерий Роснера (Rosner B., 1975, 1977).

Бесчисленное применение к одной и той же выборке одного и того же метода детектирования выпадающих наблюдений - грубая ошибка.

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

2. Здесь не может быть никакого разглагольствования о ПДК, поскольку

а) понятие ПДК неотделимо от событийного ряда, отвечающего на вопрос "Что будет, если..."
б) Существует дискретный параметр - "Класс загрязненности воды" (5 классов). Поскольку nokh упомянул, что есть загрязненные озера, надо учитывать эту информацию.
в) ГОСТ Р 59 054 - 2020 вводит дискретный параметр - "Классификация водных объектов по целям водопользования" (п. 4.2., табл. 2). И для озер хозяйственно-питьевого назначения ПДК будут иными, нежели для рыбопромысловых, etc.
г) необходимо учитывать антропогенную нагрузку на водоем (и много чего еще).

Поскольку неизвестно, что представляют собой эти 59 озер (сплошное наблюдение или выборочное обследование), попытки отделить мух от котлет "фоновое от нефонового" при такой постановке вопроса - несостоятельны. Особливо принимая во внимание, что представленные данные
- вообще-то многомерные, ergo надо проверять многомерную нормальность, применять многомерного Бокса-Кокса, детектировать многомерные выбросы и т.д.
- вообще-то пространственно-распределенные.
Диагностик
Я не детектировал последовательно каждый выброс, а совсем наоборот. Цензурировал исходную выборку, отбросив с запасом подозрительные на выброс значения.
Пусть в ней осталось n элементов. По ним находил параметры гипотетического распределения (нормального, ограниченного) и проверял максимально возможное расчетное значение крайнего (n+1)-го элемента. Если это значение оказывалась больше реального, возвращал этот элемент в выборку. И т.д. для следующего (n+2)-го.
Олег Кравец
Moderator on

Коллеги, уважайте себя и собеседников. Свалок сейчас и так в Инете полно.

Тема маленько расчищена.
Для просмотра полной версии этой страницы, пожалуйста, пройдите по ссылке.
Форум IP.Board © 2001-2025 IPS, Inc.