Integer floor calculation or type casting

[N Sanders]

Hello,

I've come across a simple procedure that I can't figure out how to express in Stan code.  

I need to access a value in an array ("A") at an index ("i") which can be approximated from the current value of a real parameter ("p") in the chain. I had planned to do this using the statement "i <- floor(p)" and then accessing "A[i]".  However, I have found that the floor function in Stan returns a real value, which cannot be used as an array index, and I see from the manual that real values cannot be demoted to integers in Stan code.

Is there another way I can accomplish this task?

(P.S. I regret that I haven't converted my model to Stan v2 yet, so I'm referring to my experience with v1.3.0)

Thank you in advance for any advice!

Sincerely,

~Nathan

[Bob Carpenter]

On 11/7/13, 12:35 AM, N Sanders wrote:
> Hello,
>
> I've come across a simple procedure that I can't figure out how to express in Stan code.
>
> I need to access a value in an array ("A") at an index ("i") which can be approximated from the current value of a real
> parameter ("p") in the chain. I had planned to do this using the statement "i <- floor(p)" and then accessing "A[i]".
> However, I have found that the floor function in Stan returns a real value, which cannot be used as an array index,
> and I see from the manual that real values cannot be demoted to integers in Stan code.
>
> Is there another way I can accomplish this task?

Not as of yet. At least not cleanly and efficiently.

You can do integer floors for indexes as follows and
it wont' be too inefficient for small values of p:

int k;
k <- 1;
while ((k + 1) < floor(p))
k <- k + 1;

Then when you're out of the loop, k = floor(p).

BUT, there's an underlying problem with converting real
values to integers in Stan, which is that it loses gradient
information. This will make sampling very inefficient.

An alternative would be to directly code using if-then if
it's a mixture model. Of course, that also loses gradient
information on p.

At the very least, you'll need to constraint p so that
the posterior is proper.

> (P.S. I regret that I haven't converted my model to Stan v2 yet, so I'm referring to my experience with v1.3.0)

The modeling language hasn't changed, only the underlying
implementation.

- Bob

[N Sanders]

Hi Bob,

Thanks very much for your reply! Your suggestion was perfect - finding the floor with a loop works well and doesn't seem to noticeably slow down sampling.

For my particular case, I think the floor calculation should not have a significant effect on the gradient calculation.  The index is just used to look up a value from a table which is applied as a small correction to the model value being compared to the data in the likelihood.  I use the floor loop to find the right index in the table, but then I linearly interpolate between the table values, so I think the gradient should still be smooth.  Also, the correction is just a local variable in the model section, so no sampled parameters are being converted to integers.  

Thanks again for your help!

Best,

~Nathan

[Michael Betancourt]

For my particular case, I think the floor calculation should not have a significant effect on the gradient calculation.  The index is just used to look up a value from a table which is applied as a small correction to the model value being compared to the data in the likelihood.  I use the floor loop to find the right index in the table, but then I linearly interpolate between the table values, so I think the gradient should still be smooth.  Also, the correction is just a local variable in the model section, so no sampled parameters are being converted to integers.  

Your correction may very well be smooth, but calculating it this way lose the gradient information of the correction so even if the posterior is calculated correction its gradient won't.  In the end the Metropolis/slice procedure will correct for any error,

but whenever the correction becomes important then the error will induce small acceptance probabilities and poor sampling.  That's the biggest concern.

[Marcus Brubaker]

Actually, the gradient is probably fine, note that he said he was doing linear interpolation.  My guess is Nathan is doing something like:

xm <- x[floor(c)];

xp <- x[ceil(c)];

x <- (c - floor(c))*xp + (ceil(c) - c)*xm;

In which case the derivative of x with respect to c is reasonably fine.  The second derivative would be wrong (always = 0) but that's another issue.

Cheers,

Marcus

--
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.

[N Sanders]

Hello again,

Sorry to bump this old thread, but it seems most relevant to a new question I have:

I am trying to take the mean value of an integer array (basically finding the expectation value of a binary vector) and am not sure how to do it.  When simply using the mean function, the parser throws the error:

no matches for function name="mean"
    arg 0 type=int[1]
available function signatures for mean:
0.  mean(real[1]) : real
1.  mean(vector) : real
2.  mean(row vector) : real
3.  mean(matrix) : real
no matches for function name="mean"
    arg 0 type=int[1]
available function signatures for mean:
0.  mean(real[1]) : real
1.  mean(vector) : real
2.  mean(row vector) : real
3.  mean(matrix) : real
expression is ill formed

This seems to be another type casting issue.  The manual states "Stan automatically promotes integer values to real values if necessary," but this does not seem to hold in this case. 

One approach I have considered is to calculate the mean manually, by summing the integers and dividing by the vector length, but I think that will lead to integer rather than float division.

This was meant to go in the transformed data block of my model, so for now I can just precompute the mean and pass it as data.  But perhaps there is a better way to do this directly in Stan.

Thanks for any advice!

Best,

~Nathan

[Daniel Lee]

One way to get around it is to wrap the int[] with a call to to_vector() in the call to mean. Want to give that a try and see if it works?



Daniel

PS. int gets promoted to real, but int[] doesn't get promoted to real[].

--
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/d/optout.

[Bob Carpenter]

Ideally, int[] would be promoted to real[] everywhere, but
we aren't quite there yet. That would create the function

real mean(int[]);

This is because we promote arguments, but we do not demote return
results. This is the right behavior here, too, because the
sample mean of a collection of ints is not necessarily an int.

- Bob

[Nathan Sanders]

Hi Daniel and Bob,

Daniel, yes, this does work - thank you! 

Bob, it sounds like there are some cases where real promotion would not always be the expected result for mean(int[]).  In that case, perhaps there could be added an explicit function for promoting integers or a small discussion of this workaround (i.e. using to_vector for the purpose of promotion to real) in the manual.

Thanks to you both for offering advice!

~Nathan

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/uqh1h158ntg/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

I just created an issue for completing the int[] to real[] promotions
in argument position.

https://github.com/stan-dev/stan/issues/830

This will not change any return types. In some cases, we have
specialized functions for int array inputs. For example, max
has two signatures:

(int[]):int
(real[]):real

so we wouldn't promote the int[] to real[] and return real,
but stick with the specialized signature that's already there.

I think I'll hold off on updating the manual. Adding the promotions
shouldn't be too difficult and there would be a lot of places to
put that hint for now.

[Nathan Sanders]

Hi Bob,

Excellent, thank you for your response to this!

Best,

~Nathan

[Daniel Rosenbloom]

Hi Nathan--

I wonder if you, or any others following this thread, figured out a faster way to cast reals to ints. I have an analogous problem that I want to solve by linear interpolation (my likelihood function is slow to compute, and so I pre-computed a lookup table).

It appears that I can successfully do inference with my interpolated likelihood, with either optimize or sampling (acceptance probs > 0.9, so it seemed to use the gradient information correctly), using the hack that Bob suggests. Does anyone know of a more direct hack to typecast without searching through a list, to speed things up? (Failing that, I could try binary search, or cubic interpolation on a smaller lookup table.)

Marcus pointed out that the second derivative is always zero -- would this end up posing any problems?

Thanks,

Daniel

[Bob Carpenter]

We're taking the nanny software approach on this one
and not implementing it on purpose because it's dangerous.
In general, it can break the gradients and therefore break
all of our gradient-based algorithms.

Is there a particular form of interpolation you need
that we might be able to build in and guarantee that
the gradients work? I know there was one built into BUGS.

Binary search is the best you'll be able to do now.

- 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.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Daniel Rosenbloom]

Hi Bob,

Thanks for clarifying. Binary search is plenty fast -- I get 70k samples per second using it, with a pretty large table (22000 rows). I'm attaching my pystan code, so that you can see how we are doing the interpolation. The interpolation is linear in x (for x<10) or in log(x) (for x>10); this was chosen since the distribution scales as LL=-2*log(x) for large x. Left and right tail behavior is also given.

The lookup table (landtable.py) is for the Landau distribution, computed using Mathematica, which has the distribution built-in (the standard form of the distribution is LandauDistribution[0,Pi/2] in Mathematica). Note that while there is an analytic approximation to this distribution (Moyal approximation, described at https://en.wikipedia.org/wiki/Landau_distribution), it is very bad outside a small interval, and it misses the long tail behavior.

Is it OK to have flat second derivative? This linear approach samples quite accurately, and I'm getting no Metropolis rejection warnings.

Thanks again!

Daniel

[Bob Carpenter]

Flat second derivative of the log density? The second
derivative w.r.t. y of log normal(y | 0, 1) is constant.

We'll need Michael's help for real geometry questions,
so may need to put this discussion in a different thread.

And thanks for posting the model; it's great to see this
kind of careful numerical coding for real problems.

- Bob

> <landau_stan.py><landtable.py>

[Glenn Goldsmith]

Very late to this conversation, but just in case it's useful to anyone else, you could #include a custom c++ function that executes the cast - something like:

template <typename T0__>

int

real_to_int_cast(const T0__& x, std::ostream* pstream__) {

    return int (x);

};

This is pretty easy in rstan, using the instructions here: https://cran.r-project.org/web/packages/rstan/vignettes/external.html. I assume it should also be easy in cmdstan, though I can't actually get it to compile property myself. Not sure about pystan.

[Ben Goodrich]

On Friday, April 21, 2017 at 1:01:16 AM UTC-4, Glenn Goldsmith wrote:

Very late to this conversation, but just in case it's useful to anyone else, you could #include a custom c++ function that executes the cast - something like:

template <typename T0__>

int

real_to_int_cast(const T0__& x, std::ostream* pstream__) {

    return int (x);

};

This is pretty easy in rstan, using the instructions here: https://cran.r-project.org/web/packages/rstan/vignettes/external.html. I assume it should also be easy in cmdstan, though I can't actually get it to compile property myself. Not sure about pystan.

For CmdStan, real_to_int() needs to be in the stan::math namespace. There is no way to do this with PyStan currently, but there is a GitHub issue about it.