tpcurves2: T-P Curves

# Problems with JAGS!

# Hi Matthieu! Thanks in advance for taking a look!

# So among the non-normal distributions I was looking at, I stumbled upon the Pearson IV distribution.
# Here is a function that replicates the Probability Density Function (PDF) in R (Extracted from Nagahara 2004)
# You'll need the "gsl" package for it to run btw!

library(gsl)

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)

}

# As you see, it has a complex gamma function in it (The "lngamma_complex" function is very handy for that), and
# as you can imagine, it was much easier for me to code it once in a separate function por repeatability.

# Now let's say that I have this data

T <- c(20:30) # temperature
P <- c(0, 5.4, 8.5, 10.2, 9.9, 11.4, 8.15, 7.65, 5, 0.5, 0) # performance

# And that I want to run a non-linear model to estimate the parameters of the pIV function based on this data

# I build the model in JAGS as usual.

library(R2jags) # I'm loading the model here just in case.

# Prepare the data
data <- list(t = T, p = P)

# Build the model
model <- function(){

  # Set the priors ( Pretty flat, I know, I don't wanna make any assumptions for now)

  mu ~ dunif(0, 100) # mu parameter, realistically it can only be positive

  tau ~ dunif(0, 50) # same thing here, in Pearson's IV this is a proxy for variance

  delta ~ dunif(-50,50) # this is a proxy for skewness so it can be negative

  b ~ dunif(4,100) # this is a proxy for kurtosis and it can only be >4

  sigma ~ dunif(0,1000)

  tau_model <- 1/sigma^2

  # Likelihood

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

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

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

  }

}

# Set the parameters to save
params <- c("mu", "tau", "delta", "b")

# And now when I fit the model...
fit <- jags(data = data, model.file = model,parameters.to.save = params,
            n.chains = 4, n.iter = 100000, n.burnin=20000)

# I calls pIV as an unknown function! What the hell? It looks like if you don't provide the full code of the function
# to the jags.model object it won't run it.

# How can an external custom function work within the jags model??

# Thanks for taking a look!