sandbox/pearsonIV.R

In ggcostoya/tpcurves2: T-P Curves

install.packages("gsl")

library(gsl)
library(PearsonDS)

lngamma_complex()
lngamma()



pIV <- function(x,mu,tau,delta,b){

  # first part top
  f_top <- exp(lngamma_complex(b + b*delta*1i))*exp(lngamma_complex(b - b*delta*1i))*tau^(2*b - 1)

  # first part bottom
  f_bot <- exp(lngamma(b))*exp(lngamma(b - 0.5))*pi^(1/2)

  # second part top
  s_top <- exp(2*b*delta*atan((x-mu)/tau))

  # second part bottom
  s_bot <- ((x - mu)^2 + tau^2)^b

  # putting it together to generate the prediction
  y <- (f_top/f_bot)*(s_top/s_bot)

  # avoid returning a complex vector
  y <- as.numeric(y)

  #return
  return(y)

}

get_mean_pIV <- function(mu,tau,delta,b){

  mean <- mu + (b*delta*tau)/(b-1)

  return(mean)
}

get_var_pIV <- function(tau,delta,b){

  var <- (tau^2/(2*b-3))*(1 + ((b*delta)/(b-1))^2)

  return(var)
}

get_skw_pIV <- function(delta,b){

  skw <- (2*b*(sqrt(2*b-3))/((b-1)*(b-2)))*(delta/(1 + ((b*delta)/(b-1))^2))

  return(skw)

}

get_kurt_pIV <- function(delta,b){

  kurt_f <- (3*(2*b-3)/(2*b-5))*((1 + ((b+2)/(b-2))*((b*delta)/(b-1))^2)/(1 + ((b*delta)/(b-1))^2))

}

get_params_pIV <- function(M,SD,S,K){

  V <- SD^(2) # transform standard deviation into variance

  b <- (9+6*S^2-5*K)/(6+3*S^2-2*K)   # get the "b" shape parameter

  tau <- (0.5)*(sqrt(V*(4*(2*b-3)-S^2*(b-2)^2)))  # get the "tau" scale parameter

  delta <- (sqrt(V)*S*(b-1)*(b-2))/(2*b*tau)   # get the "delta" shape parameter

  mu <- M - (b*delta*tau)/(b -1)   # get the "mu" location parameter

  params <- list(b = b, tau = tau, delta = delta, mu = mu)   # generate the params list

  return(params)

}

get_moments_pIV <- function(mu,tau,delta,b){

  M <- get_mean_pIV(mu = mu, tau = tau, delta = delta, b = b) # get mean

  V <- get_var_pIV(tau = tau, delta = delta, b = b) # get variance

  SD <- sqrt(V) # transform variance into standard deviation

  SKW <- get_skw_pIV(delta = delta, b = b) # get skewness

  KURT <- get_kurt_pIV(delta = delta, b = b) # get kurtosis

  moments <- list(mean = M, var = V, sd = SD, skw = SKW, kurt = KURT) # generate a list with the moments

  return(moments)

}

get_params_pIV(M = 25, SD = 5, S = 0, K = 5)

tmp <- seq(10,40,by=0.01)

mu = 25

tau = 5

delta = 0.1

b = 4

pred <- pIV(x = tmp, mu = mu, tau = tau, delta = delta, b = b)

plot(tmp, pred, type = "l")

abline(v = mu, lty = 2)

mean <- mu + (b*delta*tau)/(b - 1)

abline(v = mean, lty = 2)

max(pred)


# let's try to fit a Pearson IV model

# load jags
library(R2jags)

# generate the data
data <- gen_tpd(Topt = 25, CTmax = 30, CTmin = 20, Pmax = 10, Pmin = 0, Error = 1, Samples = 20, Degree = 3)

plot(data)

# prepare the data
jags.data <- list(t = data$Tm, p = data$P)

# build the model
jags.model <- function(){

  # Moments priors
  M ~ dunif(20,30) # mean

  SD ~ dunif(0,50) # standard deviation

  S ~ dnorm(0,5) # mean skewness is zero an it is not really skewed

  K ~ dunif(4, 100) # kurtosis needs to be positive

  # prior transformation into params

  V <- SD^(2) # transform standard deviation into variance

  b <- (9+6*S^2-5*K)/(6+3*S^2-2*K)   # get the "b" shape parameter

  tau <- (0.5)*(sqrt(V*(4*(2*b-3)-S^2*(b-2)^2)))  # get the "tau" scale parameter

  delta <- (sqrt(V)*S*(b-1)*(b-2))/(2*b*tau)   # get the "delta" shape parameter

  mu <- M - (b*delta*tau)/(b -1)   # get the "mu" location parameter

  # prior for observations

  sigma ~ dunif(0,1000) # priors for sigma

  taupred <- 1/sigma^2

  # call the function pIV within the model?
  pIV <- function(x,mu,tau,delta,b){

    # first part top
    f_top <- exp(lngamma_complex(b + b*delta*1i))*exp(lngamma_complex(b - b*delta*1i))*tau^(2*b - 1)

    # first part bottom
    f_bot <- exp(lngamma(b))*exp(lngamma(b - 0.5))*pi^(1/2)

    # second part top
    s_top <- exp(2*b*delta*atan((x-mu)/tau))

    # second part bottom
    s_bot <- ((x - mu)^2 + tau^2)^b

    # putting it together to generate the prediction
    y <- (f_top/f_bot)*(s_top/s_bot)

    # avoid returning a complex vector
    y <- as.numeric(y)

    #return
    return(y)

  }

  # Likelihood
  for(i in 1:length(t)){

    # prediction for p
    p_pred[i] <- pIV(t[i], mu, tau, delta, b)

    # observed p
    p[i] ~ dnorm(p_pred[i], taupred)
  }

}


fit <- jags.parallel(data = jags.data, model.file = jags.model,
                     parameters.to.save = c("M","SD","S","K"),
                     n.iter = 10000, n.burnin = 3300, n.chains = 3)