Определение FWHM из распределения в R

У меня есть набор данных:

 x = c(30, 23, 22, 31, 18, 16, 19, 16, 15, 23, 21, 17, 17, 27, 24, 22, 27, 32, 23, 21, 14, 19, 22, 23, 22, 19, 13, 22, 33, 25, 24, 26, 24, 24, 13, 21, 23, 23, 31, 23, 25, 20, 24, 23, 24, 16, 19, 33, 36, 29, 29, 26, 22, 23, 23, 23, 20, 27, 32, 30, 34, 37, 28, 24, 29, 32, 35, 32, 25, 28, 27, 32, 33, 25, 25, 30, 25, 25, 26, 36, 30, 37, 27, 25, 26, 34, 38, 29, 30, 32, 31, 33, 33, 31, 37, 36, 42, 32, 37, 34, 37, 38, 31, 25, 27, 27, 32, 33, 39, 38, 36, 32, 35, 33, 36, 40, 39, 32, 32, 34, 35, 38, 43, 43, 39, 30, 33, 39, 39, 46, 51, 40, 30, 37, 48, 52, 47, 58, 41, 47, 51, 55, 56, 46, 34, 47, 45, 48, 52, 66, 56, 49, 66, 71, 63, 47, 41, 50, 56, 58, 56, 56, 61, 49, 57, 58, 52, 65, 70, 75, 59, 55, 44, 48, 49, 49, 55, 61, 58, 56, 56, 62, 64, 58, 59, 64, 64, 61, 56, 62, 64, 78, 86, 80, 75, 65, 75, 66, 72, 85, 65, 70, 63, 58, 73, 83, 89, 78, 75, 79, 84, 92, 90, 80, 76, 82, 83, 82, 82, 93, 93, 79, 92, 97, 84, 85, 81, 90, 99, 106, 100, 96, 97, 101, 128, 120, 121, 109, 125, 118, 121, 118, 123, 108, 115, 123, 123, 117, 113, 124, 128, 142, 134, 120, 116, 112, 129, 149, 150, 134, 134, 133, 139, 165, 167, 173, 158, 170, 188, 186, 191, 172, 163, 176, 182, 188, 205, 203, 219, 211, 229, 242, 245, 273, 270, 285, 292, 329, 355, 358, 362, 403, 429, 480, 516, 512, 525, 544, 575, 595, 668, 622, 612, 627, 649, 620, 576, 608, 576, 545, 471, 435, 422, 416, 405, 372, 338, 299, 285, 279, 274, 251, 241, 213, 201, 197, 203, 197, 196, 189, 184, 166, 165, 161, 167, 160, 157, 139, 131, 141, 152, 143, 144, 140, 136, 148, 123, 114, 113, 109, 122, 134, 120, 103, 117, 134, 117, 106, 114, 112, 104, 86, 94, 103, 108, 98, 109, 98, 100, 108, 114, 92, 78, 83, 111, 98, 78, 80, 80, 68, 62, 76, 74, 89, 78, 85, 86, 86, 76, 71, 72, 72, 64, 76, 77, 77, 89, 76, 65, 61, 66, 68, 72, 75, 72, 67, 67, 69, 75, 65, 63, 75, 68, 65, 59, 68, 61, 60, 63, 63, 58, 63, 59, 49, 68, 55, 60, 67, 65, 69, 68, 53, 59, 64, 45, 43, 42, 48, 46, 50, 52, 41, 38, 44, 38, 51, 50, 51, 41, 40, 41, 41, 34, 41, 32, 35, 39, 52, 46, 38, 37, 39, 36, 36, 34, 41, 40, 38, 38, 47, 45, 46, 36, 40, 34, 32, 39, 41, 47, 38, 38, 33, 44, 37, 38, 30, 34, 30, 40, 43, 41, 31, 27, 39, 34, 31, 29, 29, 25, 38, 38, 33, 42, 45, 46, 42, 37, 40, 35, 50, 34, 29, 25, 30, 36, 35, 36, 35, 24, 22, 29, 29, 32, 32, 25, 32, 30, 28, 23, 28, 34, 31, 28, 30, 27, 27, 20, 25, 32, 32, 41, 28, 19, 22, 23, 20, 25, 31, 27, 24, 26, 21, 20, 25, 33, 31, 44, 31, 31, 22, 29, 29, 32, 20, 24, 26, 27, 28, 24, 16, 19, 24, 23, 28, 27, 22, 24, 18, 19, 19, 21, 26, 26, 25, 28, 28, 32, 32, 26, 23, 31, 27, 20, 18, 29, 25, 15, 23, 28, 29)

Я знаю, что данные соответствуют следующему распределению, квадратичному фону с пиком Коши:

a + bx + cx^2 + dx^3 + ex^4 + fx^5 + gx^6 + (A/pi)(gamma/((x-pos)^2 + gamma^2))

Я пытаюсь определить оптимальные параметры, соответствующие этим данным, чтобы найти параметр «гамма». Это покажет мне FWHM распределения. Я знаю только параметр pos. Это расположение пика.

Вот график этого распределения:

График значений данных

Я добился этого на python с помощью пакета lmfit, но пока не нашел способа сделать это в R. Я считаю, что для этого требуется пакет fitdistrplus, но мне не повезло. Любые идеи?


person Community    schedule 23.04.2018    source источник
comment
вы можете посмотреть здесь: stat.ethz .ch / R-manual / R-devel / library / stats / html / nls.html. Думаю, ты мог бы сделать что-нибудь.   -  person denis    schedule 23.04.2018


Ответы (1)


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

> d <- density(na.omit(x),n=1e4)
> xmax <- d$x[d$y==max(d$y)]
> x1 <- d$x[d$x < xmax][which.min(abs(d$y[d$x < xmax]-max(d$y)/2))]
> x2 <- d$x[d$x > xmax][which.min(abs(d$y[d$x > xmax]-max(d$y)/2))]
> points(c(x1, x2), c(d$y[d$x==x1], d$y[d$x==x2]), col="red")
> (FWHM <- x2-x1)
[1] 43.37897

fwhm-density

person mysteRious    schedule 23.04.2018
comment
Я не думаю, что это сработает с моими данными. Этот код работает, находя половину высоты распределения. Поскольку мои данные находятся на квадратичном фоне, они влияют на то, что будет рассматриваться как середина пика. Значение гаммы в моем распределении должно отличаться от ширины пика на половине высоты общего распределения. Это должна быть ширина козырька на половине высоты козырька. - person ; 23.04.2018