Форум врачей-аспирантов

Здравствуйте, гость ( Вход | Регистрация )

> Абсцисса пересечения двух гауссиан, для разных N
nokh
сообщение 17.02.2022 - 22:33
Сообщение #1





Группа: Пользователи
Сообщений: 1219
Регистрация: 13.01.2008
Из: Челябинск
Пользователь №: 4704



Для каждого распределения известны параметры мю и сигма, а также объём выборки (в долях единицы).
Поиск по теме дал несколько аналогичных результатов.
Например, здесь дан вывод уравнения для нахождения абсциссы точки пересечения через решение квадратного уравнения:
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 - 19:19
Эскизы прикрепленных изображений
Прикрепленное изображение
 
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
 
Открыть тему
Ответов
nokh
сообщение 19.02.2022 - 21:21
Сообщение #2





Группа: Пользователи
Сообщений: 1219
Регистрация: 13.01.2008
Из: Челябинск
Пользователь №: 4704



> 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$
сообщение 19.02.2022 - 23:49
Сообщение #3





Группа: Пользователи
Сообщений: 902
Регистрация: 23.08.2010
Пользователь №: 22694



Цитата(nokh @ 19.02.2022 - 21:21) *
Почему-то не справляется с Zn после БК: выдаёт обычную гистограмму. Попробуйте свой экселевский код, может получится?


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

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

Прикрепленные файлы
Прикрепленный файл  Копия_Data_forum_2.xls ( 211,5 килобайт ) Кол-во скачиваний: 424
 
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 

Сообщений в этой теме
- nokh   Абсцисса пересечения двух гауссиан   17.02.2022 - 22:33
- - Диагностик   Цитата(nokh @ 18.02.2022 - 03:33) по...   18.02.2022 - 01:35
|- - nokh   Цитата(Диагностик @ 18.02.2022 - 03...   18.02.2022 - 10:54
|- - Диагностик   Цитата(nokh @ 18.02.2022 - 15:54) Я ...   18.02.2022 - 12:39
- - Игорь   Цитата(nokh @ 17.02.2022 - 22:33) Дл...   18.02.2022 - 10:59
- - 100$   Цитата(nokh @ 17.02.2022 - 22:33) Ис...   18.02.2022 - 13:20
- - nokh   Спасибо всем огромное! Сегодня утром с подачи ...   18.02.2022 - 14:24
- - Диагностик   nokh, дайте данные по гистограмме, попробую вашу с...   18.02.2022 - 14:38
|- - nokh   Цитата(Диагностик @ 18.02.2022 - 16...   18.02.2022 - 18:39
|- - 100$   Цитата(nokh @ 18.02.2022 - 18:39) пр...   18.02.2022 - 18:45
|- - nokh   Цитата(100$ @ 18.02.2022 - 20:4...   18.02.2022 - 19:15
|- - 100$   Цитата(nokh @ 18.02.2022 - 19:15) Да...   18.02.2022 - 19:57
|- - Диагностик   Цитата(nokh @ 19.02.2022 - 00:15) Да...   19.02.2022 - 01:21
|- - Диагностик   Цитата(nokh @ 19.02.2022 - 00:15) Да...   19.02.2022 - 03:36
|- - 100$   Цитата(Диагностик @ 19.02.2022 - 03...   19.02.2022 - 13:01
- - Диагностик   И на это ушло 2 ч. 15 мин.? Это не важно. Важно т...   19.02.2022 - 14:28
|- - 100$   Цитата(Диагностик @ 19.02.2022 - 14...   19.02.2022 - 16:30
- - nokh   >Диагностик Конкретно здесь неоднородность не ...   19.02.2022 - 20:57
|- - Диагностик   Цитата(nokh @ 20.02.2022 - 01:57) Дл...   20.02.2022 - 01:32
- - nokh   > 100$ Японская диаграмма прикольная. Я л...   19.02.2022 - 21:21
|- - 100$   Цитата(nokh @ 19.02.2022 - 21:21) По...   19.02.2022 - 23:49
|- - nokh   Цитата(100$ @ 20.02.2022 - 01:4...   20.02.2022 - 08:08
- - Диагностик   nokh, нужно найти аномальные значения концентрации...   20.02.2022 - 09:27
|- - nokh   Цитата(Диагностик @ 20.02.2022 - 11...   20.02.2022 - 11:59
- - Диагностик   nokh, аномальные значения левого "хвоста...   20.02.2022 - 12:07
- - Диагностик   nokh,поработал со свинцом. Оказалось что концентра...   21.02.2022 - 03:37
- - Диагностик   Цитата(Диагностик @ 21.02.2022 - 08...   24.02.2022 - 06:18
|- - Диагностик   Цитата(Диагностик @ 24.02.2022 - 11...   25.02.2022 - 06:46
|- - nokh   Цитата(Диагностик @ 24.02.2022 - 08...   25.02.2022 - 13:20
|- - Диагностик   Цитата(nokh @ 25.02.2022 - 18:20) ка...   4.03.2022 - 04:28
|- - 100$   Цитата(Диагностик @ 4.03.2022 - 04:2...   4.03.2022 - 19:18
- - Диагностик   Я не детектировал последовательно каждый выброс, а...   5.03.2022 - 02:32
- - Олег Кравец   Moderator on Коллеги, уважайте себя и собеседнико...   13.03.2022 - 08:00


Добавить ответ в эту темуОткрыть тему