Самозапускающиеся функции для модели экспоненциального убывания в R - PullRequest
Новая
       75

Самозапускающиеся функции для модели экспоненциального убывания в R

1 голос
/ 21 апреля 2020

Я работаю над моделью экспоненциального затухания, где я хотел бы оценить скорость затухания. Моя текущая модель использует функцию автозапуска, SSasymp из пакета stats. Я также написал вторую модель, где я просто проверяю начальные параметры, для которых требуется пакет minpack.lm. У меня вопрос, есть ли другой способ оценить начальные параметры для перекрестной проверки функции SSasymp. Я (думаю), я понимаю, что делает код для оценки начальных параметров, но я хотел получить некоторую обратную связь о том, является ли SSasymp правильной функцией для использования с этими данными или если есть другая функция, которую я мог бы потенциально использовать. 

library(stats)
library(minpack.lm)
library(broom)
library(ggplot2)

df<-data.frame(Date=seq(1:66),
           Level=c(1438072839.75,   1397678053.5,   1358947420.5,   1313619938.25,  1269224528.25, 
1246776954.75,  1207201162.5,   1176229091.25,  1136063160, 1103721704.25,  1080591637.5,    
1048286667, 1017840460.5,   1001057052, 975815001,  943568665.5,    932026210.5,    916996593.75,    
903904288.5,    887578544.25,   871428547.5,    855417720,  843504839.25,   825835607.25,    
816060303.75,   803506361.25,   801213123,  797977217.25,   793483994.25,   780060123,  766265609.25,    
756172471.5,    746615497.5,    738002936.25,   723741644.25,   711969181.5,    696032998.5,     
686162453.25,   671953166.25,   674184571.5,    664739475,  641091932.25,   627358484.25,    
616740068.25,   602261552.25,   592440797.25,   584160403.5,    569780103.75,   556305753,   
551682927,  546535062,  537782506.5,    524251944.75,   519277188.75,   503598795,  498481312.5,     
487907885.25,   479760227.25,   474773064.75,   468246932.25,   460561701,  455266345.5,     
448451890.5,    447760119,  441236056.5,    438884417.25))

dfDecay<-nls(Level~ SSasymp(Date, Asym, R0, lrc), data = df)
dfFitted<-augment(dfDecay)
ggplot(df, aes(x=Date,y=Level))+geom_point()+  geom_line( aes(y=dfFitted$.fitted), color="red")

dfDecay2<-nlsLM(Level~b*exp(-a*Date), 
                   data = df,
                   start= list(a=.01,b=1.5e+09),
                   algorithm = "LM")
fitDecay2<-augment(dfDecay2)
ggplot(df, aes(x=Date,y=Level))+geom_point()+  geom_line( aes(y=fitDecay2$.fitted), color="red")

1 Ответ

1 голос
/ 23 апреля 2020

Относительно начальных значений:

  1. Возьмите бревна с обеих сторон и поместите в линейную модель.
  2. Параметры должны быть одинаковой величины, чтобы избежать численных проблем, поэтому используйте Level/1e9 вместо Level. Это просто изменяет единицы измерения уровня.
  3. Используя начальные значения из линейной модели, достаточно nls.

Это дает: 

fm0 <- lm(log(Level/1e9) ~ Date, df)
st <- list(a = exp(coef(fm0)[[1]]), b = -coef(fm0)[[2]])
nls(Level/1e9 ~ a * exp(-b * Date ), df, start = st)

дача: 

Nonlinear regression model
  model: Level/1e+09 ~ a * exp(-b * Date)
   data: df
     a      b 
1.3532 0.0183 
 residual sum-of-squares: 0.08055

Number of iterations to convergence: 4 
Achieved convergence tolerance: 4.023e-07
...