sandbox/trials.R

In ggcostoya/momentum:

# Vector of values provided
tseq <-seq(0,100,by=0.1)

# Predictor functions
gauss <- pred_gaussian(x, s = 10, a = 25, b = 1) # Gaussian function prediction
esgn <- pred_esgn(x, s = 45, a = 25, b = 1.9, c = 100, d = -3, e = 3)

emg <- pred_emg(tseq, s = 70, a = 25, b = 0.1, c = 1)
plot(tseq,emg,type="l")

esgn <- pred_esgn(tseq, s = 40, a = 25, b = 5, c = 1, d = 50, e =100)
plot(tseq,esgn,type="l")

w <- pred_weibull(tseq, 16.08, 38.20, 450000,71200)
plot(tseq,w,type="l")


# Plotting
plot(x,mg, type = "l", col = "grey", xlab = "Vector 'x' of provided values", ylab = "Prediction by the Gaussian function")
lines(x,emg, type = "l", col = "grey")
lines(x,esgn, type = "l", col = "grey")
lines(x,gauss, type = "l", lwd = 2)

# Data generating function

gen_tpd <- function(Topt, CTmax)

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)

}

x <- gen_base(25, 2, 1.25, 0, 0, 10, 5, 27, 23)

plot(x, type ="o")



fit_gaussian <- function(data){

  # Function to determine likelihood
  lk_gaussian <- function(s,a,b,sigma){

    x <- data$t # assign temperature variable
    y <- data$p # assign performance variable

    yp <- pred_gaussian(x,s,a,b) # predicted y

    loglik <- dnorm(y, mean = yp, sd = sigma, log = TRUE) # likelihood

    return(-sum(loglik))
  }

  # Find good estimates from the data
  start_s <- max(data$p, na.rm = T) # starting s
  start_a <- mean(data$t, na.rm = T) # starting a
  start_b <- sd(data$t, na.rm = T) # starting b
  low_a <- min(data$t, na.rm = T) # low a
  up_a <- max(data$t, na.rm = T) # Upper a

  # Actual fit
  fit <- mle(lk_gaussian, method = "L-BFGS-B",
             start = list(s = start_s, a = start_a, b = start_b, sigma = 10),
             lower = c(1, low_a, 0.1, 0.1),
             upper = c(100, up_a, 100, 100))

  # Extract fit coefficients
  coefs <- as.data.frame(t(coef(fit)))

  # Get the AIC score
  aic <- AIC(fit)

  # Join coefficients and AIC score
  table <- cbind(coefs, aic)

  return(table)
}


fitaic_gaussian(x)
fitaic_emg(x)
fitaic_weibull(x)

gline <- pred_gaussian(tseq, 9.5139, 25, 1.7102)
emgline <- pred_emg(tseq,20.4677,24.3733,1.6,0.69)
wline <- pred_weibull(tseq,7.57,25,4e6,5e5)

plot(x, pch = 19)
lines(tseq,gline, type = "l", lwd = 2)
lines(tseq,emgline, type ="l", lwd =2, col ="red")
lines(tseq,wline,type="l",lwd =2, col = "green")