Как подогнать к этим данным кривую Гаусса?

Я новичок в R и пытаюсь получить кривую, соответствующую этим данным разброса, которые дают мне кривую Гаусса. Я был бы очень признателен за любую помощь. Данные:

library(tidyverse)
MK20 <- tribble(~X.Intensity,    ~Average,
             0.400,  0.0000000,
             0.463,  0.0000000,
             0.536,  0.000000,
             0.621,  0.0000000,
             0.719,  0.0000000,
             0.833,  0.0000000,
             0.965,  0.0000000,
             1.120,  0.0000000,
             1.290,  0.0000000,
             1.500,  0.0000000,
             1.740,  0.0000000,
             2.010,  0.0000000,
             2.330,  0.0000000,
             2.700,  0.0000000,
             3.120,  0.0000000,
             3.620,  0.0000000,
             4.190,  0.0000000,
             4.850,  0.0000000,
             5.610,  0.0000000,
             6.500,  0.0000000,
             7.530,  0.0000000,
             8.720,  0.0000000,
             10.100,  0.0000000,
             11.700,  0.0000000,
             13.500,  0.0000000,
             15.700,  0.0000000,
             18.200,  0.0000000,
             21.000,  0.0000000,
             24.400,  0.0000000,
             28.200,  0.0000000,
             32.700,  0.0000000,
             37.800,  0.0000000,
             43.800,  0.7023333,
             50.700,  3.3700000,
             58.800,  7.3933333,
             68.100, 11.4666667,
             78.800, 14.3666667,
             91.300, 15.4000000,
             106.000, 14.5000000,
             122.000, 12.0000000,
             142.000,  8.6433333,
             164.000,  5.2200000,
             190.000,  2.4500000,
             220.000,  0.7580000,
             255.000,  0.1306667,
             295.000,  0.0000000,
             342.000,  0.0000000,
             396.000,  0.0000000,
             459.000,  0.0000000,
             531.000,  0.0000000,
             615.000,  0.0000000,
             712.000,  0.0000000,
             825.000,  0.0000000,
             955.000,  0.0000000,
             1110.000,  0.0000000,
             1280.000,  0.0000000,
             1480.000,  0.0000000,
             1720.000,  0.0000000,
             1990.000,  0.0000000,
             2300.000,  0.0000000,
             2670.000,  0.0000000,
             3090.000,  0.0000000,
             3580.000,  0.0000000,
             4150.000,  0.0000000,
             4800.000,  0.0000000,
             5560.000,  0.0000000,
             6440.000,  0.0000000,
             7460.000,  0.0000000,
             8630.000,  0.0000000)

Код, который я использую для построения:

plot(log10(MK20$X.Intensity), MK20$Average, col=1, xlim=c(-0.5,4), ylim=c(0,20), xlab="Log(Average diameter)", ylab="Intensity", xaxt='n')

Я использую функцию minor.tick.axis для добавления второстепенных делений на логарифмическую ось x. Я хочу добавить к этим данным кривую Гаусса (которая лучше всего подходит). Я попытался добавить type='l' на график, но кривая не была гладкой, и мне не нужна кривая, которая обязательно касается каждой точки данных, а только та, которая подходит лучше всего.

Извините, если решение простое, но я не смог его понять.


r gaussian best-fit-curve
person Manps    schedule 10.12.2019    source источник
comment
Кажется, это дубликат связанного целевого вопроса. Я рекомендую взглянуть на ответы на вопрос о цели обмана. Ни один из ответов, опубликованных ниже, не показывает, как соответствовать нормальному распределению данных.   -  person Maurits Evers    schedule 10.12.2019

Ответы (2)


arrow_upward
0
arrow_downward

Кривая, которая касается каждой точки, наверняка лучше всего соответствует вашим данным. :)

Кроме того, вы можете попытаться включить сглаженную кривую, например.

plot(log10(MK20$X.Intensity), MK20$Average, col=1, xlim=c(-0.5,4), ylim=c(0,20), 
     xlab="Log(Average diameter)", ylab="Intensity", xaxt='n', type='n')
lines(lowess(MK20$Average ~ log10(MK20$X.Intensity), f=0.3))

Вы можете изменить параметр f= между (0 и 1), чтобы изменить уровень сглаживания.

Вот как выглядит результат с f = 0,3.

введите здесь описание изображения

person Otto Kässi    schedule 10.12.2019

arrow_upward
0
arrow_downward

В этом случае мы не можем использовать обычный подход fitdistr для подбора нормального распределения, потому что у нас нет исходных данных. Похоже, что столбец «Среднее» представляет собой некоторую оценку плотности. Если бы это был pdf, то он должен быть интегрирован в 1, но это не так.

f <- approxfun(x = log10(MK20$X.Intensity), y= MK20$Average)
integrate(f, lower = log10(0.4), upper = log10(8630))

#6.142134 with absolute error < 0.00043

Таким образом, мы могли бы превратить это в PDF, уменьшив его примерно на 6,14, а затем попытаться найти среднее значение и стандартное отклонение, соответствующие этому PDF.

Вот первая попытка простой подгонки по Гауссу. Сначала я выбрал среднее значение 2 (посмотрев, где плотность была наибольшей), коэффициент масштабирования k = 6,14 (значение интеграла), а затем поиграл с sd, пока не было разумного соответствия.

m=2
s=0.15
k=6.14
x_seq = seq(1,3,length.out = 100)
df <- tibble(x_seq = x_seq, dens = dnorm(x_seq, m, s))


MK20 %>% 
  mutate(log_intensity = log10(X.Intensity)) %>% 
  ggplot(aes(log_intensity, Average/k)) +
  geom_point() +
  geom_line(data = df, aes(x_seq, dens))

введите здесь описание изображения

Затем я использовал optimx, чтобы подобрать 3 параметра (k = коэффициент масштабирования, m = среднее значение, s = стандартное отклонение), минимизировав сумму квадратов между подгонкой и данными.

Целевая функция (сумма квадратов разностей между подгонкой и данными)

f <- function(x) {
  k = x[1]
  m = x[2]
  s = x[3]
  MK20 %>% 
  mutate(log_intensity = log10(X.Intensity)) %>%
  mutate(fit = dnorm(log_intensity, m, s)) %>% 
  summarise(sum((fit - Average/k)^2)) %>% pull
}

Используйте optimx для поиска параметров (минимизация суммы квадратов). Начальные значения параметров берутся из подгонки на глаз.

library(optimx)    
optimx(par = c(6.14, 2, 0.15), fn = f )

#k = 6.294696 m = 1.971488 s= 0.1583936

Давайте перестроим с подходящими параметрами

# points for a gaussian
x_seq = seq(1,3,length.out = 100) 
df <- tibble(x_seq = x_seq, dens = dnorm(x_seq, m, s))


MK20 %>% 
  mutate(log_intensity = log10(X.Intensity)) %>% 
  ggplot(aes(log_intensity, Average/k)) +
  geom_point() +
  geom_line(data = df, aes(x_seq, dens))

введите здесь описание изображения

person Tony Ladson    schedule 10.12.2019
comment
Это не подходит. Вы строите нормальное распределение для выбранного набора параметров (среднее значение и дисперсия). Подгонка — это процедура вывода для получения оценок параметров, зависящих от данных. Можете ли вы пересмотреть свой ответ, чтобы показать фактическое соответствие? - person Maurits Evers; 10.12.2019
comment
@MauritsEvers мои параметры зависят от данных - я вручную выбрал параметры, чтобы подобрать их на глаз. Однако вы правы, возможно и уместно сделать лучше. Теперь добавили числовую процедуру подгонки. - person Tony Ladson; 10.12.2019
comment
Выбор параметров вручную для подгонки на глаз — это не то, как сообщество статистики обычно интерпретирует вопрос о том, как подогнать модель к данным. На самом деле, это может быть очень опасно, так как глаз может плохо судить (см., например, Anscombe квартет). Не поймите неправильно, но это не очень хороший ответ. Во-первых, это свидетельствует о плохой статистической практике. Кроме того, в своем редактировании вы не даете никаких объяснений (параметров, кода, алгоритма и т. д.), демонстрирующих использование optimx. [...] - person Maurits Evers; 11.12.2019
comment
[продолжение] Канонический способ подгонки нормального распределения к данным состоит в использовании MASS:: fitdistr или fitdistrplus::fitdist, и в этом случае ответ на вопрос ОП становится однострочным (см. Цель обмана). Stack Overflow стремится создать всеобъемлющий каталог хороших вопросов и ответов, которые будут полезны людям, столкнувшимся с похожими проблемами. По причинам, указанным выше, я не думаю, что ваш ответ отвечает этим требованиям. - person Maurits Evers; 11.12.2019
comment
@MauritsEvers, возможно, я неправильно истолковываю вопрос, но я не думаю, что подход fitdistr здесь сработает, потому что у нас нет исходных данных. У нас есть некоторая оценка плотности. В любом случае, я отредактировал свой ответ, чтобы объяснить, откуда я. Я также последовал вашему совету, чтобы объяснить, что делает optimx. - person Tony Ladson; 13.12.2019
comment
Я согласен с вами в том, что fitdistr здесь не сработает. Я удалил свой отрицательный голос. Я все еще думаю, что это не очень хороший ответ (помимо того факта, что это дублирующий вопрос). В нем используется нестандартная терминология (более канонический и традиционный статистический подход определил бы функцию правдоподобия dnorm(x, mean, sd), а затем использовал бы mle для получения оценок MLE для mean и sd вместо того, чтобы по существу заново изобретать колесо и минимизировать сумму квадратов). между подгонкой и данными). Моя самая большая проблема по-прежнему связана со всем бизнесом на глаз. Учит очень плохим практикам. - person Maurits Evers; 16.12.2019