What is required to implement df < 2? In particular many real-world time series exhibit df around 1.7-1.8.

@fritzo That condition is only used for moment-matching, where we need to match covariance of two student-ts. For example, given a, b, we want to find c such that St(b, 0, c * Id) ~ St(a, 0, Id). I believe that we can relax that condition as long as we have a better approximation.

In https://arxiv.org/abs/1703.02428 section 4.3 and https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=8822764 , the authors use KL for such approximation, hence df > 2 is not required. However, using KL requires numerical methods to minimizing KL(St(b, 0, c * Id), St(a, 0, Id)), which is non-trivial. :(

I think that we can leverage tail behaviors of St(b, 0, c * Id) and St(a, 0, Id) to find c such that both distributions have the same 95% confidence interval (or maybe have the same icdf(0.95))... instead of using moment matching or kl as above. But I haven't come up with a formula yet.

When df<2 the df-power moments still exist, so we could try to match those at least for df >= 1

mu = E[x]
sigma = (E|x-mu|^df)^(1/df)

Good idea!! I'll try to find the formula and match one of the moments < df.

This new paper https://arxiv.org/abs/1912.01607 computes general moments for Student's t-distribution. There are many choices for moment order, especially for multidimensional cases. I'll try to pick up a simple one from them.

@fritzo Unfortunately, the implementation is only useful for sequential filtering. :( The reason is the .condition op does not asscociate with .add, that is x + xy.condition(obs) != (x.event_pad(right=...) + xy).condition(y) (on the left hand side, condition will increase dof then add will decrease it to preserve tail; on the right hand side, we will always have dof increased). A StudentTHMM can be implemented using StudentT op but parallel scan will be invalid. Do you think that this op is useful if we can't do parallel scan?

Do you think that this op is useful if we can't do parallel scan?

Well the main purpose of Gaussian and GammaGaussian was to implement parallel-scan algorithms, so maybe StudentT ops are not worth implementing.

What do you think about an alternative approach to StudentTHMM inference:
Following StableHMM, we could implement a StudentTHMM with .rsample() but no .log_prob() method, and implement auxiliary variable reparameterizers StudentTReparam and StudentTHMMReparam using scale-mixtures of Normal and GaussiannHMM, respectively. During inference we would learn/sample (completely in parallel) the auxiliary Gamma-distributed "outlyingness" variables for each time step and observation step; then conditioned on "outlyingness" we could parallel-scan evaluate GaussianHMM.log_prob() or filter. Inference would then be exact using SVI, HMC, NUTS, etc. and would work for both latent and observed HMM observations. I think we would want to implement the MultivariateStudentT versions separately.

re not worth implementing: Agreed! I'll move some cleanups into a separate PR and leave the implementation here in case it is useful for benchmarking in the future.

re alternative approach to StudentTHMM inference: It is reasonable and should be easier to have than StableHMM. I think that it is also more useful than GammaGaussian, whose posteriors of the last state is more Gaussian-like for long time series with high dimensional outputs. I'll follow this approach instead.

I'll follow this [StudentTHMM] approach instead.

Great! I think the code should be straight forward. I see the biggest open research question being:
"What kind of amortized guides are best for estimating the latent scale parameters"? I think we can share guide patterns across StableHMM and StudentTHMM.

I get your point. I'll think more about that question...

The reference for this pull request is interesting. So is there a way to fit a t smoother with Pyro? I am coming from the PyFlux library where, for example, I can fit a dynamic linear model with student’s t noise. PyFlux is no longer being actively developed. I would like to switch to Pyro to fit these types of models if possible.

Hi @matthewfieger , one thing we are working on is a StudentTHMM together with an auxiliary variable reparameterizer StudentTReparam similar to StableHMMReparam. This would allow StudentT processes to be represented as conditional Gaussian processes, permitting smoothing via SVI or HMC. This approach would be more accurate and more expensive than the moment-matching approach proposed in this PR.

EDIT here is the alternative line of work: #2254