Errors from rinvgamma

[DRJP]

I am getting error messages from a minimal BUGS model when using dinvgamma. I get no such messages when I use dgamma instead.

With dgamma...

test1 <- nimbleCode ({
    hyperSh <- hyperShape
    hyperSc <- hyperScale
    Shp      ~ dgamma(shape=hyperSh, scale=hyperSc) 
})
Inits <- list(hyperShape=1.0, hyperScale=1.0, Shp=1.0)
mod1  <- nimbleModel(test1, inits=Inits)
simNodes <- mod1$getDependencies(c("hyperShape","hyperScale"), downstream=TRUE, self=FALSE)
mod1$simulate(nodes=simNodes)   ## Works as expected
mod1$calculate()                              ## Works as expected

and with dinvgamma

test4 <- nimbleCode ({
    hyperSh <- hyperShape
    hyperSc <- hyperScale
    Shp      ~ dinvgamma(shape=hyperSh, scale=hyperSc) 
})
Inits <- list(hyperShape=1.0, hyperScale=1.0, Shp=1.0)
mod4  <- nimbleModel(test4, inits=Inits)
simNodes <- mod4$getDependencies(c("hyperShape","hyperScale"), downstream=TRUE, self=FALSE)
mod4$simulate(nodes=simNodes)
## Error in rinvgamma(1, shape = model$hyperSh[1], rate = model$lifted_d1_over_hyperSc[1]) : unused argument (rate = model$lifted_d1_over_hyperSc[1])

Swapping from using scale to using rate has no effect.

I'm using the latest CRAN version of nimble on Ubuntu.

[danielturek]

@DRJP Only noting that your minimally reproducible example (using rinvgamma) runs fine for me (no errors, and output appears correct), using nimble version 0.6-8 (installed from CRAN), running on OSX. I took a quick look at the codebase, and the error you're reporting, as best I can tell, doesn't seem like it should happen. So I'm at a bit of a loss.

@paciorek Able to try this out in a linux environment?

[DRJP]

Hi Daniel. Thank you for the rapid reply. If I write a wrapper function for the distribution I can bypass the error.

dInvGamma <- nimbleFunction (
    run = function(x     = double(0),
                   Shape = double(0, default=1), 
                   Scale = double(0, default=1), 
                   log   = integer(0, default = 0)) {
        returnType(double(0))
        logDensX = dgamma(x, shape=Shape, scale=Scale, log=TRUE) 
        if (log) return(logDensX) else return(exp(logDensX)) 
    }
)

rInvGamma <- nimbleFunction (
    run = function(n = integer(0, default=1),
                   Shape = double(0, default=1), 
                   Scale = double(0, default=1)) {
        returnType(double(0))
        if (n!=1)
            stop("Only implimented for n=1")
        x = rgamma(n=1, shape=Shape, scale=Scale)    
        return(x)
    }
)

registerDistributions(list(dInvGamma = list(
                               BUGSdist = "dInvGamma(Shape,Scale)",
                               discrete = FALSE,
                               range    = c(0, Inf), 
                               pqAvail  = FALSE)))

test5 <- nimbleCode ({
    hyperSh <- hyperShape
    hyperSc <- hyperScale
    Shp      ~ dInvGamma(Shape=hyperShape, Scale=hyperScale) 
})
Inits <- list(hyperShape=1.0, hyperScale=1.0, Shp=1.0)
mod5  <- nimbleModel(test5, inits=Inits)
mod5$plotGraph()
(simNodes <- mod5$getDependencies(c("hyperShape","hyperScale"), downstream=TRUE, self=FALSE)) ##
mod5$simulate(nodes=simNodes)
mod5$calculate(nodes=simNodes)
mod5$Shp

This appears to run fine.

[DRJP]

Here's my session info...

> sessionInfo()
R version 3.2.3 (2015-12-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.3 LTS

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_GB.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_GB.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] igraph_1.1.2   MCMCpack_1.4-0 MASS_7.3-45    coda_0.19-1    nimble_0.6-8  

loaded via a namespace (and not attached):
 [1] lattice_0.20-33    codetools_0.2-14   grid_3.2.3         R6_2.2.2          
 [5] MatrixModels_0.4-1 magrittr_1.5       SparseM_1.77       Matrix_1.2-3      
 [9] tools_3.2.3        parallel_3.2.3     compiler_3.2.3     pkgconfig_2.0.1   
[13] mcmc_0.9-5         quantreg_5.33  

[DRJP]

This seems to be a conflict with package MCMCpack. In a new session things work fine if I don't load that package.

Incidentally, I only use MCMCpack now because of the Nimble message "rdirch only handles n = 1 at the moment".

[danielturek]

@DRJP Nice detective work, thanks for following-up.

[JimmyJHickey]

I had the same issue with the LaplacesDemon package which also has a dinvgamma function. After detaching the package with detach("package:LaplacesDemon") I was able to use the dinvgamma function within the nimble code.

Is there any way to specify that the nimble function should be used rather than having to detach other packages? When using the standard R syntax for this: Shp ~ nimble::dinvgamma(shape=hyperSh, scale=hyperSc), I get an error Error in if (dist %in% c("T", "I")) { : the condition has length > 1.

[paciorek]

Thanks, this is helpful to see this case - we're not set up to handle the namespace resolution operation.

I'll take a look to see how we can better handle internally to make sure we are using nimble's version of a density even if others are available in a given session.

[paciorek]

So more generally, we are prone to masking of any R functions used in nodeFunctions, including calculate, simulate, getParam, etc. during uncompiled execution, and similarly in uncompiled execution of nimbleFunctions. E.g., if a user creates their own sqrt, that would affect any use of getParam that uses sqrt and any use of sqrt in a nimbleFunction.

This is because the model and node functions (and user-defined nimbleFunctions) are in the user's global environment and not the nimble namespace. So R looks up through the search path, looking in the global environment first and then in package namespaces where the order of lookup depends on the order a user loaded packages. This is obviously prone to errors/lack of reproducibility.

If we wanted to be really safe we would include namespace resolution (generally stats:: but sometimes nimble::) for any R function in the various nodeFunctions and user-defined nimbleFunctions and then strip it out during compilation.

We could also plan to address this in nCompiler/new nimbleModel and leave as is for now.

Or we might just consider this just how nimble works given R's scoping/lookup rules. But I don't think that feels right to me given how we define what it means to use the nimble language -- we have specific C++ translations for non-user-defined functions invoked by a user and therefore probably want that to the case in R too.

@perrydv @danielturek Thoughts?

[paciorek]

Also looping @kenkellner in on this given that namespace resolution came up in the context of macros.

One other note here is the perhaps obvious comment that we do by default do model$calculate uncompiled when building a model.

[perrydv]

@paciorek I agree with the idea of deferring this to the deeper renovations we have underway. It seems that if we have support for namespace resolution as a goal from the start we will be able to build it in comprehensively. At the moment although it is obviously clunky, it is not causing a large number of problems and there are workarounds.

[paciorek]

The one caveat is how we want to handle the question of scoping in terms of macros.

Perhaps @kenkellner can open a new NCT issue that reminds us of the need in terms of scoping/namespace/packaging questions.