R Нелинейный метод наименьших квадратов (nls) Подгонка модели

Я пытаюсь подогнать информацию из функции G моих данных к следующему математическому режиму: y = A / ((1 + (B ^ 2) * (x ^ 2)) ^ ((C + 1) / 2 )). Форму этого графика можно увидеть здесь:

http://www.wolframalpha.com/input/?i=y%20=%201/%20%28%281%20%2b%20%282%5E2%29%2a%28x%5E2%29%29%5E%28%282%2b1%29/2%29%29

Вот простой пример того, что я делал:

data(simdat)

library(spatstat)

simdat.Gest <- Gest(simdat) #Gest is a function within spatstat (explained below)

Gvalues <- simdat.Gest$rs

Rvalues <- simdat.Gest$r

GvsR_dataframe <- data.frame(R = Rvalues, G = rev(Gvalues))

themodel <- nls(rev(Gvalues) ~ (1 / (1 + (B^2)*(R^2))^((C+1)/2)), data = GvsR_dataframe, start = list(B=0.1, C=0.1), trace = FALSE)

«Gest» - это функция из библиотеки spatstat. Это функция G или функция ближайшего соседа, которая отображает расстояние между частицами на независимой оси в зависимости от вероятности обнаружения ближайшей соседней частицы на зависимой оси. Таким образом, он начинается при y = 0 и достигает точки насыщения при y = 1.

Если вы построите simdat.Gest, вы заметите, что кривая имеет s-образную форму, что означает, что она начинается с y = 0 и заканчивается при y = 1. По этой причине я уважал векторные G-значения, которые являются зависимыми переменные. Таким образом, информация находится в правильной ориентации, чтобы соответствовать указанной выше модели.

Вы также можете заметить, что я автоматически установил A = 1. Это потому, что G (r) всегда насыщается на 1, поэтому я не стал сохранять его в формуле.

Моя проблема в том, что я продолжаю получать ошибки. В приведенном выше примере я получаю эту ошибку:

Error in nls(rev(Gvalues) ~ (1/(1 + (B^2) * (R^2))^((C + 1)/2)), data = GvsR_dataframe,  : 
  singular gradient

Я также получаю эту ошибку:

Error in nls(Gvalues1 ~ (1/(1 + (B^2) * (x^2))^((C + 1)/2)), data = G_r_dataframe,  : 
  step factor 0.000488281 reduced below 'minFactor' of 0.000976562

Я понятия не имею, откуда взялась первая ошибка. Второе, однако, я считаю, произошло потому, что я не выбрал подходящие начальные значения для B и C.

Я надеялся, что кто-нибудь поможет мне выяснить, откуда произошла первая ошибка. Кроме того, как наиболее эффективно выбрать начальные значения, чтобы избежать второй ошибки?

Спасибо!


r model plot model-fitting nls
person MikeZ    schedule 27.06.2012    source источник
comment
Когда вы набираете simdat.Gest <- Gest(simdat), вы говорите нам, что у вас есть функция с именем Gest. Но вы не дали его нам. Я не думаю, что было бы так сложно создать тестовый набор данных с использованием rnorm, хотя в идеале нам должны были бы дать, скажем, первые двадцать строк, но нам действительно нужно знать, что делает Gest.   -  person IRTFM    schedule 28.06.2012
comment
Вы также используете в этой формуле другое имя ('Gvalues'), чем используется в data.frame ('G')   -  person IRTFM    schedule 28.06.2012
comment
Это была моя ошибка, я сейчас отредактировал пост. Gest - это функция из библиотеки spatstat. Gest - это функция ближайшего соседа, которая отображает расстояние между частицами на независимой оси в зависимости от вероятности обнаружения ближайшей соседней частицы на зависимой оси. Таким образом, он начинается при y = 0 и достигает точки насыщения при y = 1.   -  person MikeZ    schedule 28.06.2012
comment
Кроме того, я пробовал использовать nls.lm, и это тоже меня огорчило.   -  person MikeZ    schedule 05.07.2012

Ответы (1)


arrow_upward
3
arrow_downward

Как уже отмечалось, ваша проблема, скорее всего, связана с начальными значениями. Вы можете использовать две стратегии:

  1. Используйте грубую силу, чтобы найти начальные значения. См. Пакет nls2, чтобы узнать, как это сделать.
  2. Постарайтесь получить разумное предположение о начальных значениях. В зависимости от ваших значений можно линеаризовать модель.

G = (1 / (1 + (B^2)*(R^2))^((C+1)/2))

ln(G)=-(C+1)/2*ln(B^2*R^2+1)

Если B ^ 2 * R ^ 2 большой, это становится ок. ln (G) = - (C + 1) * (ln (B) + ln (R)), что является линейным.

Если B ^ 2 * R ^ 2 близко к 1, это прибл. ln (G) = - (C + 1) / 2 * ln (2), что является постоянным.

(Пожалуйста, проверьте наличие ошибок, это было поздно ночью из-за футбольного матча.)

Изменить после предоставления дополнительной информации: данные выглядят так, как будто они соответствуют кумулятивной функции распределения. Если крякает, как утка, скорее всего, это утка. И фактически ?Gest заявляет, что CDF оценивается.

library(spatstat)
data(simdat)
simdat.Gest <- Gest(simdat)
Gvalues <- simdat.Gest$rs
Rvalues <- simdat.Gest$r
plot(Gvalues~Rvalues)

#let's try the normal CDF
fit <- nls(Gvalues~pnorm(Rvalues,mean,sd),start=list(mean=0.4,sd=0.2))
summary(fit)
lines(Rvalues,predict(fit))
#Looks not bad. There might be a better model, but not the one provided in the question.
person Roland    schedule 28.06.2012
comment
Это ОЧЕНЬ полезно, спасибо! Я могу изучить другие модели и попытаться оптимизировать ту, с которой работал ранее, но это отличная отправная точка. - person MikeZ; 28.06.2012
comment
Роланд, Моделирование моей информации с помощью CDF оказалось намного проще, но не предоставляет столько полезной информации, как я надеялся. Функция в моем первоначальном вопросе, смог бы я смоделировать ее с данными, предоставила бы более полезную информацию. Знаете ли вы, какие стратегии я мог бы использовать, чтобы найти источник ошибок, о которых я говорил выше? - person MikeZ; 05.07.2012
comment
Попытка подобрать модель, которая просто не очень хорошо описывает форму данных, обычно создает проблемы оптимизации. По сути, это то, о чем вам сообщают сообщения об ошибках. Выбор модели нелинейной регрессии должен (я бы почти сказал, что должен) основываться на физических или математических соображениях. Чем вы оправдываете свой выбор модели? Я подозреваю, что вы пытаетесь решить какую-то неизвестную проблему, которую можно было бы лучше решить с помощью других инструментов. - person Roland; 05.07.2012