1

I'm using the BayesianTools package to run some MCMCs on data, but my specific project requires me to take some outputs of one MCMC and convert it into a prior for a new MCMC. Priors need to be defined as single input functions as far as I'm aware in BayesianTools, so what I'm trying to do is write a function that takes the posteriors from one round of MCMCs and generates a probability density function as a result. I tried doing this using the code below, but it seems that R doesn't save dens inside of the function probabilityDensityFunction as I had hoped it would. I need to generate 16 of these objects and my life would be so much easier if I could automatically generate them in a manner like this instead of manually coding for each one, so any help is apreciated.

Thanks!

probabilityDensityFunction_generator <- function(dens) {
  # This function is to generate a separate function called the probability density function. 
  #It takes objects of type density() as an imput and returns a function which estimates the probability density. 
  
  probabilityDensityFunction <- function(x){
    if(x < dens$x[1]){
      out = 0
    } else if(x>dens$x[length(dens$x)]){
      out = 0
    } else{
      for(i in 1:length(dens$x)){
        if(x<dens$x){
          out = dens$y[i]
        } else{
          break
        }
      }
    }
    return(out)
  }
  
  return(probabilityDensityFunction)
}
andrew brian
  • 11
  • 2
  • 2
    It's very possible - there's a chapter on the topic in the free Advanced R book, [Function Factories](https://adv-r.hadley.nz/function-factories.html). I'd suggest reading that and doing the accompanying examples. – Gregor Thomas Jan 11 '23 at 14:39
  • Remove `probabilityDensityFunction <- ` and remove the last line (`return`, all of it). – Rui Barradas Jan 11 '23 at 14:41
  • So what exactly is wrong with this function? What errors do you get when you use it? R could capture variables in the closure when returning new function. It's easier to help you if you include a simple [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) with sample input and desired output that can be used to test and verify possible solutions. – MrFlick Jan 11 '23 at 14:54
  • The main problem I see in your code is that you never evaluate `dens` in `probabilityDensityFunction_generator`, so lazy evaluation means it might change by the time you call the function you're returning. Add a line `force(dens)` before creating the new function and that'll be fine. – user2554330 Jan 11 '23 at 16:45

1 Answers1

1

Here is a way.
The main error is assign the function to an object, there's no need for that, you should return the function creation.

There is another problem, the code to get where in dens$x is the new x, base function findInterval does it automatically. But it returns zero if it is below the minimum, so assign NA to those y's.

And due to R's lazy evaluation, dens must be evaluated in the function generator. This is done with force. If not when dens changes so will the all functions it created. See below that when applied to the same vector newx, f and g graphs are completely different.

probabilityDensityFunction_generator <- function(dens) {
  # This function is to generate a separate function called the probability density function. 
  # It takes objects of type density() as an imput and returns a function which estimates the probability density. 
  force(dens)
  function(x){
    i <- sapply(x, findInterval, dens$x)
    is.na(i) <- i == 0
    dens$y[i]
  }
}

set.seed(2023)
X <- rnorm(1e3)
d <- density(X)
f <- probabilityDensityFunction_generator(d)

Y <- rchisq(1e3, df = 3)
e <- density(Y)
g <- probabilityDensityFunction_generator(e)

xlim <- range(range(d$x), range(e$x))
ylim <- range(range(d$y), range(e$y))
plot(d, xlim = xlim, ylim = ylim)
lines(e)
#
newx <- seq(xlim[1], xlim[2], 0.2)
points(newx, f(newx), col = "red", pch = 3)
points(newx, g(newx), col = "blue", pch = 3)

Created on 2023-01-11 with reprex v2.0.2

Rui Barradas
  • 70,273
  • 8
  • 34
  • 66