Calling convention for _GaussianBase.log_likelihood(data)? Missing transpose?

[slinderman]

Hey Matt,
I'm a bit confused about the calling convention for _GaussianBase.log_likelihood(). My understanding is that the multivariate Gaussian should have a Dx1 mean vector mu and a DxD covariance matrix sigma. I am calling log_likelihood with an NxD matrix where each row is a vector I'd like to evaluate the likelihood of, and I expected an Nx1 vector of likelihoods to be returned. However, if I make the call with N=1 then I get back a Dx1 vector where each entry is the same, and if I do this with N=2 it crashes due to invalid matrix sizes.

Looking into the code, it seems that the problem arises from a broadcast issue:
x = np.nan_to_num(x).reshape((-1,D)) - mu (distributions.py:361)
broadcasts to a DxD matrix. I think it should be
x = np.nan_to_num(x).reshape((-1,D)) - mu.T instead, since mu is a Dx1 vector.

Making this change fixes the problem in my case, but I'm afraid I might be using the interface incorrectly. Please let me know.
Thanks!
Scott

[mattjj]

Actually I keep mean vectors to have shape tuple (D,). That is, I don't keep a singleton dimension around to make it a column with shape (D,1).

So transposing wouldn't do anything to my mean vectors. It might be possible to have both conventions work by adding that change; however, it might break things elsewhere.

Would it work for you to follow the mean-vectors-are-vectors-not-matrices convention? (I actually think having proper 1D vectors is a nice advantage numpy has over Matlab.)

[mattjj]

Just to be concrete about the mean vectors:

In [1]: from pybasicbayes import distributions as dst

In [2]: a = dst.Gaussian().max_likelihood(randn(100,2))

In [3]: a.mu
Out[3]: array([ 0.02605143, -0.08394119])

In [4]: a.sigma
Out[4]:
array([[ 0.93461397,  0.06652411],
       [ 0.06652411,  1.06895761]])

[slinderman]

Thanks for clarifying Matt! My misconception about the shape of mu actually went back a few layers in my code. I managed to resample a GaussianFixedCov and get a 2D mean. This occurred when I called resample with a data vector of shape (N,). When I call it with a data vector of shape (N,1) it works fine.

In [3]: hyp = {'mu_0' : np.array([0]), 'lmbda_0' : np.eye(1), 'sigma' : np.eye(1)}
In [4]: g = GaussianFixedCov(**hyp)
In [5]: data = np.random.randn(100)

In [6]: g.resample(data)
In [7]: g.mu
Out[7]: array([[-0.19293857]])

In [8]: g.resample(data[:,None])
In [9]: g.mu
Out[9]: array([ 0.04490252])

This is only a problem in the 1D case, but it might be worth checking? I've fixed my code to use 1D mean vectors, and NxD data arrays and everything works fine again!

[mattjj]

Yeah the 1D case can be a problem. I have separate classes for those cases (the ScalarGaussian ones) but having separate 1D classes seems like an unfortunate state of affairs in some respects. For now, that's just the state of affairs!