# Analysis API

**URL:** https://discourse.mc-stan.org/t/analysis-api/3486
**Category:** Developers
**Created:** [March 6, 2018, 3:04pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486 "2018-03-06T15:04:44Z")
**Posts on this page:** 20
**Page:** 1

<div class="post-metadata">

### Author: ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)
#### Post date: [March 6, 2018, 3:04pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486/1 "2018-03-06T15:04:44Z")

</div>

In order to avoid the extensive code duplication that is currently spread across the interfaces we need the core Stan API to start providing analysis routes, including the calculation of effective sample size, rhat, divergences, E-FMI, etc. This will require some reorganization of the `stan` source.

First and foremost we need to decide on the exact format for I/O, in particular how sample output must be input to the API and how messages and error status is returned by the API. Much of this can be cribbed from the recent discussion about exposing C++ effective sample size functions .

Then we need to reorganize the `stan` source. Namely we need all analysis code collected together (in particular, effective sample size and rhat functionality should be removed from the chains class as the OO structure is not needed). The combined functionality then needs to be exposed via appropriate services.

Given the partial steps that are being discussed we might as well do this systematically.

---

<div class="post-metadata">

### Author: ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)
#### Post date: [March 6, 2018, 6:48pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486/2 "2018-03-06T18:48:59Z")

</div>

I’ll add non-MCMC analysis to this system. Now we have only ADVI, but we may have something else later, but all distributional approximations can use the new Pareto diagnostic. Given draws theta^{(s)} from any distributional approximation q, and log(q(theta^{(s)})) and log(p(theta^{(s)})) we can compute Pareto k (khat) and effective sample size (N\_eff) estimates. These should be reported in every interface similarly as Rhat and N\_eff are reported for HMC/NUTS.

---

<div class="post-metadata">

### Author: ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)
#### Post date: [March 6, 2018, 9:20pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486/3 "2018-03-06T21:20:41Z")

</div>

Yes, once we’ve established conventions for the routes we can add others for other algorithms. We will, however, have to be careful to either enforce the necessary inputs or have well-defined error messages (and return codes) when the inputs are not well-formed for a particular analysis.

---

<div class="post-metadata">

### Author: ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [March 7, 2018, 10:37pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486/4 "2018-03-07T22:37:18Z")

</div>

I think that’s a good idea. I also think the core implementations could go into stan-dev/math, with only the interfaces and I/O dealt with at the Stan level.

I’m most concerned with getting the C++ APIs in better shape for all these functions. We’ve been discussing the effective sample size calcs in another thread and thought `vector<const double*>` would be general enough to avoid big copies and accomodate both `std::vector` array types and `Eigen::Matrix` vector types.

Currently, I/O isn’t handled at the Stan level—it’s just delegated to the interfaces. RStan has the ability to write out output that CmdStan’s analysis tools can read and also write out data that CmdStan can read. It can also read CmdStan data since that’s just an R format. Are you going to be suggesting restarting the JSON/protobuf discussions?

---

<div class="post-metadata">

### Author: ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)
#### Post date: [March 7, 2018, 11:45pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486/5 "2018-03-07T23:45:03Z")

</div>

No, I am suggesting that we define an I/O convention for how the interfaces have to wrap the samples and diagnostic info from a fit object into std::vectors or POD arrays to be analyzed by routes and then what they expect in return (probably strings and error codes if not numerical results like ESS and rhat).  
Bonus points if the convention minimized excess copying on the C++ side.

None of this would be exposed to users – the interfaces would call these routes internally.

I am hesitant to put this stuff into math for the same reason I don’t think the Welford stuff really belongs there – the math library is essentially an autodiff library, not a general purpose math library. This stuff would be just as well in a separate auxiliary math library.

---

<div class="post-metadata">

### Author: ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [March 8, 2018, 1:03am UTC](https://discourse.mc-stan.org/t/analysis-api/3486/6 "2018-03-08T01:03:23Z")

</div>

> [@betanalpha](#):
>
> the math library is essentially an autodiff library, not a general purpose math library

That’s a nice way to think about it. Effective sample size estimators and R-hat are more like the sampling algorithms—subject to change as we learn more. So it makes sense to leave those with the algorithms. Plus, I think it makes sense to have you manage them.

If we ever get around to adding derivatives to FFT, that’ll go in the math lib.

Welford’s algorithm is one of those oddball things where it’s at an API level of an accumulator rather than being like our ordinary math functions. So I’m OK with that going either way. I could see some of the algorithms that do need to be in the math library, like calculating sample variance, needing it, though. If that happens, then whatever they use in common should be in the math lib.

> [@betanalpha](#):
>
> I am suggesting that we define an I/O convention for how the interfaces have to wrap the samples and diagnostic info from a fit object into std::vectors or POD arrays to be analyzed by routes and then what they expect in retur

OK, that’s what we’re already doing. See [this pull request for n\_eff](https://github.com/stan-dev/stan/pull/2477).

I’ll reverse link that back to here, too.

---

<div class="post-metadata">

### Author: ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)
#### Post date: [March 8, 2018, 2:30am UTC](https://discourse.mc-stan.org/t/analysis-api/3486/7 "2018-03-08T02:30:27Z")

</div>

> [@Bob\_Carpenter](#):
>
> OK, that’s what we’re already doing. See this pull request for n\_eff.

Yes, this is what incited my post. We shouldn’t do this piecemeal but put something systematic together on which we can add all of the other diagnostic calculations.

---

<div class="post-metadata">

### Author: ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [March 10, 2018, 12:56am UTC](https://discourse.mc-stan.org/t/analysis-api/3486/8 "2018-03-10T00:56:58Z")

</div>

Do you think we won’t be able to use `vector<double*>` everywhere and wrap in `Map<Matrix<...>>` where necessary for Eigen?

---

<div class="post-metadata">

### Author: ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)
#### Post date: [March 11, 2018, 9:23pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486/9 "2018-03-11T21:23:06Z")

</div>

But how would we get the variable/diagnostic x draw x chain shape into a single vector? If we used vectors of vectors then we couldn’t guarantee memory locality for map. We’re probably have to do a copy at least once to put everything together, but let’s figure out the signature that will make it easiest for the interfaces to pack everything together and for the C++ to unpack into Eigen objects.

---

<div class="post-metadata">

### Author: ![ariddell](https://avatars.discourse-cdn.com/v4/letter/a/57b2e6/32.png) [@ariddell](https://discourse.mc-stan.org/u/ariddell)
#### Post date: [March 12, 2018, 12:24pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486/10 "2018-03-12T12:24:07Z")

</div>

Chains are stored in separate holders in Python and R. And the variable/diagnostic \times draw is stored using memory-contiguous column-major-ordered holders (e.g., numpy array in Python) by the interfaces so those can be referenced with a `double*`.

---

<div class="post-metadata">

### Author: ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [March 12, 2018, 3:18pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486/11 "2018-03-12T15:18:57Z")

</div>

> [@betanalpha](#):
>
> how would we get the variable/diagnostic x draw x chain shape into a single vector?

Is there a place where we operate on a `Chain x Draw` object any other way than as a vector for each chain? I don’t see a need for `Variable x Chain x Draw`.

> [@betanalpha](#):
>
> the signature that will make it easiest for the interfaces to pack everything together and for the C++ to unpack into Eigen objects.

If you have `vector<double*>`, it’s easy to map each entry to an Eigen vector or row vector. So if the answer to the above is that all the underlying operations use Eigen vectors, this is the most economical way to pass things that doesn’t implicate Eigen on the interface side, which is what Allen was requesting.

---

<div class="post-metadata">

### Author: ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)
#### Post date: [March 12, 2018, 5:35pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486/12 "2018-03-12T17:35:45Z")

</div>

For the ESS calculations as their are implemented we’d need matrices with shape draw x integration – so can the interfaces guarantee that the info from each chain points to contiguous memory? If so then we’d be looking at a signature like

```
error_code_type diagnostic_function(std::vector<double*>&, 
                                    std::vector<std::string>&, 
                                    int n_var, int n_draw, writer& writer);

```

---

<div class="post-metadata">

### Author: ![syclik](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/syclik/32/6_2.png) [@syclik](https://discourse.mc-stan.org/u/syclik)
#### Post date: [March 12, 2018, 6:14pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486/13 "2018-03-12T18:14:29Z")

</div>

> [@betanalpha](#):
>
> can the interfaces guarantee that the info from each chain points to contiguous memory?

I think we could, or at least request and assume that it does.

> [@betanalpha](#):
>
> If so then we’d be looking at a signature like
> 
> error\_code\_type diagnostic\_function(std::vector\<double\*\>&,  
> std::vector\<std::string\>&,  
> int n\_var, int n\_draw, writer& writer);

What’s the motivation for the `std::string`, `n_var`, and `writer`? Can it not just return a value for most diagnostics?

---

<div class="post-metadata">

### Author: ![syclik](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/syclik/32/6_2.png) [@syclik](https://discourse.mc-stan.org/u/syclik)
#### Post date: [March 12, 2018, 7:44pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486/14 "2018-03-12T19:44:52Z")

</div>

For example, I think ess could have the signature:

```
double ess(std::vector<const double*>& x, int size);

```

And I see why the writers would be useful, so maybe:

```
double ess(std::vector<const double*>& x, int size, writer& error, writer& warning, writer& log);

```

For this to really work, the interfaces need to store each variable within a chain in a contiguous block of memory. I think we could require that.

---

<div class="post-metadata">

### Author: ![ariddell](https://avatars.discourse-cdn.com/v4/letter/a/57b2e6/32.png) [@ariddell](https://discourse.mc-stan.org/u/ariddell)
#### Post date: [March 12, 2018, 9:00pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486/15 "2018-03-12T21:00:28Z")

</div>

I’m not sure we need the writers for ess. They add a lot of complexity.

Could we use exceptions instead? Or return a std::pair with the first  
value being the result and the second value being a string containing  
the error message (or empty if no error). (This is a pattern Go uses.)

---

<div class="post-metadata">

### Author: ![syclik](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/syclik/32/6_2.png) [@syclik](https://discourse.mc-stan.org/u/syclik)
#### Post date: [March 13, 2018, 2:21am UTC](https://discourse.mc-stan.org/t/analysis-api/3486/16 "2018-03-13T02:21:10Z")

</div>

Yes, was already thinking that exceptions could handle exceptional events. I don’t think ess has any output other than the return (we don’t need any writers), but other diagnostic functions might?

---

<div class="post-metadata">

### Author: ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)
#### Post date: [March 13, 2018, 11:12am UTC](https://discourse.mc-stan.org/t/analysis-api/3486/17 "2018-03-13T11:12:14Z")

</div>

> [@syclik](#):
>
> Yes, was already thinking that exceptions could handle exceptional events. I don’t think ess has any output other than the return (we don’t need any writers), but other diagnostic functions might?

Nothing that couldn’t be returned and _then_ passed to a specialized writer method.

---

<div class="post-metadata">

### Author: ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [March 13, 2018, 9:14pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486/18 "2018-03-13T21:14:29Z")

</div>

> [@betanalpha](#):
>
> the info from each chain points to contiguous memory

Yes, the `double*` for each chain would be contiguous, but they wouldn’t be contiguous with each other.

> [@ariddell](#):
>
> Could we use exceptions instead? Or return a std::pair with the first value being the result and the second value being a string containing the error message

I’d much prefer exceptions if that’s workable from the interface perspectives.

If there’s an error, the value won’t make sense. Returning both the value and error message is then clunky as only one will ever get used. This error code return is also the tradition in C with functions with signatures of the form `int foo(..., result_t* result)` where the return is a return code indicating error or not and the function setting `result` with the actual value calculated.

---

<div class="post-metadata">

### Author: ![ariddell](https://avatars.discourse-cdn.com/v4/letter/a/57b2e6/32.png) [@ariddell](https://discourse.mc-stan.org/u/ariddell)
#### Post date: [March 13, 2018, 10:10pm UTC](https://discourse.mc-stan.org/t/analysis-api/3486/19 "2018-03-13T22:10:12Z")

</div>

Exceptions are well-supported in Python (Cython) and R (Rcpp). The Stan  
library already uses them in a few places.

---

<div class="post-metadata">

### Author: ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)
#### Post date: [March 14, 2018, 1:20am UTC](https://discourse.mc-stan.org/t/analysis-api/3486/20 "2018-03-14T01:20:32Z")

</div>

Unfortunately the conversation has diverged a bit here. Let me try to bring it back a bit.

1. We do not throw exceptions to the interfaces for the algorithm API. This thread is not about redoing the API but adding additional routes for common diagnostic calculations. If there is interest in changing the client relationship and how error information is propagated then a new thread should be started.

2. Diagnostic functions return various statuses depending on the thresholds we set (and possibly modifiable by the user). Okay, Warning, and Error messages do not indicate that the function has failed but rather that the input fit information is good or suspect. Exceptions would be more appropriate for failures to calculate the diagnostics.

3. The question introduced in this thread is designing a uniform interface for all diagnostic functions that is sufficiently robust for current (and immediate future diagnostics) that is compatible with the algorithm API.

[Next page](https://discourse.mc-stan.org/t/analysis-api/3486.md?page=2)
