Добрый день!
Прошу помочь с выбором метода анализа.
Данные - диаметры сферических частиц, размер выборки около 1000 шт
Задача анализа ответить на вопрос: - размеры частиц представляют собой мономодальное распределение или это смесь двух или более типов частиц.
Визуально частицы можно разделить на маленькие шарики и лепешки большего размера.
Заранее благодарю,
kmuranov
Полная версия этой страницы: Мономодальное или бимодальное распределение?
Цитата(kmuranov @ 3.05.2012 - 13:27) 
Добрый день!
Прошу помочь с выбором метода анализа.
Данные - диаметры сферических частиц, размер выборки около 1000 шт
Задача анализа ответить на вопрос: - размеры частиц представляют собой мономодальное распределение или это смесь двух или более типов частиц.
Визуально частицы можно разделить на маленькие шарики и лепешки большего размера.
Заранее благодарю,
kmuranov
Прошу помочь с выбором метода анализа.
Данные - диаметры сферических частиц, размер выборки около 1000 шт
Задача анализа ответить на вопрос: - размеры частиц представляют собой мономодальное распределение или это смесь двух или более типов частиц.
Визуально частицы можно разделить на маленькие шарики и лепешки большего размера.
Заранее благодарю,
kmuranov
Я в таких случаях ограничивался исключительно гистограммой распределения, т.к. если распределение не унимодальное, а иное, то это очевидно. Хотя потом можно подтянуть и какие-то статистические критерии. Но начинать в любом случае нужно с гистограммы распределения. Из пакетов могу посоветовать бесплатный PAST: http://folk.uio.no/ohammer/past/
Можно просто построить гистограмму по предварительно выделенному столбцу данных (Plot - Histogram) и посмотреть ядерную (kernel) плотность. В качестве числа интервалов (Bins) при 1000 наблюдений можно задать 25 или даже больше (нужно нажимать Enter после изменения). Если выявится би- или полимодальность, то можно попробовать разделить смесь распределений в разделе Model - Mixture analysis. Для смеси нормальных распределений по View numbers можно посмотреть параметры разделённых распределений (среднее и стандартное отклонение). Они вычисляются по достаточно продвинутому EM-алгоритму.
Если ваши частицы образуются в результате дробления, то распределение может быть унимодальным, но сильно скошенным - примерно логарифмически нормальным. Тогда данные можно предварительно прологарифмировать, а уже потом прогнать по описанной схеме.
Цитата(nokh @ 3.05.2012 - 15:48)
Я в таких случаях ограничивался исключительно гистограммой распределения, т.к. если распределение не унимодальное, а иное, то это очевидно. Хотя потом можно подтянуть и какие-то статистические критерии. Но начинать в любом случае нужно с гистограммы распределения. Из пакетов могу посоветовать бесплатный PAST: http://folk.uio.no/ohammer/past/
Можно просто построить гистограмму по предварительно выделенному столбцу данных (Plot - Histogram) и посмотреть ядерную (kernel) плотность. В качестве числа интервалов (Bins) при 1000 наблюдений можно задать 25 или даже больше (нужно нажимать Enter после изменения). Если выявится би- или полимодальность, то можно попробовать разделить смесь распределений в разделе Model - Mixture analysis. Для смеси нормальных распределений по View numbers можно посмотреть параметры разделённых распределений (среднее и стандартное отклонение). Они вычисляются по достаточно продвинутому EM-алгоритму.
Если ваши частицы образуются в результате дробления, то распределение может быть унимодальным, но сильно скошенным - примерно логарифмически нормальным. Тогда данные можно предварительно прологарифмировать, а уже потом прогнать по описанной схеме.
Можно просто построить гистограмму по предварительно выделенному столбцу данных (Plot - Histogram) и посмотреть ядерную (kernel) плотность. В качестве числа интервалов (Bins) при 1000 наблюдений можно задать 25 или даже больше (нужно нажимать Enter после изменения). Если выявится би- или полимодальность, то можно попробовать разделить смесь распределений в разделе Model - Mixture analysis. Для смеси нормальных распределений по View numbers можно посмотреть параметры разделённых распределений (среднее и стандартное отклонение). Они вычисляются по достаточно продвинутому EM-алгоритму.
Если ваши частицы образуются в результате дробления, то распределение может быть унимодальным, но сильно скошенным - примерно логарифмически нормальным. Тогда данные можно предварительно прологарифмировать, а уже потом прогнать по описанной схеме.
Спасибо за быстрый ответ!
С построения гистограмм я и начал - похоже, что есть три пика. Но какой метод использовать дальше - я заткнулся. Использую Statistica и Attestat. Если посоветуете методы в этих пакетах буду благодарен. А как попробовать выявить связь формы и размера? Похоже она есть!
Для уточнения задачи - частицы это олигмерный белок - альфа-кристаллин.
С уважением,
kmuranov
Цитата(nokh @ 3.05.2012 - 14:48)
Я в таких случаях ограничивался исключительно гистограммой распределения, т.к. если распределение не унимодальное, а иное, то это очевидно. Хотя потом можно подтянуть и какие-то статистические критерии. Но начинать в любом случае нужно с гистограммы распределения. Из пакетов могу посоветовать бесплатный PAST: http://folk.uio.no/ohammer/past/
Можно просто построить гистограмму по предварительно выделенному столбцу данных (Plot - Histogram) и посмотреть ядерную (kernel) плотность. В качестве числа интервалов (Bins) при 1000 наблюдений можно задать 25 или даже больше (нужно нажимать Enter после изменения). Если выявится би- или полимодальность, то можно попробовать разделить смесь распределений в разделе Model - Mixture analysis. Для смеси нормальных распределений по View numbers можно посмотреть параметры разделённых распределений (среднее и стандартное отклонение). Они вычисляются по достаточно продвинутому EM-алгоритму.
Если ваши частицы образуются в результате дробления, то распределение может быть унимодальным, но сильно скошенным - примерно логарифмически нормальным. Тогда данные можно предварительно прологарифмировать, а уже потом прогнать по описанной схеме.
Можно просто построить гистограмму по предварительно выделенному столбцу данных (Plot - Histogram) и посмотреть ядерную (kernel) плотность. В качестве числа интервалов (Bins) при 1000 наблюдений можно задать 25 или даже больше (нужно нажимать Enter после изменения). Если выявится би- или полимодальность, то можно попробовать разделить смесь распределений в разделе Model - Mixture analysis. Для смеси нормальных распределений по View numbers можно посмотреть параметры разделённых распределений (среднее и стандартное отклонение). Они вычисляются по достаточно продвинутому EM-алгоритму.
Если ваши частицы образуются в результате дробления, то распределение может быть унимодальным, но сильно скошенным - примерно логарифмически нормальным. Тогда данные можно предварительно прологарифмировать, а уже потом прогнать по описанной схеме.
По моему еще полезно не одну реализацию ядерного сглаживания получить, а бутстреп учинить
Имея несколько сот кривых сразу станет ясно сколько пиков имеют право на жизнь.Топикстартер, давайте Ваш датасет в тред. Посмотрим что в нем есть
Цитата(p2004r @ 3.05.2012 - 22:42)
По моему еще полезно не одну реализацию ядерного сглаживания получить, а бутстреп учинить Имея несколько сот кривых сразу станет ясно сколько пиков имеют право на жизнь.
Топикстартер, давайте Ваш датасет в тред. Посмотрим что в нем есть
Имея несколько сот кривых сразу станет ясно сколько пиков имеют право на жизнь.Топикстартер, давайте Ваш датасет в тред. Посмотрим что в нем есть
Посылаю в txt.
Спасибо,
kmuranov
Цитата(kmuranov @ 4.05.2012 - 11:39)
Посылаю в txt.
Спасибо,
kmuranov
Спасибо,
kmuranov
Код
## Собственно эксперимент по накомлению 10000 перевыборок с возвращением и построением плотности распределения
xxxx<-replicate(10000,
density(sample(data$size,
size=length(data[,1]),
replace=TRUE),
from=8,
to=28)$y)
line.means<-rowMeans(xxxx) # вычисляем среднее для каждой из 512 точек оценки плотности распределения
## вычисляем положение точек в которых оценивалась плотность распределения
x <- density(sample(data$size, size=length(data[,1]), replace=TRUE), from=8, to=28)$x
## рисуем
plot(x, line.means, pch=".")
## считаем границы доверительного 95% интервала
line.high<-sapply(1:512, function (i) {sort(xxxx[i,], decreasing = TRUE)[0.025*10000]})
line.low<-sapply(1:512, function (i) {sort(xxxx[i,], decreasing = TRUE)[0.975*10000]})
lines(x, line.low, col="green")
lines(x, line.high, col="red")
xxxx<-replicate(10000,
density(sample(data$size,
size=length(data[,1]),
replace=TRUE),
from=8,
to=28)$y)
line.means<-rowMeans(xxxx) # вычисляем среднее для каждой из 512 точек оценки плотности распределения
## вычисляем положение точек в которых оценивалась плотность распределения
x <- density(sample(data$size, size=length(data[,1]), replace=TRUE), from=8, to=28)$x
## рисуем
plot(x, line.means, pch=".")
## считаем границы доверительного 95% интервала
line.high<-sapply(1:512, function (i) {sort(xxxx[i,], decreasing = TRUE)[0.025*10000]})
line.low<-sapply(1:512, function (i) {sort(xxxx[i,], decreasing = TRUE)[0.975*10000]})
lines(x, line.low, col="green")
lines(x, line.high, col="red")
ну и они полностью разделены
Цитата(kmuranov @ 3.05.2012 - 17:24)
... Использую Statistica и Attestat. Если посоветуете методы в этих пакетах буду благодарен.
Я их тоже использую, но для других задач, там нет разделения смесей распределений. Другие методы здесь вряд ли уместны, но к сожалению доказать не графически, а статистически, что мод именно 3 (или 5-6) я тоже не знаю как.
Цитата(kmuranov @ 3.05.2012 - 17:24)
А как попробовать выявить связь формы и размера? Похоже она есть!
Выявлять, собственно, нечего, т.к. размеры разных форм не трансгрессируют. Как уже сказал p2004r они полностью разделены. Можно точно указать до каких размеров идёт шар и с какого начинается тороид. Если нужен статистический критерий, то подойдёт даже критерий знаков.
>p2004r
PAST для трёх распределений выдаёт параметры указанные на картинке. Хотя, похоже, первый кластер тоже неоднородный и состоит из 2 или даже скорее 3 подгрупп. Нашёл пакет для R, который должен делать аналогичное разделение и даже больше: http://www.math.mcmaster.ca/peter/mix/mix.html . Если у вас есть время, интересно было бы сопоставить результаты пакетов.
Самое забавное во всей этой истории- это то , что фраза
представляет собой ответ на вопрос
Если "лепешки" не являются агломератом (несколькими слипшимися "шариками"), то и без возни с распределениями понятно, что в детерминированных науках (физика, химия) одинаковые причины не могут приводить к разным последствиям (форме частиц). Следовательно, существуют, как минимум, две причины, влияющие на пространственную форму частиц. Эти две причины формируют две статистические совокупности. Можно, конечно, на радостях объявить их смесью распределений, только что это дает?
Цитата
Визуально частицы можно разделить на маленькие шарики и лепешки большего размера.
представляет собой ответ на вопрос
Цитата
...размеры частиц представляют собой мономодальное распределение или это смесь двух или более типов частиц.
Если "лепешки" не являются агломератом (несколькими слипшимися "шариками"), то и без возни с распределениями понятно, что в детерминированных науках (физика, химия) одинаковые причины не могут приводить к разным последствиям (форме частиц). Следовательно, существуют, как минимум, две причины, влияющие на пространственную форму частиц. Эти две причины формируют две статистические совокупности. Можно, конечно, на радостях объявить их смесью распределений, только что это дает?
Цитата(nokh @ 4.05.2012 - 20:33)
>p2004r
PAST для трёх распределений выдаёт параметры указанные на картинке. Хотя, похоже, первый кластер тоже неоднородный и состоит из 2 или даже скорее 3 подгрупп. Нашёл пакет для R, который должен делать аналогичное разделение и даже больше: http://www.math.mcmaster.ca/peter/mix/mix.html . Если у вас есть время, интересно было бы сопоставить результаты пакетов.
PAST для трёх распределений выдаёт параметры указанные на картинке. Хотя, похоже, первый кластер тоже неоднородный и состоит из 2 или даже скорее 3 подгрупп. Нашёл пакет для R, который должен делать аналогичное разделение и даже больше: http://www.math.mcmaster.ca/peter/mix/mix.html . Если у вас есть время, интересно было бы сопоставить результаты пакетов.
Увы только в понедельник смогу продолжить
Мне представляется что "лепешка" в зависимости от того под каким углом видна даст различный размер в проекции. И даже если все "лепешки" равны получится некое распределение. Что то типа задачи расчета длинны нити в клубке по наблюдаемому сечению клубка.
В случае шарообразной формы проецирование дает одинаковую картинку.
Цитата(p2004r @ 5.05.2012 - 01:44)
Увы только в понедельник смогу продолжить
Мне представляется что "лепешка" в зависимости от того под каким углом видна даст различный размер в проекции. И даже если все "лепешки" равны получится некое распределение. Что то типа задачи расчета длинны нити в клубке по наблюдаемому сечению клубка.
В случае шарообразной формы проецирование дает одинаковую картинку.
Мне представляется что "лепешка" в зависимости от того под каким углом видна даст различный размер в проекции. И даже если все "лепешки" равны получится некое распределение. Что то типа задачи расчета длинны нити в клубке по наблюдаемому сечению клубка.
В случае шарообразной формы проецирование дает одинаковую картинку.
Полагаю, что если диаметр тороида измеряется в разных проекциях, то это приведёт к унимодальному распределению. А раз в распределении лепёшек два пика, значит это разные классы лепёшек.
Попробовал пакет сам. В результатах есть различия, но я пока не разбирался с чем они связаны: с несколько иной группировкой (здесь разбивал на 17 классов, а получилось 18), с различиями алгоритма или настройкой типа распределения. Пока просто сделал. Как смог
!!!Переименовал выложенный файл в belok.txt
> belok<-read.table("data/belok.txt", h=T)
> attach(belok)
> library(mixdist)
> szgr<-mixgroup(belok[,1],breaks=c(0,seq(11.8,24.6,0.8),26))
> szgr
X count
1 11.8 180
2 12.6 147
3 13.4 156
4 14.2 122
5 15.0 5
6 15.8 60
7 16.6 67
8 17.4 60
9 18.2 53
10 19.0 67
11 19.8 23
12 20.6 28
13 21.4 31
14 22.2 38
15 23.0 31
16 23.8 34
17 24.6 29
18 Inf 12
> plot(szgr)
Смотрим на получившуюся гистограмму и задаём примерные центры кластеров.
> fitclaster<-mix(szgr,mixparam(c(13,18,22.5),.5),"gamma",mixconstr(consigma="CCV"))
> summary(fitclaster)
Parameters:
pi mu sigma
1 0.5332 12.45 1.029
2 0.2832 17.14 1.416
3 0.1836 22.22 1.835
Standard Errors:
pi.se mu.se sigma.se
1 0.01533 0.04778 0.03441
2 0.01511 0.10971 NA
3 0.01322 0.16570 NA
Analysis of Variance Table
Df Chisq Pr(>Chisq)
Residuals 11 107.99 < 2.2e-16 ***
Получается, что модель с тремя кластерами плохо приближает реальность. Полагаю, что это связано в первую очередь со смешанным характером первого кластера, который состоит примерно из 3 групп близких по размеру объектов.
> plot(fitclaster)
Понял смысл не всех параметров команд, которые срисовал. Ну да бог с ними. Важнее результат и картинка. Подскажите пожалуйста, какие настройки нужно добавить в последнюю команду plot, чтобы по оси Y влез весь рисунок, а также чтобы между большими метками обеих шкал сделать коротенькие засечки маленьких меток.
Цитата(nokh @ 8.05.2012 - 12:33)
Понял смысл не всех параметров команд, которые срисовал. Ну да бог с ними. Важнее результат и картинка. Подскажите пожалуйста, какие настройки нужно добавить в последнюю команду plot, чтобы по оси Y влез весь рисунок, а также чтобы между большими метками обеих шкал сделать коротенькие засечки маленьких меток.
у plot есть параметр задающий лимиты по осям (?plot.default ---> ylim=c(y1,y2) ), но в данном случае работает метод из пакета и боюсь результат будет деструктивным.
тоже касается меток на осях... надо посмотреть в процедуру рисования определенную в пакете, там явно параметры берутся из графика гистограммы а не графика плотности... обычно рассчитывают максимумы и рисуют пустой график с нужными осями, потом на него в нужном порядке выводят все остальное.
Дорогие участники обсуждения, спасибо всем!
То, что хотел, то обрёл!
С уважением,
kmuranov
То, что хотел, то обрёл!
С уважением,
kmuranov
Цитата(nokh @ 8.05.2012 - 12:33)
Полагаю, что если диаметр тороида измеряется в разных проекциях, то это приведёт к унимодальному распределению. А раз в распределении лепёшек два пика, значит это разные классы лепёшек.
случай когда у нас "палочки" выглядит так
Код
densityplot(cos(runif(10000, # число случайных "палочек"
min=0, # угол нормали палочки плашмя
max=pi/2))) # угол палочки стоящей вертикально
min=0, # угол нормали палочки плашмя
max=pi/2))) # угол палочки стоящей вертикально
случай когда у нас окружность, это эквивалентно -- две перпендикулярных друг другу "палочки" с (как следует из перпендикуляности) полностью независимыми углами к нормали. Описанная вокруг проекции окружность будет ограничена максимальной проекцией.
Код
densityplot(pmax(cos(runif(10000,min=0, max=pi/2)),
cos(runif(10000,min=0, max=pi/2))))
cos(runif(10000,min=0, max=pi/2))))
ну и неожиданно когда у нас сплюснутая окружность в отношении 1 к 2 короткий диаметр к длинному
Код
densityplot(pmax(1*cos(runif(10000,min=0, max=pi/2)),
2*cos(runif(10000,min=0, max=pi/2))))
2*cos(runif(10000,min=0, max=pi/2))))
Все построено из предположения нулевой толщины объекта.
Круто! Способностями к математическому моделированию не обладаю, видимо остаётся только поверить, что распределение тороидов бимодальное Как человек с цитогенетическим прошлым имел дело с тороидами и их размерами в виде эритроцитов и эритробластов (микроядерный тест). Там распределение унимодальное, но объект специфический: при высыхании препарата клетки распределяются преимущественно в плоскости стекла. Возможно, если (1) силы, разворачивающие объекты в пространстве, отсутствуют или слабы, а (2) измерения осуществляются автоматом, то распределение получится как у вас.
Как человек с цитогенетическим прошлым имел дело с тороидами и их размерами в виде эритроцитов и эритробластов (микроядерный тест). Там распределение унимодальное, но объект специфический: при высыхании препарата клетки распределяются преимущественно в плоскости стекла. Возможно, если (1) силы, разворачивающие объекты в пространстве, отсутствуют или слабы, а (2) измерения осуществляются автоматом, то распределение получится как у вас.
Цитата(nokh @ 21.05.2012 - 22:03)
Способностями к математическому моделированию не обладаю, видимо остаётся только поверить, что распределение тороидов бимодальное
Хотите верьте, а хотите не верьте, но у Вас, судя по исходным данным, распределение вырожденное. А у вырожденного распределения никаких модальностей быть не может по определению.
Цитата(nokh @ 21.05.2012 - 20:03)
Круто! Способностями к математическому моделированию не обладаю, видимо остаётся только поверить, что распределение тороидов бимодальное Как человек с цитогенетическим прошлым имел дело с тороидами и их размерами в виде эритроцитов и эритробластов (микроядерный тест). Там распределение унимодальное, но объект специфический: при высыхании препарата клетки распределяются преимущественно в плоскости стекла. Возможно, если (1) силы, разворачивающие объекты в пространстве, отсутствуют или слабы, а (2) измерения осуществляются автоматом, то распределение получится как у вас.
Как человек с цитогенетическим прошлым имел дело с тороидами и их размерами в виде эритроцитов и эритробластов (микроядерный тест). Там распределение унимодальное, но объект специфический: при высыхании препарата клетки распределяются преимущественно в плоскости стекла. Возможно, если (1) силы, разворачивающие объекты в пространстве, отсутствуют или слабы, а (2) измерения осуществляются автоматом, то распределение получится как у вас.Конечно у меня грубая модель (в ней только два диаметра и нулевая толщина). Окружность при точном описании должна дать дельта функцию например. Вот например как качественно распределение проекции слегка неточной сферы должно выглядеть.
Код
normal1 <-runif(1000,min=0, max=pi/2);
normal2 <-runif(1000,min=0, max=pi/2);
densityplot(pmax(cos(normal1), 0.8*cos(normal1-pi/2), 1.3*cos(normal2), 0.8*cos(normal2-pi/2) ))
normal2 <-runif(1000,min=0, max=pi/2);
densityplot(pmax(cos(normal1), 0.8*cos(normal1-pi/2), 1.3*cos(normal2), 0.8*cos(normal2-pi/2) ))
Свободное от ограничений распределение наверное при поточной флоуметрии должно получаться.
Ремарка к графикам. R по умолчанию делает графики далёкие от того, как они должны выглядеть в публикациях. Я видел коды, которые позволяют это исправить - они громоздкие. Полагаю, что для большинства людей (кроме может самых ярых фанатов TEX) редактировать графику удобнее напрямую. Я нашёл замечательный векторный редактор для научной графики от автора эконометрического пакета matrixer. Нахожу очень удобным забирать графику в буфер в пакетах со слабыми графическими возможностями (типа PAST) или сложным редактором (типа R) и доводить до ума в нём.
http://tpx.sourceforge.net/
http://tpx.sourceforge.net/
Цитата(nokh @ 22.05.2012 - 16:27)
Ремарка к графикам. R по умолчанию делает графики далёкие от того, как они должны выглядеть в публикациях. Я видел коды, которые позволяют это исправить - они громоздкие. Полагаю, что для большинства людей (кроме может самых ярых фанатов TEX) редактировать графику удобнее напрямую. Я нашёл замечательный векторный редактор для научной графики от автора эконометрического пакета matrixer. Нахожу очень удобным забирать графику в буфер в пакетах со слабыми графическими возможностями (типа PAST) или сложным редактором (типа R) и доводить до ума в нём.
http://tpx.sourceforge.net/
http://tpx.sourceforge.net/
Ну статистика фор виндовс (или там эксель) делает графики вообще далекие от печати , не то что от публикаций
и наука как то движется В R у графики просто гора параметров настройки можно сделать практически все из самого R.
Я бы предостерёг от редактирования графиков после их получения. Причина проста: что то поменяется в данных (ну там опечатка найдется) и получить повторно график можно будет только кропотливым повторным редактированием. А если все сделать изнутри R (ну или связкой R + какой pgf), то повторить расчет на исправленных данных можно столько раз сколько понадобится.
Хотя конечно есть случаи когда очень хочется что то просто подвинуть "визуально". Или как модно лепить внутрь рисунка всякие подрисунки и диаграммы...
Для просмотра полной версии этой страницы, пожалуйста, пройдите по ссылке.