How to define function in OpenBUGS

Question

In my Hierarchical Bayesian Model, in the data layer y[i] <- M(t,u,v)+ N[0,sigma_y] where M is a complicated function. I have generated the output of this M by some numerical simulation for a large grid of parameters (t,u,v).

In order to do MCMC on this grid i need to define a function which looks up the table of the Model M's output for a given input (t,u,v).

But in OpenBUGS documentation, I could find only standard distributions like dnorm etc for defining y[i].

Isn't it possible to use a general function for it? Or, Is there any other way to do it?

Update: My model is as shown below.

$i$ represents the measured data, and $j$ represent the groups.

$y_i \sim M(t,u_{j[i]},v_{j[i]}) + \mathcal{N}[0,\sigma_i]$

In the group level, $(u_j,v_j) \sim \mathcal{N}[(\mu_{u},\mu_v), \Sigma]$

The function $M$ does not have a simple analytical form to express in Bugs. But if i can define a general function to look up a table, then for each step in MC chain, I can lookup the value of deterministic relation $M$ for a given t,u,v parameter.

My goal is to estimate, $\mu_u$, $\mu_v$,$\Sigma$ and $t$ which best fits the the given $y_i$ and $\sigma_i$ data.

I am attaching a diagram of the model below.

Answer 1

In your question you are assigning a value to y[i], not describing how the y's are distributed. If I understand you correctly, you want to describe the deterministic/predicted value Yhat using the function M(t,u,v), with the stochastic/sampling process determined by a normal distribution (i.e. the actual observations). This can be accomplished by:

 y[i] ~ dnorm(Yhat[i], tau.sigma_y)
 Yhat[i] <- M(t,u,v)

In more general terms:

$$y_i \sim \mathcal{N}(\hat{y}_i,\sigma^2_y) \\ \hat{y}_i = M(t,u,v)$$ You then need to feed the correct values from the function M to each observation (through Yhat), but that can be accomplished in several different ways. I don't know the specifics of your function, but as an made-up example case a function M(t,u,v) could be written as:

Yhat[i] <- exp(beta1*t[i]) - (beta2 * u[i] * abs(v[i]))/(beta3 + u[i])

Here values of u, v and t for observation i is fed to Bugs as data and the parameters beta1, beta2 and beta3 are estimated (and you need to supply initial values of these). Or do you have a function M that truly cannot be specified in Bugs?

You can also "pre-calculate" the values of M for each observation (from the t,v & u of the observation) and supply this as data to the model as:

Yhat[i] <- M[i] # or skipping this step since it's redundant 

This means that Yhat[i] is directly calculated from your simulation, and the step is redundant. You could just as well give Yhat as data (directly in Openbugs or e.g. called from R through R2Openbugs). Without knowing the exact purpose of your modelling exercise it is difficult to know what you want to accomplish with the MCMC model.

If you need more flexibility to create user-defined functions, maybe to make external calls, you should look at STAN, where I think this is possible to a larger extent (I haven't used STAN though).