Подгонка экспериментальных данных указывает на различные кумулятивные распределения с использованием R

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

Итак, мне нужно подобрать кумулятивное распределение с некоторой функцией (функция с двумя/тремя параметрами). Это кажется довольно простой задачей, но я уже некоторое время жужжу над этим.

Позвольте мне показать вам, каковы мои переменные:

    x=c(0.01,0.011482,0.013183,0.015136,0.017378,0.019953,0.022909,0.026303,0.0302,0.034674,0.039811,0.045709,0.052481,0.060256,0.069183,0.079433,0.091201,0.104713,0.120226,0.138038,0.158489,0.18197,0.20893,0.239883,0.275423,0.316228,0.363078,0.416869,0.47863,0.549541,0.630957,0.724436,0.831764,0.954993,1.096478,1.258925,1.44544,1.659587,1.905461,2.187762,2.511886,2.884031,3.311311,3.801894,4.365158,5.011872,5.754399,6.606934,7.585776,8.709636,10,11.481536,13.182567,15.135612,17.378008,19.952623,22.908677,26.30268,30.199517,34.673685,39.810717,45.708819,52.480746,60.255959,69.183097,79.432823,91.201084,104.712855,120.226443,138.038426,158.489319,181.970086,208.929613,239.883292,275.42287,316.227766,363.078055,416.869383,478.630092,549.540874,630.957344,724.43596,831.763771,954.992586,1096.478196)
    y=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00044816,0.00127554,0.00221488,0.00324858,0.00438312,0.00559138,0.00686054,0.00817179,0.00950625,0.01085188,0.0122145,0.01362578,0.01514366,0.01684314,0.01880564,0.02109756,0.0237676,0.02683182,0.03030649,0.0342276,0.03874555,0.04418374,0.05119304,0.06076553,0.07437854,0.09380666,0.12115065,0.15836926,0.20712933,0.26822017,0.34131335,0.42465413,0.51503564,0.60810697,0.69886817,0.78237651,0.85461023,0.91287236,0.95616228,0.98569093,0.99869001,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999)

Это график, где я установил ось X как журнал: Raw data

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

# sigmoid function definition
sigmoid = function(params, x) {
  params[1] / (1 + exp(-params[2] * (x - params[3])))
}

# fitting code using nonlinear least square
fitmodel <- nls(y~a/(1 + exp(-b * (x-c))), start=list(a=1,b=.5,c=25))

# get the coefficients using the coef function
params=coef(fitmodel)

# asigning to y2 sigmoid function
y2 <- sigmoid(params,x)

# plotting y2 function
plot(y2,type="l")

# plotting data points
points(y)

Это привело меня к хорошим результатам подбора (я не знаю, как это измерить). Но когда я смотрю на график функции подгонки Сигмуда, я не понимаю, почему форма S теперь происходит в диапазоне значений x от 40 до 7 (если смотреть на форму S, она должна быть в значениях x от 40 до 7). 10 до 200).

Исходные данные

Поскольку я не мог объяснить такое поведение, я подумал о том, чтобы попробовать применить уравнение Вейбулла для подгонки, но пока не могу заставить код работать.

Подводить итоги:

  1. Ты хоть представляешь, почему сигмоид дает мне такую ​​странную настройку?
  2. Знаете ли вы какое-нибудь лучшее уравнение с двумя или тремя параметрами для этого подхода?
  3. Как я мог определить качество подгонки? Что-то вроде r^2?

r sigmoid statistics curve-fitting
person numb    schedule 30.05.2017    source источник
comment
Он отображает индекс массива, потому что вы не указываете значение x. Попробуйте plot(x, y2,type="l") и points(x,y).   -  person Lyngbakr    schedule 30.05.2017
comment
@Lyngbakr Спасибо. Это решает мой первый вопрос. Я набрал plot(x,y,type="l",log="x"), чтобы лучше видеть эту S-образную кривизну. Но это лишь подтверждает, что посадка выглядит не очень.   -  person numb    schedule 30.05.2017
comment
Моей первоначальной мыслью было попробовать a + b * tanh(x/c), но это тоже дает паршивые результаты...   -  person Lyngbakr    schedule 30.05.2017

Ответы (2)


arrow_upward
1
arrow_downward 

# Data
df <- data.frame(x=c(0.01,0.011482,0.013183,0.015136,0.017378,0.019953,0.022909,0.026303,0.0302,0.034674,0.039811,0.045709,0.052481,0.060256,0.069183,0.079433,0.091201,0.104713,0.120226,0.138038,0.158489,0.18197,0.20893,0.239883,0.275423,0.316228,0.363078,0.416869,0.47863,0.549541,0.630957,0.724436,0.831764,0.954993,1.096478,1.258925,1.44544,1.659587,1.905461,2.187762,2.511886,2.884031,3.311311,3.801894,4.365158,5.011872,5.754399,6.606934,7.585776,8.709636,10,11.481536,13.182567,15.135612,17.378008,19.952623,22.908677,26.30268,30.199517,34.673685,39.810717,45.708819,52.480746,60.255959,69.183097,79.432823,91.201084,104.712855,120.226443,138.038426,158.489319,181.970086,208.929613,239.883292,275.42287,316.227766,363.078055,416.869383,478.630092,549.540874,630.957344,724.43596,831.763771,954.992586,1096.478196),
           y=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00044816,0.00127554,0.00221488,0.00324858,0.00438312,0.00559138,0.00686054,0.00817179,0.00950625,0.01085188,0.0122145,0.01362578,0.01514366,0.01684314,0.01880564,0.02109756,0.0237676,0.02683182,0.03030649,0.0342276,0.03874555,0.04418374,0.05119304,0.06076553,0.07437854,0.09380666,0.12115065,0.15836926,0.20712933,0.26822017,0.34131335,0.42465413,0.51503564,0.60810697,0.69886817,0.78237651,0.85461023,0.91287236,0.95616228,0.98569093,0.99869001,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999))

# sigmoid function definition
sigmoid = function(x, a, b, c) {
  a * exp(-b * exp(-c * x))
}

# fitting code using nonlinear least square
fitmodel <- nls(y ~ sigmoid(x, a, b, c), start=list(a=1,b=.5,c=-2), data = df)

# plotting y2 function
plot(df$x, predict(fitmodel),type="l", log = "x")

# plotting data points
points(df)

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

Я использовал функцию функции Гомперца и этот пост в блоге объясняет, почему R² должен не может использоваться с нелинейными посадками и предлагает альтернативу.

person Lyngbakr    schedule 30.05.2017
comment
Спасибо. Я не знал о существовании функции Гомперца. Это действительно хорошо подходит даже для других моих наборов данных. Что касается R ^ 2, я также где-то читал, что его рекомендуется использовать для нелинейной регрессии, и эта ссылка, которую вы отправили со стандартной ошибкой регрессии, кажется лучшим способом. - person numb; 31.05.2017
comment
Кстати, есть также пакет R специально для сигмоида функции. - person Lyngbakr; 31.05.2017
comment
Потрясающий! В этом пакете есть несколько моделей, которые я не смог найти в growthmodels. Вы мне действительно очень помогли!! - person numb; 31.05.2017
comment
Рад, что смог помочь! (И спасибо за предупреждение о моделях роста — я не знал об этом пакете.) - person Lyngbakr; 31.05.2017
comment
Что ж, после первоначального успеха у меня возникла проблема с нашим подходом из-за начальных значений предположения при попытке использовать разные наборы данных. Эта ошибка продолжала показывать Ошибка в nls(y ~ sigmoid::Gompertz(x, a, b, c), start = list(a = -1, b = 1, сингулярный градиент... Итак , после некоторых исследований я нашел этот пакет drc package Который спас мне жизнь, до сих пор. - person numb; 01.06.2017

arrow_upward
0
arrow_downward

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

Код выглядит следующим образом для указанного набора данных:

df <- data.frame(x=c(0.01,0.011482,0.013183,0.015136,0.017378,0.019953,0.022909,0.026303,0.0302,0.034674,0.039811,0.045709,0.052481,0.060256,0.069183,0.079433,0.091201,0.104713,0.120226,0.138038,0.158489,0.18197,0.20893,0.239883,0.275423,0.316228,0.363078,0.416869,0.47863,0.549541,0.630957,0.724436,0.831764,0.954993,1.096478,1.258925,1.44544,1.659587,1.905461,2.187762,2.511886,2.884031,3.311311,3.801894,4.365158,5.011872,5.754399,6.606934,7.585776,8.709636,10,11.481536,13.182567,15.135612,17.378008,19.952623,22.908677,26.30268,30.199517,34.673685,39.810717,45.708819,52.480746,60.255959,69.183097,79.432823,91.201084,104.712855,120.226443,138.038426,158.489319,181.970086,208.929613,239.883292,275.42287,316.227766,363.078055,416.869383,478.630092,549.540874,630.957344,724.43596,831.763771,954.992586,1096.478196),
       y=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.00044816,0.00127554,0.00221488,0.00324858,0.00438312,0.00559138,0.00686054,0.00817179,0.00950625,0.01085188,0.0122145,0.01362578,0.01514366,0.01684314,0.01880564,0.02109756,0.0237676,0.02683182,0.03030649,0.0342276,0.03874555,0.04418374,0.05119304,0.06076553,0.07437854,0.09380666,0.12115065,0.15836926,0.20712933,0.26822017,0.34131335,0.42465413,0.51503564,0.60810697,0.69886817,0.78237651,0.85461023,0.91287236,0.95616228,0.98569093,0.99869001,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999,0.99999999))


library(drc)
fm <- drm(y ~ x, data = df, fct = G.3()) #The Gompertz model G.3()
plot(fm)

#Gompertz Coefficients and residual standard error 

summary(fm)

Сюжет после подгонки

person numb    schedule 01.06.2017