Problem with coding a new probability distribution in Greta

Hi everyone,

I am trying to use the generalized normal distribution (https://en.wikipedia.org/wiki/Generalized_normal_distribution) in a Greta model. As this distribution is not available, I have tried to follow the only example I have met so far, which is here:https://mdscheuerell.github.io/gretaDFA/[https://mdscheuerell.github.io/gretaDFA/].

The probelm I have is that when I run the Greta model I have the following error:
Error in self$tf_distrib(parameters, dag)$log_prob(x) :
object ‘tf’ cannot be found

Any idea on how to change my code? Here is the code I wrote:

library (R6)
distrib <- .internals$nodes$constructors$distrib
distribution_node <- .internals$nodes$node_classes$distribution_node
check_dims <- .internals$utils$checks$check_dims
as.greta_array <- .internals$greta_arrays$as.greta_array
is_scalar <- .internals$utils$misc$is_scalar
fl <- .internals$utils$misc$fl

function to call

gnorm <- function (mu, sigma, beta, dim = 1) {
distrib(“gnorm”, mu, sigma, beta, dim)
}

gnorm_distribution <- R6Class (
“gnorm_distribution”,
inherit = distribution_node,
public = list(

initialize = function (mu, sigma, beta, dim) {
  
  # (dim)
  if (length(dim) > 1 ||
      dim <= 0 ||
      !is.finite(dim) ||
      dim != floor(dim)) {
    
    stop ("dim must be a scalar positive integer, but was: ",
          capture.output(dput(dim)),
          call. = FALSE)
    
  }
  
 

  # check if sigma is scalar
  # (sigma)
  if (length(sigma) > 1) {
    
    stop ("sigma must be a scalar positive integer, but was: ",
          capture.output(dput(sigma)),
          call. = FALSE)
    
  }
  
  

  # check if mu is scalar
  # (mu)
  if (length(mu) > 1) {
    
    stop ("mu must be a scalar positive integer, but was: ",
          capture.output(dput(mu)),
          call. = FALSE)
    
  }
  
  

  # check if beta is scalar
  # (beta)
  if (length(beta) > 1) {
    
    stop ("beta must be a scalar positive integer, but was: ",
          capture.output(dput(beta)),
          call. = FALSE)
    
  }
  sigma <- as.greta_array(sigma)
  mu <- as.greta_array(mu)
  beta <- as.greta_array(beta)
  
  self$bounds <- c(0, Inf)
  super$initialize("gnorm", dim, truncation = c(0, Inf))
 self$add_parameter(sigma, "sigma")
 self$add_parameter(mu, "mu")
 self$add_parameter(beta, "beta")

},

tf_distrib = function (parameters, dag) {

  mu <- parameters$mu
  beta <- parameters$beta
  sigma <- parameters$sigma

  # log pdf(x | i, M, sigma)
  log_prob = function (x) {
    tf$log(beta) - tf$log(fl(2)*sigma)-tf$lgamma(1/beta) - tf$pow(tf$abs((x-mu)/sigma), (beta))
  }

  list(log_prob = log_prob, cdf = NULL, log_cdf = NULL)

},

# no CDF for discrete distributions
tf_cdf_function = NULL,
tf_log_cdf_function = NULL

)
)

Sincerely,

Frédéric Gosselin

Hi Frederic,
You can check @nick description on how to add new distributions to greta in this old pull request. The easiest would be to fork the repository, add your function code there and then load this package in R via devtools::install_github.

I tried this out with a zero-inflated poisson: https://github.com/lionel68/greta/commit/6aa4080e2979baf3e97952d0fd23241d73c91d4f, you can check there where the function code should go. I also recommend to add so very simple tests to check that the distribution function you created leads to similar density estimations than established function in other R packages. You can check the different tests in the test/ subfolder on the repo.

Once you’ve done this best would be to open a pull request to the main repo so that your implementation can be checked and made available to the rest of the world.

Hope this helps.
Lionel

Thank you very much Lionel. Very helpful.

I will try to follow your example (yet, I am not a Github expert so may take a little time).

All the best,

Frédéric

Dear all,
I have followed the initial line of my program. The problem is that some times (in sme R sessions), the program works and not in others.

Here is below a program easier to replicate.

Sincerely,

Frédéric

#Packages
library(coda)

Sys.setenv(RETICULATE_PYTHON = “C:/Users/frederic.gosselin/Miniconda3/envs/r-tensorflow”)
reticulate::use_condaenv(“r-tensorflow”)
library(greta)

library (R6)
distrib <- greta::.internals$nodes$constructors$distrib
distribution_node <- greta::.internals$nodes$node_classes$distribution_node
check_dims <- greta::.internals$utils$checks$check_dims
as.greta_array <- greta::.internals$greta_arrays$as.greta_array
is_scalar <- .internals$utils$misc$is_scalar
fl <- .internals$utils$misc$fl

function to call

gnorm <- function (mu, sigma, beta0, dim = 1) {
distrib(“gnorm”, mu, sigma, beta0, dim)
}

Leung & Drton distn for prior of diag(Z)

gnorm_distribution <- R6Class (
“gnorm_distribution”,
inherit = distribution_node,
public = list(

initialize = function (mu, sigma, beta0, dim) {
  

  sigma <- as.greta_array(sigma)
  mu <- as.greta_array(mu)
  beta0 <- as.greta_array(beta0)
 dim <- check_dims(mu, sigma, beta0, target_dim = dim)

 super$initialize("gnorm", dim)
 self$add_parameter(sigma, "sigma")
 self$add_parameter(mu, "mu")
 self$add_parameter(beta0, "beta0")

},

tf_distrib = function (parameters, dag) {

  mu <- parameters$mu
  beta0 <- parameters$beta0
  sigma <- parameters$sigma

  # log pdf(x | i, M, sigma)
  log_prob = function (x) {
  
   tf$log(beta0)- tf$log(fl(2)*sigma)-tf$lgamma(fl(1)/beta0) - tf$pow(tf$abs((x-mu)/sigma), (beta0))
         }

  list(log_prob = log_prob, cdf = NULL, log_cdf = NULL)

},

# no CDF for discrete distributions
tf_cdf_function = NULL,
tf_log_cdf_function = NULL

)
)

set.seed(1)
Y<-rnorm(100)
sigma<-normal(0,2,truncation=c(0,Inf))
beta0<-normal(2,2,truncation=c(0,Inf))
mu<-normal(0,1)
Y<-as_data(Y)
distribution(Y)<-gnorm(mu,sigma,beta0)
m<-model(mu,sigma,beta0)

####MO154writing4
set.seed(1)
time1bc<-system.time(test1bc<-greta::mcmc(m,sampler = hmc(Lmin = 10, Lmax = 50)))