2015-09-29 5 views
0

Я пытаюсь решить двухкомпонентную модель распада в R, используя функцию nls, но столкнувшись с ошибками. Уравнение:проблемы с nls в R

two-component model

где т время, Ctot является С1 + С2, а p1 и p2 известны пропорции Ctot.

мои данные (дд) является:

> head(dd,n=15) 
     t Ctot 
1 0.00 6.62 
2 0.33 6.45 
3 0.50 6.38 
4 0.67 6.44 
5 0.83 6.38 
6 1.00 6.39 
7 1.17 6.35 
8 1.33 6.33 
9 1.50 6.33 
10 1.67 6.28 
11 1.83 6.17 
12 2.00 6.11 
13 2.17 6.07 
14 2.33 5.89 
15 2.50 5.86 

Использование NLS Я пробовал:

p1 <- 0.3 
p2 <- 0.7 
z <- nls(Ctot~(p1*C1*(exp(-k1*t)))+(p2*C2*(exp(-k2*t))), data=dd, start=list(C1=6, C2=0.1, k1=0.01, k2=0.01)) 

Однако я получаю:

z <- nls(Ctot~(p1*C1*(exp(-k1*t)))+(p2*C2*(exp(-k2*t))), data=dd, start=list(C1=6, C2=0.1, k1=0.01, k2=0.01)) 
Error in numericDeriv(form[[3L]], names(ind), env) : 
    Missing value or an infinity produced when evaluating the model 

Я был бы признателен, если кто-то имеет предложения!

+0

Может быть, что оптимальные значения '' k1' и k2' близки друг другой, чтобы модель с одним термином была бы лучше. –

+0

@ G.Grothendieck Да, но это двухкомпонентная модель распада. – squishy

+0

Если k1 = k2, то это не так. –

ответ

1

Данные кажутся довольно ограниченными и явно неполными, поскольку это только голова. Если мы делаем какие-то данные методы тестирования ... и выйти из запутанных p1 и p2:

t=seq(0, 20, by=.3) 
Ctot = 3 * exp(-1 * t) + 4 * exp(-5*t) 
# following hte example on gnm::gnm's help page: 

saved.fits <- list(); library(gnm) 
for (i in 1:10) { 
     saved.fits[[i]] <- suppressWarnings(gnm(Ctot ~ Exp(1 + t, inst = 1) + 
                 Exp(1 + t, inst = 2), 
        verbose=FALSE))} 
plot(Ctot~t) 
lines(saved.fits[[3]]$fitted~t) 
lines(saved.fits[[3]]$fitted~t,col="red") 

Я не был знаком с пакетом ГНСА и поэтому в конечном итоге прочитав первые несколько разделов, а затем обработанным 2 пример компоновки данных компонента в виньетке: https://cran.r-project.org/web/packages/gnm/vignettes/gnmOverview.pdf. Большинство припадков будет, как и ожидалось, но некоторые найдут локальный максимум в вероятности того, что не является глобальной макс:

> saved.fits[[1]]$coefficients 
        (Intercept) Exp(. + t, inst = 1).(Intercept) 
        1.479909e-12      1.098612e+00 
      Exp(1 + ., inst = 1).t Exp(. + t, inst = 2).(Intercept) 
        -1.000000e+00      1.386294e+00 
      Exp(1 + ., inst = 2).t 
        -5.000000e+00 
attr(,"eliminated") 
[1] 0 
> exp(saved.fits[[1]]$coefficients[4]) 
Exp(. + t, inst = 2).(Intercept) 
           4 
> exp(saved.fits[[1]]$coefficients[2]) 
Exp(. + t, inst = 1).(Intercept) 
           3 

 Смежные вопросы

  • Нет связанных вопросов^_^