sandbox/trials_nls.R

tseq <- seq(0,100,by=0.1)

gen_base <- function(Topt, Tb50, Tb80, Skw50, Skw80, Pmax, Pmin, CTmax, CTmin){

  # temperature
  t <- c(CTmin,
         (Topt - Tb50/2 + Skw50),
         (Topt - Tb80/2 + Skw80),
         Topt,
         (Topt + Tb80/2 + Skw80),
         (Topt + Tb50/2 + Skw50),
         CTmax)

  # performance
  p <- c(Pmin,
         mean(c(Pmax,Pmin)),
         ((Pmax-Pmin)*0.8 + Pmin),
         Pmax,
         ((Pmax-Pmin)*0.8 + Pmin),
         mean(c(Pmax,Pmin)),
         Pmin)

  # bind datasets
  tpd <- data.frame(t,p)

  return(tpd)

}

pred_weibull <-function(x,s,a,b,c){

  cs <- (c-1)/c
  xab <- (x-a)/b
  y <- s*cs^((1-c)/c)*((xab + (cs^(1/c)))^(c-1))*exp(-(xab + (cs^(1/c))^(c-1))^c)+cs

  return(y)

}

tpd <- gen_base(25, 2, 1, 0.25, 0.15, 10, 5, 27, 23)
tpd <- rbind(tpd,tpd+0.1,tpd+0.05,tpd-0.05,tpd-0.1)

plot(tpd,pch = 19)

start_s <- max(tpd$p, na.rm = T)
start_a <- mean(tpd$t[tpd$p == max(d$p, na.rm = TRUE)])
start_b <- max(tpd$t, na.rm = TRUE) - min(tpd$t, na.rm = TRUE)
start_c <- 4

start_vals <- c(s = start_s, a = start_a, b = start_b, c = start_c)


low_s <- 1
low_a <- min(tpd$t, na.rm = TRUE)
low_b <- 0.1
low_c <- 0.1

low_limit <- c(low_s,low_a,low_b,low_c)

up_s <- 50 # times 50 is more than enough
up_a <- max(tpd$t, na.rm = TRUE)
up_b <- Inf
up_c <- Inf

up_limit <-c(up_s, up_a, up_b, up_c)


A <- nls(p ~ pred_weibull(x = t, s, a, b, c),
    data = tpd,
    start = start_vals,
    lower = low_limit,
    upper = up_limit,
    algorithm = "port",
    trace = TRUE)

emg <- pred_emg(tseq,18.766,24.994,1.747,0.1)
weibull <- pred_weibull(tseq, 11.460, 24.435, 8.762, 6.109)
plot(tpd,pch=19)
lines(tseq, weibull, type = "l", lwd = 2)
lines(tseq,emg,type="l", lwd = 2, col = "red")

AIC(fit_emg(tpd))
AIC(fit_weibull(tpd))
ggcostoya/momentum documentation built on Feb. 14, 2021, 6:12 p.m.