Model components don't add up

Thanks for this great piece of software! I have noticed something that I would like to bring up here. When I do this:

model, components = results.rv.evaluate(inst, return_components=True)

the components don't always add up to the model, which I naïvely would expect.

After some digging, I understand what happens: The function selects a random set of samples, then it evaluates the model and the components for each sample, and finally it computes the median of all the samples at each time point. The median is computed separately for the full model and the individual components, meaning that e.g. the output for the "trend" component will not necessarily correspond to the trend component included in the full model. In other words, the median is not a linear operator.

Of couse, this behavior is not necessarily wrong, but it could lead to unexpected results when subtracting e.g. the trend component from the full model (it did for me). My first suggestion is to explain this to the user in the docstring. My second suggestion is to do some sort of "argmedian" on the full model, such that the component values at each time would correspond to the same model, but that probably comes with its own problems.

This is related to issue #40, suggesting an option to return the actual model samples, which would allow the user to handle this the way he/she prefers. That might be the best solution!

nespinoza / juliet

Model components don't add up #56