Baeysian PCA on Poisson data

[pkese]

Hi,

I'm trying to build a PCA model to fit poisson data. My use case is nearly the same as the PCA tutorial, except that my observed data is Poisson distributed.

For example instead of Variable.GaussianFromMeanAndPrecision in https://github.com/dotnet/infer/blob/master/src/Tutorials/BayesianPCA.cs#L102 my case would contain something like
data[observation, feature] = Variable.Poisson(events * T[observation, feature]).

My problem is that Variable.Poisson expects mean rates in Log space, so the above does not compile.
The error code is saying Gaussian is not assignable from Gamma for result of method PoissonOp.MeanAverageLogarithm (interesting thing regarding this error is that my code doesn't contain any gamma distribution - they are all gaussian - yet error report complains about gamma).

On the other hand, if I transform my principal components into log space and do a
data[observation, feature] = Variable.Poisson(events * Variable.Exp( T[observation, feature] )) then the code compiles, but I can't get correct results because MatrixMultiply doesn't do the right thing -- in log-space, weights and factors should be added rather than multiplied.

I'd be very thankful for any hints on how to approach this.

[tminka]

The T variable generated by this model can be negative. Therefore you must have some kind of link function to turn it into a Poisson rate. See section 3.2.2 of "Likelihood-based approaches to modeling the neural code. Bayesian brain: Probabilistic approaches to neural coding".

[pkese]

Understood.

I would nevertheless like to avoid non-linearities...

With

T = Variable.MatrixMultiply(Z, W)
observed_data = Variable.Poisson(events * Variable.Exp( T ))

Poisson gets positive values but then Z and W get skewed. I'd need something like

T = Variable.MatrixMultiply(Variable.Exp(Z), Variable.Exp(W))
observed_data = Variable.Poisson(events * T)

...so that I could reconstruct correct Z and W.
But it doesn't appear that MatrixMultiply would accept anything but Gaussian.

I've also tried Gamma with Poisson and Beta with Binomial, but again MatrixMultiply isn't happy with either Beta or Gamma.

Is there a workaround for that? There's been some hints to NegativeBinomial, but I'm not sure how to approach that...

[tminka]

If you expand the MatrixMultiply as explicit pairwise sums and products and use VariationalMessagePassing, then it will allow you to use Gamma-distributed Z and W.