Options for discrete/custom priors?

[Ryan Burner]

Hi All,

A quick question about priors. I've got a model where I specify informative priors for covariate values that are estimated quantities in the model. I do this because I know something about the values of these covariates but want to incorporate the uncertainty around those covariate estimates.

For each covariate value I have essentially a vector of possible values, where each possible values occurs one or more times (in proportion to its probability of being the true value). Currently I specify the prior for each covariate using the beta distribution that best fits the relevant vector. In some cases the beta distribution is very similar to that vector, but in other cases they are quite different.

Is there a way I can specify that vector of values as the prior directly, without converting it to a recognized, parameterized distribution? I could also of course convert the vector to e.g., a dataframe with a row for each value and a column that gave the probabilities of each.

Thanks!

Ryan    

[Daniel Turek]

Ryan, as the subject of your email suggests, it sounds a lot like you want a custom-written discrete prior distribution.  This is no problem at all in nimble.  Do you have any experience writing nimbleFunctions, or your own custom distributions using nimbleFunctions?

Doing what you describe would be straightforward, where the (custom written) discrete prior (density function) returns some particular value from the vector of probabilities, depending on the value of the random variate (x).  This would also be done on the log scale, so it would, usually, return the log of one of these probabilities.  Furthermore, you should write this density function to return 0 (-Inf on the log scale) if the x value is not among the set of possible values.

Another important point - you'd need to manually register your distribution with nimble (using the function registerDistributions), specifically in order to provide the argument "discrete = FALSE".  Importantly, this will induce the MCMC to sample these as discrete variates.

See Chapter 12 Creating user-defined BUGS distributions and functions of the nimble user manual.  If you need more help getting started, feel free to ask.

Cheers,

Daniel

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/b83043f4-d700-4659-bd0a-b65d95c49207n%40googlegroups.com.

[Ryan Burner]

Hi Daniel,

Thanks for the quick reply! Ok, yes, I can see how that would work after reading Chapter 11, although I might end up with questions about a bit of the syntax (and you meant 'discrete = TRUE', right?). But one other issue - this distribution varies for each of several thousand site_covariate combinations (i.e., it isn't consistent across sites for a given covariate). So, I'd actually need to define thousands of named custom distributions. I assume I could define them using a loop and some kind of consistent naming system, but then a) would that huge number of distributions bog down my model, and b) is there an efficient way to then assign those custom prior distributions to the relevant model parameters?

Currently to specify the priors using the appropriate beta distribution for each site_covariate I just have a site x covariate x distribution parameter array in model constants that contains the beta distribution parameters for each site and covariate, so it is easy to assign these priors in the model in a loop, but I can't envision how I'd do that when the names varied unless Nimble model syntax would allow e.g., some kind of paste(get( MyCustomDistName[i]  )) function that cycled through a list of distribution names?

Thanks!

Ryan  

[Daniel Turek]

Ryan, I'm sorry this is very brief.  But to answer your questions:

1. Yes, I certainly meant "discrete = TRUE", my apologies.

2. Make the custom distribution take an additional few arguments (for example, i = site number, and j = covariate number).  Then, in the internal logic of the distribution, it can use the i and j values to access the correct dimensions / values for that particular distribution.  Also, there are a few ways to provide the distributional information (array of values and probabilities) to the custom distribution.  For now, try providing them as an argument (an N-dimensional array) to the distribution.  There's also a new mechanism we can try out, to embed this information inside the custom distribution itself.  I'll help you flesh that out in due course, but I'm sorry I can't right now.

Cheers, keep in touch,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/7ab09e76-af60-47d2-94f4-5f67badac07cn%40googlegroups.com.

[Ryan Burner]

Hi Daniel,

No worries, thank you so much for the quick replies, I think I mostly get it now! All we need is a function that returns the correct log(probability) for any given value for any specified covariate. So I can easily write a single function that does that. I do run into some trouble registering the distribution though (below), as I can't find an example of a discrete custom distribution for that function.

On the discrete question though - I have some choice in how narrow I make the bins in the discrete distribution. All of the ~500 values I'll use to define the prior for a given site_covariate fall between 0 and 1, so it might make sense to use a bin width of ~~0.01 or 0.02. Any more precisions is unnecessary. But then I think I should (as below) let the function round the input value to the nearest multiple of that bin width, right, because most of the time that number will be close to a lookup value in my array but not exact, so I wouldn't want to return 0 for every value that had more than 2 decimal places. But then if I'm rounding, do I still need to call it a discrete distribution? 

I'm also curious if any of the ways I wrote the function will be likely to make the model run slowly (it is a beast even with the built-in dbeta() priors.

Thank you again!

Ryan    

#create dummy probability array (real version would be included in model constants)
probArray=array(round(runif(101*1000*5*5),2),dim=c(101,1000,5,5))
rownames(probArray)=c((0:100)/100)

#My original priors.....

for (x in 1:nsites){

  for (y in 1:ncovariates){

    for(z in 1:nyears) {

covs[x, y, z] ~ dbeta(shp1[x, y, z], shp2[x, y, z])

#will instead become something like....

covs[x,y,z] ~ myDiscreteDist(covs[x,y,z], x, y, z, probabilityArray) #I assume that feeding 'covs[x, y, z]' into the function is the way to give it the correct number from the mcmc?

     }

   }

} 

#where....
myDiscreteDist <- nimbleFunction(
  run <- function (c, x, y, z,  probabilityArray ) {
    rowIndex <- match(as.character(round(c,2)), rownames(probArray)) #will parts of this get too slow?
    if (is.na(rowIndex)==FALSE){
      pr <- probArray[rowIndex , x, y, z]
      prLog <- log(pr)
    } else {prLog <- 0}
    return(prLog)
    returnType(double(2))
  }
)

#but then I get stuck here because this doesn't quite work....
registerDistributions(list(
  myDiscreteDist= list(
    BUGSdist = "myDiscreteDist(c, x, y, z, probabilityArray)",
    discrete = TRUE
  )
 )

)

[Daniel Turek]

Ryan, I made some changes to your code.  I'll provide the fixed code below, and a short description of changes I made below.

First, quick note again, please read section 12.2 of the User Manual.  It covers all of this, and more.  It seems you may have followed the part about nimbleFunctions, but not the part about writing custom distributions with nimbleFunctions.

For custom distributions, the name of the density function needs to begin with "d".  I changed yours to "dmyDiscreteDist".

The "run" function needs to be defined using "=", not the assignment operator "<-".  You're not actually defining "run" inside the nimbleFunction() function, but rather, you're providing a function object as the argument named "run" to the nimbleFunction() function.

The returnType is a scalar (a probability, or a log-probability).  So the returnType() statement has to be "returnType(double(0))", to indicate it's returning a scalar.

For a custom distribution, the first argument to the "run" function *must* be called "x" and in the function, this first "x" argument always assumes the value of the LHS (lefthand side) of the stochastic declaration using the custom distribution.  (here, the value of "covs[i,j,k]")  So I added this first argument "x", and changed the other arguments (and your loop indices int eh model code) to i, j, k.

In your model code, the probabilityArray being passed into the distribution requires full indexing.  In model declarations you (almost) always need to provide full indexing on array variables.  So, with full indexing on the array, it needs to look more like:

covs[i,j,k] ~ myDiscreteDist(i, j, k, probabilityArray[1:101, 1:1000, 1:5, 1:5])

In practice, the values 101, 1000, 5, and 5 would be better represented as things like "nsites", "ncovariates", and "nyears", but I don't know which is which, so for the indexing, I used the numbers 101, 1000, 5, and 5.

(side note: without a full model code, etc, I can't actually build / test your model, so some things may not be exactly right).

You can't use "as.character", "rownames", or "match" in the (compilable part of) a nimbleFunction).  So those need to go, and there's no sense in adding rownames to probArray.

You can, however, use the "which" function inside the run function, which I think will help accomplish what you want.

You also need to add the (necessary) "log" argument.  This dictates whether it returns the probability, or the log-probability.  I've added this.

The arguments to the "run" function need type declarations, telling what (scalar/vector/array) type they are.  I added those, as what seems correct.

See changes in the code below, as well as a placeholder for you to "add logic here", which again, should use "x", "probabilityArray", and also the which() function will be handy.

Keep in touch, let me know what doesn't make sense, and take a look at section 12.2 of the manual.

Cheers,
Daniel

#My original priors.....
for(i in 1:nsites){
    for(j in 1:ncovariates){
        for(k in 1:nyears) {
            ####covs[x, y, z] ~ dbeta(shp1[x, y, z], shp2[x, y, z])

            #will instead become something like....

            covs[i,j,k] ~ myDiscreteDist(i, j, k, probabilityArray[1:101, 1:1000, 1:5, 1:5])
        }
    }
}

dmyDiscreteDist <- nimbleFunction(
    ##run <- function (c, x, y, z,  probabilityArray ) {
    run = function(x = double(0), i = double(0), j = double(0), k = double(0), probabilityArray = double(4), log = integer(0, default = 0)) {        
        ######rowIndex <- match(as.character(round(c,2)), rownames(probArray)) #will parts of this get too slow?
        ######if (is.na(rowIndex)==FALSE){
        ######    pr <- probArray[rowIndex , x, y, z]
        ######    prLog <- log(pr)
        ######} else {prLog <- 0}
        ######returnType(double(2))
        ## add logic here, using x, i, j, k, and probabilityArray,
        ## to set "prob" to be either a value from probabilityArray, or 0.
        returnType(double(0))
        if(log) return(log(prob)) else return(prob)

    }
)



#but then I get stuck here because this doesn't quite work....
registerDistributions(list(

    dmyDiscreteDist = list(
        BUGSdist = "dmyDiscreteDist(c, x, y, z, probabilityArray)",

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/0b419780-7be0-496c-8eba-ca2492010018n%40googlegroups.com.

[Ryan Burner]

Hi Daniel,

Ok, I think I've got it! I haven't tested yet but will soon. Note that I will bin the discrete distribution in 101 bins (i.e., 0, 0.01, 0.02.....1.0), hence the 1:101 indexing of the first dimension of the array. 

To lookup the correct probability for x, then, I can simply (when 0 <= x <= 1) round x to two decimal places, multiple by 100, and it becomes the row number of the array. So I avoid e.g., 'which()' that way. Maybe the which() solution, paired with a first column of potential x values, would be more general though but I don't think that mattes in this case. 

Will check back in after testing.

Thanks again,

Ryan

#Priors

for(i in 1:nsites){
  for(j in 1:ncovariates){
    for(k in 1:nyears) {
     

      covs[i,j,k] ~ myDiscreteDist(i, j, k,

                                   probabilityArray[1:101, 1:nsites,
                                                    1:ncovariates, 1:nyears])
    }
  }
}

#Distribution function
dmyDiscreteDist <- nimbleFunction(

  run = function(x = double(0), i = double(0), j = double(0), k = double(0),
                 probabilityArray = double(4), log = integer(0, default = 0)) {        

    if (x>=0 & x<=1) {
      prob <- probabilityArray[(round(x,2)*100)+1,i,j,k]
    } else {
        prob <- 0

      }
   
    returnType(double(0))
    if(log) return(log(prob)) else return(prob)
   
  }
)

#Registering distribution
registerDistributions(list(
  dmyDiscreteDist = list(
    BUGSdist = "dmyDiscreteDist(i, j, k, probabilityArray)",

    discrete = TRUE
  )
 )
)

[Ryan Burner]

Hi Daniel,

Ok, I've almost got this working but am stumped by an error message. Updated code is below, and I've moved some of the indexing outside of the custom distribution function, which seems to speed things up a bit. But whether I keep all the indexing inside the function or do some of it outside I get the same message:

Error in model$checkBasics() :  

Dimension of distribution argument(s) 'probArrayVec' does not match required dimension(s) for the distribution 'dmyDiscreteDistb'. Necessary dimension(s) are 0. 

This is in spite of my 'probVec = double(1)' ( or probArray = double(4) when indexing is inside function). 

Thanks for the help!

Ryan

    ##COVARIATE PRIORS

    #loop over all relevant COVARIATES
    for(mcov in 1:nmiss_cov) {
      #loop over all relevant STOPS
      for (st in stops_miss_any_year[[mcov]]) {

        #define prior distributions by STOP_YEAR
        #loop over all relevant YEARS
        for (yr in yrs_miss_cov_stops[[mcov]][[st]]) {
 
          covs[st,mcov,yr] ~ dmyDiscreteDistb(probVec=
                                                                                  probArray[1:101,
                                                                                                           st,
                                                                                                     mcov,
                                                                                                            yr]
                                           )

        } #years
      } #stops
    } #covariates

####### Create custom distribution ####
#Distribution function
dmyDiscreteDistb <- nimbleFunction(
  run = function(x = double(0), #x is the required name for the left side value
                  probVec  = double(1), log = integer(0, default = 0)) {        

    if (x>=0 & x<=1) {

      prob <-  probVec[round(x*100)+1]

    } else {
      prob <- 0
    }
   
    returnType(double(0))
    if(log) return(log(prob)) else return(prob)
   
  }
)


#Registering distribution
registerDistributions(list(

  dmyDiscreteDistb = list(
    BUGSdist = "dmyDiscreteDistb(probVec)",

    discrete = TRUE
  )
)
)

[Daniel Turek]

Ryan, good question, I think you're only missing one thing.  The distribution definition looks good, but when you manually register the distribution, you also need to tell registerDistributions about any non-scalar arguments. Try adding this one line, to specify the "types" of the arguments:

#Registering distribution
registerDistributions(list(
  dmyDiscreteDistb = list(
    BUGSdist = "dmyDiscreteDistb(probVec)",

    types    = c("value = double(0)", "probVec = double(1)"),

    discrete = TRUE
  )
)
)

And hopefully that should get it working.  Otherwise, nice job with the custom discrete prior.

Cheers,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/316a6817-0a1d-4cc9-904b-65121d4d1d7an%40googlegroups.com.

[Ryan Burner]

Hi Daniel,

That was it, thank you, works now! 

Next challenge though - I'm getting these warning messages below during runMCMC(), and the vast majority of my covariate chains are just stuck on their initial values, although a few chains run off into enormous values (~1 million, which should all be zero probability) so still some issues. 

Warning: slice sampler reached maximum number of contractions for 'covs[43, 3, 1]'. Current parameter value is 0.45.

Maybe something with my initial values? I think I'm setting all initial values as the most probable value for each covariate, but it is possible I'm messing that up. Will hopefully have some time to sort that out tomorrow and will report back if it fixes anything. These discrete priors can make a pretty narrow target - often all of the prior probability will fall into a range of values from e.g. 0.3 - 0.5. It seems like most of the 'standard' distributions would have long tails of gradually changing probability to 'guide' the mcmc to the most likely values, whereas I'm lacking that. I could artificially create it of course if that was something people do.  I don't know how the mcmc jump size (not the right term I'm sure) adapts to this kind of situation....   

Thanks again!

Ryan

[Daniel Turek]

Ryan, are you certain you're using valid initial values for the model, prior to beginning the MCMC?  For example, if you run model$calculate(), does that return a real (usually negative) number?  If, instead, it returns NaN, NA, or -Inf, that could be the root of your problem, being that the model is not properly initialized, and the MCMC never has a valid foothold (location in parameter space) from which to begin exploration.

If the model is properly initialized (model$calculate() returns a real number), then there could potentially be an issue with the implementation of your custom discrete prior, which is not "playing nicely" with the slice sampler, which I assume is assigned to update those nodes.  That's something we could take a closer look at, if necessary, but the initialization is the first thing to check.

Cheers,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/29144499-5c0e-4829-a971-613aba3dbed1n%40googlegroups.com.