Ниже мы запускаем nls (с использованием слегка измененной модели) и nlxb (из nlsr), но nlxb останавливается до схождения. Несмотря на эти проблемы, обе они, тем не менее, дают результаты, которые хорошо визуально соответствуют данным. Эти проблемы предполагают наличие проблем с самой моделью, поэтому в разделе Другое, руководствуясь выходными данными nlxb, мы показываем, как исправить модель, давая подмодель исходной модели, которая легко соответствует данным с помощью как nls, так и nlxb, а также хорошо подходят. В конце раздела Примечания мы предоставляем данные в воспроизводимой форме.
nls
Предполагая, что установка воспроизводимо показана в примечании в конце, переформулируйте проблему для алгоритма nls plinear, определив правую матрицу, столбцы которой умножают каждый из линейных параметров, a1, a2, a3 и d, соответственно. plinear не требует начальных значений для упрощающих настройку. Он сообщит о них как .lin1, .lin2, .lin3 и .lin4 соответственно.
Чтобы получить начальные значения, мы использовали более простую модель без группировки и поиск по сетке по b от 1 до 10 и c также от 1 до 10, используя nls2 в одноименном пакете. Мы также обнаружили, что nls по-прежнему выдает ошибки, но при использовании abs в формуле, как показано, все работает до завершения.
Проблемы с моделью предполагают, что с ней существует фундаментальная проблема, и в разделе «Другое» мы обсуждаем, как ее исправить.
xx <- c(x, x, x)
yy <- c(y1, y2, y3)
# startingi values using nls2
library(nls2)
fo0 <- yy ~ cbind(b ^ abs(xx - c), 1)
st0 <- data.frame(b = c(1, 10), c = c(1, 10))
fm0 <- nls2(fo0, start = st0, alg = "plinear-brute")
# run nls using starting values from above
g <- rep(1:3, each = length(x))
fo <- yy ~ cbind((g==1) * b ^ abs(xx - c[g]),
(g==2) * b ^ abs(xx - c[g]),
(g==3) * b ^ abs(xx - c[g]),
1)
st <- with(as.list(coef(fm0)), list(b = b, c = c(c, c, c)))
fm <- nls(fo, start = st, alg = "plinear")
plot(yy ~ xx, col = g)
for(i in unique(g)) lines(predict(fm) ~ xx, col = i, subset = g == i)
fm
давая:
Nonlinear regression model
model: yy ~ cbind((g == 1) * b^abs(xx - c[g]), (g == 2) * b^abs(xx - c[g]), (g == 3) * b^abs(xx - c[g]), 1)
data: parent.frame()
b c1 c2 c3 .lin1 .lin2 .lin3 .lin4
1.997 0.424 1.622 1.074 0.680 0.196 -0.532 9.922
residual sum-of-squares: 133
Number of iterations to convergence: 5
Achieved convergence tolerance: 5.47e-06
(продолжение после сюжета)
nlsr
С nlsr это будет сделано так. Поиск по сетке для начальных значений не требовался, и добавление abs также не требовалось. Значения b и d кажутся похожими на решение nls, но другие коэффициенты отличаются. Визуально оба решения кажутся подходящими для данных.
С другой стороны, из столбца JSingval мы видим, что якобиан имеет недостаточный ранг, из-за чего он останавливался и не производил значений SE, и сходимость сомнительна (хотя этого может быть достаточно, учитывая, что визуально график, не показанный, выглядит как хорошо подходит). Мы обсуждаем, как это исправить, в разделе «Другое».
g1 <- g == 1; g2 <- g == 2; g3 <- g == 3
fo2 <- yy ~ g1 * (a1 * b ^ (xx - c1) + d) +
g2 * (a2 * b ^ (xx - c2) + d) +
g3 * (a3 * b ^ (xx - c3) + d)
st2 <- list(a1 = 1, a2 = 1, a3 = 1, b = 1, c1 = 1, c2 = 1, c3 = 1, d = 1)
fm2 <- nlxb(fo2, start = st2)
fm2
давая:
vn: [1] "yy" "g1" "a1" "b" "xx" "c1" "d" "g2" "a2" "c2" "g3" "a3" "c3"
no weights
nlsr object: x
residual sumsquares = 133.45 on 153 observations
after 16 Jacobian and 22 function evaluations
name coeff SE tstat pval gradient JSingval
a1 3.19575 NA NA NA 9.68e-10 4097
a2 0.64157 NA NA NA 8.914e-11 662.5
a3 -1.03096 NA NA NA -1.002e-09 234.9
b 1.99713 NA NA NA -2.28e-08 72.57
c1 2.66146 NA NA NA -2.14e-09 10.25
c2 3.33564 NA NA NA -3.955e-11 1.585e-13
c3 2.0297 NA NA NA -7.144e-10 1.292e-13
d 9.92363 NA NA NA -2.603e-12 3.271e-14
Мы можем рассчитать SE, используя nls2 в качестве второго этапа, но это все еще не решает проблему со всем, что предлагают сингулярные значения.
summary(nls2(fo2, start = coef(fm2), algorithm = "brute-force"))
давая:
Formula: yy ~ g1 * (a1 * b^(xx - c1) + d) + g2 * (a2 * b^(xx - c2) + d) +
g3 * (a3 * b^(xx - c3) + d)
Parameters:
Estimate Std. Error t value Pr(>|t|)
a1 3.20e+00 5.38e+05 0.0 1
a2 6.42e-01 3.55e+05 0.0 1
a3 -1.03e+00 3.16e+05 0.0 1
b 2.00e+00 2.49e-03 803.4 <2e-16 ***
c1 2.66e+00 9.42e-02 28.2 <2e-16 ***
c2 3.34e+00 2.43e+05 0.0 1
c3 2.03e+00 8.00e+05 0.0 1
d 9.92e+00 4.42e+05 0.0 1
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.959 on 145 degrees of freedom
Number of iterations to convergence: 8
Achieved convergence tolerance: NA
Другой
Когда у nls возникают проблемы с подгонкой модели, это часто предполагает, что что-то не так с самой моделью. Немного поигравшись с этим, руководствуясь столбцом JSingval в выходных данных nlsr выше, который предполагает, что c параметры или d могут быть проблемой, мы обнаруживаем, что если мы исправим все значения параметров c на 0, то модель будет легко соответствовать при достаточно хорошем начальные значения, и он по-прежнему дает низкую остаточную сумму квадратов.
library(nls2)
fo3 <- yy ~ cbind((g==1) * b ^ xx, (g==2) * b ^ xx, (g==3) * b ^ xx, 1)
st3 <- coef(fm0)["b"]
fm3 <- nls(fo3, start = st3, alg = "plinear")
давая:
Nonlinear regression model
model: yy ~ cbind((g == 1) * b^xx, (g == 2) * b^xx, (g == 3) * b^xx, 1)
data: parent.frame()
b .lin1 .lin2 .lin3 .lin4
1.9971 0.5071 0.0639 -0.2532 9.9236
residual sum-of-squares: 133
Number of iterations to convergence: 4
Achieved convergence tolerance: 1.67e-09
который, как указывает следующая аниова, сравним с fm сверху, несмотря на то, что у него на 3 параметра меньше:
anova(fm3, fm)
давая:
Analysis of Variance Table
Model 1: yy ~ cbind((g == 1) * b^xx, (g == 2) * b^xx, (g == 3) * b^xx, 1)
Model 2: yy ~ cbind((g == 1) * b^abs(xx - c[g]), (g == 2) * b^abs(xx - c[g]), (g == 3) * b^abs(xx - c[g]), 1)
Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
1 148 134
2 145 133 3 0.385 0.14 0.94
Мы можем повторить fm3, используя nlxb, вот так:
fo4 <- yy ~ g1 * (a1 * b ^ xx + d) +
g2 * (a2 * b ^ xx + d) +
g3 * (a3 * b ^ xx + d)
st4 <- list(a1 = 1, a2 = 1, a3 = 1, b = 1, d = 1)
fm4 <- nlxb(fo4, start = st4)
fm4
давая:
nlsr object: x
residual sumsquares = 133.45 on 153 observations
after 24 Jacobian and 33 function evaluations
name coeff SE tstat pval gradient JSingval
a1 0.507053 0.005515 91.94 1.83e-132 8.274e-08 5880
a2 0.0638554 0.0008735 73.11 4.774e-118 1.26e-08 2053
a3 -0.253225 0.002737 -92.54 7.154e-133 -4.181e-08 2053
b 1.99713 0.002294 870.6 2.073e-276 -2.55e-07 147.5
d 9.92363 0.09256 107.2 3.367e-142 -1.219e-11 10.26
Примечание
Предполагаемый ввод ниже такой же, как и в вопросе, за исключением того, что мы дополнительно установили начальное число, чтобы сделать его воспроизводимым.
set.seed(123)
my_model <- function(x, a, b, c, d) a * b ^ (x - c) + d
x <- seq(0, 10, 0.2)
b <- 2; d <- 10 # shared
a1 <- 1; c1 <- 1
y1 <- my_model(x, a = a1, b = b, c = c1, d = d) + rnorm(length(x))
a2 <- 2; c2 <- 5
y2 <- my_model(x, a = a2, b = b, c = c2, d = d) + rnorm(length(x))
a3 <- -2; c3 <- 3
y3 <- my_model(x, a = a3, b = b, c = c3, d = d) + rnorm(length(x))
person
G. Grothendieck
schedule
22.08.2020
