Non-Gaussian likelihoods

I'm creating this issue as a place to discuss supporting non-Gaussian likelihoods.

There's been a fair amount of discussion on this for Poisson / Negative Binomial likelihoods (#337). This comes up in problems with small count data. I think a Negative Binomial likelihood would be a very high-value add to the package.

#1492 provides a working example of a Gamma likelihood. As discussed there, the application is for positive and positively-skewed time series.

I see two main ways of getting this functionality in place.

**Handling it like trend**
Selecting different likelihoods is related to the current functionality of selecting different trends. In the current release there are two trends: piecewise linear or piecewise logistic. #1466 just landed and adds a third (flat). In the Stan code we pass in an indicator for which trend to use:
https://github.com/facebook/prophet/blob/16e632a6958bc1fbfcc67ed8628ba8c972df15db/python/stan/unix/prophet.stan#L98
and then compute the trend using the appropriate function:
https://github.com/facebook/prophet/blob/16e632a6958bc1fbfcc67ed8628ba8c972df15db/python/stan/unix/prophet.stan#L118-L124
Then, anywhere in the Py code where the trend is computed, we similarly use a switch to select the correct function. Here is an example, and #1466 is a great example of all of the places where the trend function is used:
https://github.com/facebook/prophet/blob/16e632a6958bc1fbfcc67ed8628ba8c972df15db/python/fbprophet/forecaster.py#L1315-L1323

We could do a similar thing with the likelihood. For Gaussian and Negative Binomial this would be fairly clean. There could be a switch on the likelihood in the Stan code here:
https://github.com/facebook/prophet/blob/16e632a6958bc1fbfcc67ed8628ba8c972df15db/python/stan/unix/prophet.stan#L136
and then in the Py code, the likelihood doesn't play a huge role. In fact the only place in the Py code where the Gaussian assumption shows up is when estimating uncertainty intervals here:
https://github.com/facebook/prophet/blob/16e632a6958bc1fbfcc67ed8628ba8c972df15db/python/fbprophet/forecaster.py#L1466
For NB, that would just have to be replaced with simulating from the NB distribution.

However, there are three main issues with this approach:

- Different likelihoods may have different parameters. Gaussian has the noise variance. Poisson has no parameter. Negative Binomial has the dispersion parameter phi (https://mc-stan.org/docs/2_20/functions-reference/nbalt.html). #1492 shows the differences for Gamma. Trying to override the existing Gaussian noise parameter,
https://github.com/facebook/prophet/blob/16e632a6958bc1fbfcc67ed8628ba8c972df15db/python/stan/unix/prophet.stan#L132
may not be reasonable due to changes in valid ranges / reasonable priors.

- The posterior mean might require distribution-specific computation. With the Gaussian likelihood we have something like `y ~ Normal(trend(t) + seasonality(t) + regressors(t), sigma)` so we can compute `yhat` (the posterior mean) directly as `trend(t) + seasonality(t) + regressors(t)`. The same would be true for the Poisson or Negative Binomial likelihoods. But it isn't true for the Gamma likelihood, where the mean is alpha/beta. I wonder if the Gamma distribution could be reparameterized in such a way to get around this, but that would probably require implementing the reparameterized distribution in Stan and so it'd be nice to have a cleaner solution that doesn't require that.

- Trying out custom link functions would require making changes in several places throughout the code. The actual amount of code would probably be relatively small, but it would take a lot of effort for someone to be sure that they'd gotten all of the right places. #1492 is along the lines of what would be required, as is #1466.

**Creating a Model class**

I think the alternative would be to add a whole bunch more structure to how the model is defined and make it so the available modeling options (primarily trend and likelihood) could be customized just by implementing a new model class (let's call it Model here). The Model class would define the trend and likelihood. All of the interaction in the Prophet class with either the trend or the link function would be handled by calling this class in a way that has a uniform API for all trends/likelihoods. This would handle #705 at the same time as allowing for custom likelihoods. For instance, switch statements like this:
https://github.com/facebook/prophet/blob/16e632a6958bc1fbfcc67ed8628ba8c972df15db/python/fbprophet/forecaster.py#L1315-L1323
would be replaced by something like
```
trend = self.model.evaluate_trend(df, self.params, self.changepoints_t)
```
Code like
https://github.com/facebook/prophet/blob/16e632a6958bc1fbfcc67ed8628ba8c972df15db/python/fbprophet/forecaster.py#L1216-L1219
might be replaced by something like
```
df2['yhat'] = self.model.evaluate_yhat(df2['trend'] * (1 + df2['multiplicative_terms']) + df2['additive_terms'])
```

The main advantage of this is that it would enable trying out custom trends or link functions by implenting code in a single place. But there are some challenges too that I haven't thought through fully enough yet:
- It will be a significant effort to scope out this class and figure out where to draw the line between what stays in the main Prophet class and what is relegated to the Model class and thus made customizable. It seems like computing the linear regressor feature matrix seems would stay in the Prophet class, as would handling future trend uncertainty (which means the piecewise trend formulation would be required for any trend).
- The definitions of model components may depend on the likelihood. The Gamma likelihood in #1492 is a good example of this, where there is now a trend for the alpha parameter and a separate trend for the beta parameter. This will have major implications throughout the stack, including for things like plotting where we assume that there is always a `trend` component that can be plotted.

**What to do next**
I'd be interested to hear people's thoughts on how to make the model likelihood and/or trend easily customizable.

My most immediate interest in changing the likelihood is to add a Negative Binomial likelihood. I think this could be done pretty easily with the switch approach, and the code change required would be even less than #1466. Since this is one of the issues that comes up the most frequently I'm inclined to take the most immediate path on this and just get it done that way. It could also be a nice test case that ensures there aren't any other considerations from swapping out the likelihood that I've missed.

If we want to extend to any additional likelihoods, then I think making some type of an abstraction to produce a clean interface between Prophet and the likelihood would be necessary. The results in #1492 for the Gamma likelihood are quite compelling, and I think makes a perfect example of what needs to be supported by such an abstraction.