Implementing complicated models (defining custom functions?)

[Chris Rose]

Hi

I am working on a problem where the model with which I'm trying perform inference is very complicated. In brief, the model describes the spatial variation of a measurable quantity over three (though possibly more) dimensions, and the model itself can be though of as being piecewise, with very many pieces. While I think Stan has the necessary functions and control structures to implement the model, the model file would be quite horrific to read. That is, of course, unless it's possible to define functions within the Stan language; from reading the manual, it doesn't seem like this is possible (am I wrong?).

What options are available to me if I want to use Stan? Can I define my own functions within the Stan language? Can I do so (sanely) by writing some kind of extension in C++?

Thanks in advance,

— Chris

[Bob Carpenter]

There's been a lot of discussion on this list recently about
what you need to do to interface a C++ function to Stan.

Short answer: If you can write a C++ function where parameters
are instantiated as stan::agrad::var and data as double/int
depending on type, then you can interface pretty easily.

This means all routines acting on parameters need to also be
defined for stan::agrad::var. You can use our stan::math and
stan::agrad and stan::prob libs, but it's usually a lot of work.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.

[Chris Rose]

Thanks Bob. I searched the mailing list before posting, but perhaps I found material that was too in-depth and dismissed it. Are you able to provide a link or two to relevant discussions?

Is it possible or likely that the Stan language will develop to allow users to define their own functions? It seem the more procedural nature of the Stan language versus the more declarative BUGS/JAGS language would lend itself to this. Or, would this encourage users to formulate models that are far too complex to reasonably fit?

— Chris

[Bob Carpenter]

I interpreted you to mean defining your own functions in C++
and then interfacing them with Stan. I'd start here if that's
your goal: https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan

If you mean do we have plans to add a way to define functions within
the Stan language, yes, that's hopefully something we'll roll out this
year. The problem's going to be that it won't be much use if the
functions can't be independently tested / debugged.

- Bob

> stan-users+...@googlegroups.com <javascript:>.
> > For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.

[Marcus Brubaker]

Hi Chris,

The short answer is that no, you can't define custom functions easily in Stan.  If you're willing to hack the C++ (as Bob alluded to) you can add your own functions that way.  At a really high level, how to do this is described here: https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan  If you only wanted the function for yourself (i.e., not planning to submit it to the main repository), you can skip all the testing.

What you seemed mostly to be asking about was a way, in a Stan model file, to define a new function in the Stan modelling language.  This is not possible right now but is something we've talked about adding since Stan was first released.

That said, if you're willing to use the C preprocessor yourself, there is a way you can do this in your own.  Basically, use C preprocessor macros in a file, lets say called foo.stancp

#define FOOFUNC(a,sigma) a ~ normal(0,sigma)

model {

  FOOFUNC(x,sigma1);

  FOOFUNC(y,sigma1+sigma2);

...

}

and then, when you're building the model using CmdStan, call:

> cpp foo.stancp foo.stan

> stanc foo.stan 

Where the first call uses the C preprocessor on it's own.  If you're using RStan or PyStan this is still possible, but would be a little more complicated.

Cheers,

Marcus

[Chris Rose]

Thanks, Bob, that's useful. And it would be great to see user-definable functions (and some sensible way to test and debug them) in Stan.

— Chris

[Chris Rose]

Thanks Marcus. I did think of using some kind of macro substitution, but completely overlooked the use of the C/C++ preprocessor for this — neat idea.

— Chris

[Bob Carpenter]

On 1/23/14, 5:37 PM, Chris Rose wrote:
> Thanks, Bob, that's useful. And it would be great to see user-definable functions (and some sensible way to test and
> debug them) in Stan.

User-definable functions are on the short-term to-do list.

I'm imagining a syntax like this for functions:

simplex foo(vector x, real[,] y, int N) {
... local var defs ...
... statements ...
return result;
}

for a function that takes a vector, a 2D array of reals,
and an integer as input and returns a simplex.

And then for relational calls that could be used to
define submodels, something like:

void linear_regression(matrix x, vector beta, real<lower=0> sigma) {
sigma ~ cauchy(0,2.5);
beta ~ cauchy(0,2.5);
y ~ normal(x * beta, sigma);
}

So you'd still have access to the log probabilities where
appropriate.

We still have to work out the details.

- Bob

[Chris Rose]

One comment on your proposed syntax for defining functions:

In general, I like the option to use named function arguments rather than argument order carrying this semantic information. This is something that I think R does well, particularly in the way it allows optional arguments to be defined with default values. I would not argue that argument names must be specified when calling functions, and that users who prefer using argument order to specify the meaning of the arguments could do so. I find named arguments make code more readable and lead to fewer bugs (e.g., no confusing standard deviation for precision when specifying a normal, for example). I note that named parameters is already used in Stan for specifying limits on values.

So, to adopt your example, what about the following:

// Define functions as you suggested:

void linear_regression(matrix x, vector beta, real<lower=0> sigma) { 
      sigma ~ cauchy(0,2.5); 
      beta ~ cauchy(0,2.5); 
      y ~ normal(x * beta, sigma); 
    }

// Optionally call functions using named arguments:

linear_regression(beta = 0, x = my_matrix, sigma = 0.5);

Best,

Chris