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))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.