Can't instantiate HalfT prior

[niwaka-ame]

Hi,

If I run this code:

from GPy.core.parameterization.priors import HalfT
HalfT(1, 3)

it gives the following error:

Traceback (most recent call last):                                                           
  File "<stdin>", line 1, in <module>                                                        
  File "/home/y/anaconda3/lib/python3.12/site-packages/GPy/core/parameterization/priors.py"$
    o = super(Prior, cls).__new__(cls, A, nu)                                                
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                                                
TypeError: object.__new__() takes exactly one argument (the type to instantiate)      

My version of GPy is 1.13.2.

I think this error is related to #503 and there might be a simple fix, but i'm not confident enough to give it a reliable fix (for research purpose). Could anyone from inside/outside the team help? Many thanks.

Edit: added version of GPy.

[MartinBubel]

Hi @niwaka-ame
Thanks for reporting. I'm busy this week but I will take a look and propose a fix soon.

Best
Martin

[niwaka-ame]

Hi @MartinBubel, thanks very much!

I actually tried to give a temporary workaround (following #503) and I sampled from HalfT(1,1).

However, half of the samples are exactly 0 - which is different from samples from Scipy's HalfCauchy. The number of zeros doesn't change much if I sample from HalfT(1,3) or others.

It might be worth checking when you fix the issue - it could be an issue with my workaround or a bug in the algorithm itself.

Code:

from GPy.core.parameterization.priors import HalfT, Prior
import weakref
import GPy
import numpy as np

class MyHalfT(HalfT):
    def __new__(cls, A, nu):  # Singleton:
        if cls._instances:
            cls._instances[:] = [instance for instance in cls._instances if instance()]
            for instance in cls._instances:
                if instance().A == A and instance().nu == nu:
                    return instance()
        newfunc = super(Prior, cls).__new__
        if newfunc is object.__new__:
            o = newfunc(cls)
        else:
            o = newfunc(cls, A, nu)
        # o = super(Prior, cls).__new__(cls, A, nu)
        cls._instances.append(weakref.ref(o))
        return cls._instances[-1]()

xs = np.random.normal(0, 1, size=(10, 1))
ys = np.random.normal(0, 1, size=(10, 1))
k1 = GPy.kern.RBF(input_dim=1, variance=1.0, lengthscale=1.0)
m = GPy.models.GPRegression(xs, ys, kernel=k1)
m.kern.variance.set_prior(MyHalfT(1, 1))  # halfT is problematic
m.kern.lengthscale.set_prior(LogGaussian(0, 2))  # loggaussian is fine
n = int(5e3)
samples_var = np.zeros(n)
samples_len = np.zeros(n)
for j in range(n):
    m.randomize()
    samples_var[j] = m.kern.variance[0]
    samples_len[j] = m.kern.lengthscale[0]
print(np.sum(samples_var == 0.))

comparing with halfcauchy is scipy:

from scipy.stats import halfcauchy
r = halfcauchy.rvs(size=n)
print(np.sum(r == 0.))

Edit: fixed typo in the halfcauchy code

[MartinBubel]

Hi @niwaka-ame
thanks for the update. I just gave it a try myself and it seems like I have to dig a bit deeper into this. I found that there is no test on the HalfT prior and some others.

The fix you are proposing looks good and is also part of some other prior objects. I still have to check the implementation of the prior in general to see why the result using the scipy prior deviates.