How to perform basic operations in rstan?

Hey all,
for a simple Multinomial-Dirichlet model I want to normalize data, but iteratively (i.e. in every iteration). However, I am already failing to tell R/Rstan data / constant appropriately. Since the multinomial_lpmf is defined for int’s only, I want to use the int operators descibes here: 2.1 Integer-valued arithmetic operators | Stan Functions Reference.
R however complains immediately that “operator” is unknown. In summary: How to apply simple arithmetic operations?

``````// The input data is a vector 'y' of length 'N'.
data {
int<lower=0> N;
int<lower=0> G;
int counts[N,G];
vector<lower=0>[G] alpha;
real r;
real Tstart;
}

// The parameters accepted by the model. Our model
parameters {
simplex[G] theta;
}
transformed parameters {
real T = Tstart;
T = r*T;
}

// the model to be estimated
model {
for(i in 1:N){

target += dirichlet_lpdf(theta | alpha);
target += multinomial_lpmf(int operator%/%(counts[i], T) | theta);
//void print (T);
}
}
``````

I have attached a csv that contains the data for counts. r and Tstart shall define a geometric sequence with T[j] = r^j * Tstart (or less costly by saying T = r * T in every iteration).

Bonus question: What about void print (T)? Can I print something within the execution of compiled stan files?
rna_data.csv (1.2 KB)

BG Jonas

The `int operator%/%` notation in the documentation is just to indicate that `%/%` is an operator, rather than a function to be called. In other words, you just want to use `x %/% y`.

As for the model, the integer division operator is only compatible when both inputs are integers:

``````int %/% int
``````

(i.e., type `int`), however in your model you’re attempting to call:

``````int[] %/% real
``````

Where `counts[i]` is one-dimensional array of integers (`int[]`), and `T` is a real.

Also, I’m not entirely sure what you’re attempting when you refer to ‘normalising at each iteration’. It looks like at every iteration, you want to divide the integer responses for each observation (i.e., `counts[i]`), by `Tstart * r`). But none of these quantities are parameters, and so would not change at all during sampling. It would likely provide a pretty good speed-up of your model if you just pre-computed this division and passed it as data

First of all, thanks for the loightning-fast answer.

Regarding `x %/% y` … I thought so. But even with `counts[i] %/% 10` I get an “PARSER EXPECTED: <expression>” error. If it was typical normalization, you were of course right. But the idea is simulated anealing with a decrasing normalization to mitigate steep gradients (not for this particular data set though). Hence, I could round the T to an integer if the issue above is clarified.

When you call `counts[i] %/% 10` you’re still not using a valid signature, unfortunately the `%/%` operator is not vectorised, so you’ll have to loop over each integer in `counts[i]`:

``````model {
target += dirichlet_lpdf(theta | alpha);

for(i in 1:N) {
int counts_normalised[G];
for(g in 1:G) {
counts_normalised[g] = counts[i, g] %/% 10;
}
target += multinomial_lpmf(counts_normalised | theta);
}
}
``````

Note that I’ve also moved the `dirichlet_lpdf` outside of the loop, as including it within the loop is essentially weighting the posterior `N` times.

Just to also quickly clarify, the key thing here is that `T` is constant throughout all iterations (in this model at least), so there is no benefit to performing the division in the `transformed parameters`/`parameters` blocks.

If this more of a dummy model, and you’re intending to construct `T` using something declared in the `parameters` block, then you’re likely to run into a different set of issues. Taking a continuous parameter and discretising it to return an integer would result in discontinuities in the respective parameters, and result in seriously poor sampling performance

Forgot to include a reference for the issues with sampling! There’s some more information in this section of the manual: 3.7 Step-like functions | Stan Functions Reference

Okay, this helps a lot! I thought it is vectorized. Is there any other option to avoid the loop?
E.g. `normalized = floor(count ./ T)` and defining `count` as `matrix` before?
This way temperature could be a `real` defined in `parameters` and I do not have to discretize it (referring to your warning; here count and normalized represent data and I would say it does not affect the gradients/performance). Below you see what I got so far (btw I still get the error `PARSER EXPECTED: <expression>`).

One last question until I would happily close this topic: How can I create an iterative (more precisely geometric) parameter adjustment? For instance, can I initialize `real temperature = Tstart` and then declare in every iteration something like `temperature *= r`? If not, is there an accessible iteration variable, e.g. `i`, that counts how many samples have been drawn? It is so embarassing to ask this, but I am not even able to define `foo += 1` myself in stan.

``````data {
int<lower=0> N;
int<lower=0> G;
int counts[N,G];
vector<lower=0>[G] alpha;
real<lower=0, upper=1> r;
int<lower=1> Tstart;
}

// The parameters accepted by the model.
parameters {
simplex[G] theta;
}

transformed parameters {
//matrix[N,G] normalized;
//normalized = to_matrix(counts) / T;
real temperature = r * Tstart;
}

// the model to be estimated
model {

for(i in 1:N) {
int normalized[G];
for(g in 1:G) {
normalized[g] = counts[i,g] %/% 1;
}
target += multinomial_lpmf(normalized | theta);
}
target += dirichlet_lpdf(theta | alpha);
}
``````

Unfortunately there is no way to access the iteration number within a Stan model.

For your model, the simplest approach to constructing the normalised integer outcomes is to perform the `count %/% (Tstart * r) ` division in R and pass the result as an additional argument to your stan model

Thanks for all your support. I will leave it like that.
Just to make this clear: For the model gescribed above, normalizing beforehand might work. But I require iteration-specific normalization. Your suggestion would mean:

1. Compile the stan code
2. Normalize in R
3. Run one iteration using rstan
4. Do 2. and 3. roughly 10.000 times
``````real geometric_rng(real temperature,real r){
and apply it iteratively on the parameter `temperature` in `generated quantities`.