Negative correlation SineBivariateVonMises

I am using the SineBivariateVonMises distribution and would like to represent a negative correlation between my variables. When sampling from the distribution I can do so with a negative correlation parameter, when doing inference (MCMC with NUTS) I get an error, that no valid initial parameters can be found (hard setting the correlation to -1).
Is this an intended behaviour?
In the documentation this paper is linked:

Probabilistic model for two dependent circular variables Singh, H., Hnizdo, V., and Demchuck, E. (2002)

Looking into it, it defines the correlation parameter \lambda on: -\inf < \lambda < \inf

Thanks in advance!

Hi @mjm, could you try to see log probability of that variable is NaN or inf? You can get a trace which contains value and distribution at each site with

trace = numpyro.handlers.trace(numpyro.handlers.seed(model, rng_seed=0)).get_trace(*args, **kwargs)

where args, kwargs are arguments/keyword arguments of your model.

Hi @fehiepsi,

thanks again for you help. I have my problems with finding the log probability in the trace. But I recreated the error in a small example and got get_trace() to run. Here is a rundown:

Sampling from SineBivariateVonMises

data = SineBivariateVonMises(0, 0, 1, 1, -1).sample(PRNGKey(42), (1000, ))


    config={"phi_loc": CircularReparam(), 
            "psi_loc": CircularReparam(),
def min_example(data_2d):
    phi_loc = sample('phi_loc', VonMises(0, 10))
    psi_loc = sample('psi_loc', VonMises(0, 10))
    phi_conc = sample('phi_conc', Beta(1, 1))
    psi_conc = sample('psi_conc', Beta(1, 1))
    conc = -sample('conc', HalfNormal(1))
    depInd = SineBivariateVonMises(phi_loc, psi_loc, 70 * phi_conc, 70 * psi_conc, conc)
    obs = sample('obs', depInd, obs=data_2d)

Run model

rng_key = PRNGKey(0)

num_warmup, num_samples = 100, 200

# Run NUTS.
kernel = NUTS(min_example)
mcmc = MCMC(
), data)
posterior_samples = mcmc.get_samples()

leads to:


RuntimeError: Cannot find valid initial parameters. Please check your model again.

I can also provide the full call stack.


trace = numpyro.handlers.trace(numpyro.handlers.seed(min_example, rng_seed=0)).get_trace(data)


File .../env/lib/python3.9/site-packages/numpyro/distributions/, in Distribution.sample(self, key, sample_shape)
    238 def sample(self, key, sample_shape=()):
    239     """
    240     Returns a sample from the distribution having shape given by
    241     `sample_shape + batch_shape + event_shape`. Note that when `sample_shape` is non-empty,
    248     :rtype: numpy.ndarray
    249     """
--> 250     raise NotImplementedError


This error seems to be due to the Reparam statements of the model. If I remove them, I get trace to work, but the model gives a warning for bad performance when using continous variables for circular parameters.

Trace output with removed circular reparam

              {'type': 'sample',
               'name': 'phi_loc',
               'fn': <numpyro.distributions.directional.VonMises at 0x7fe43c24ab20>,
               'args': (),
               'kwargs': {'rng_key': DeviceArray([2718843009, 1272950319], dtype=uint32),
                'sample_shape': ()},
               'value': DeviceArray(0.29627275, dtype=float32),
               'scale': None,
               'is_observed': False,
               'intermediates': [],
               'cond_indep_stack': [],
               'infer': {}}),
              {'type': 'sample',
               'name': 'psi_loc',
               'fn': <numpyro.distributions.directional.VonMises at 0x7fe3f87c6880>,
               'args': (),
               'kwargs': {'rng_key': DeviceArray([1278412471, 2182328957], dtype=uint32),
                'sample_shape': ()},
               'value': DeviceArray(-0.03310037, dtype=float32),
               'scale': None,
               'is_observed': False,
               'intermediates': [],
               'cond_indep_stack': [],
               'infer': {}}),
              {'type': 'sample',
               'name': 'phi_conc',
               'fn': <numpyro.distributions.continuous.Beta at 0x7fe4204706a0>,
               'args': (),
               'kwargs': {'rng_key': DeviceArray([4104543539, 3483300570], dtype=uint32),
                'sample_shape': ()},
               'value': DeviceArray(0.2734751, dtype=float32),
               'scale': None,
               'is_observed': False,
               'intermediates': [],
               'cond_indep_stack': [],
               'infer': {}}),
              {'type': 'sample',
               'name': 'psi_conc',
               'fn': <numpyro.distributions.continuous.Beta at 0x7fe3f8457220>,
               'args': (),
               'kwargs': {'rng_key': DeviceArray([1194623263, 2038155241], dtype=uint32),
                'sample_shape': ()},
               'value': DeviceArray(0.6324667, dtype=float32),
               'scale': None,
               'is_observed': False,
               'intermediates': [],
               'cond_indep_stack': [],
               'infer': {}}),
              {'type': 'sample',
               'name': 'conc',
               'fn': <numpyro.distributions.continuous.HalfNormal at 0x7fe4205215b0>,
               'args': (),
               'kwargs': {'rng_key': DeviceArray([2205739499, 3850766070], dtype=uint32),
                'sample_shape': ()},
               'value': DeviceArray(0.6486334, dtype=float32),
               'scale': None,
               'is_observed': False,
               'intermediates': [],
               'cond_indep_stack': [],
               'infer': {}}),
              {'type': 'sample',
               'name': 'obs',
               'fn': <numpyro.distributions.directional.SineBivariateVonMises at 0x7fe3f824cd60>,
               'args': (),
               'kwargs': {'rng_key': None, 'sample_shape': ()},
               'value': DeviceArray([[-0.7632046 , -2.5725281 ],
                            [-1.4259993 ,  0.9532311 ],
                            [ 0.84518456,  0.53240395],
                            [ 0.89588714, -0.06103492],
                            [ 0.50887036, -0.34725928],
                            [-1.3985238 ,  0.53051925]], dtype=float32),
               'scale': None,
               'is_observed': True,
               'intermediates': [],
               'cond_indep_stack': [],
               'infer': {}})])

Now I am a bit confused as to where I find the probability. Sorry I never used trace so far.
Thank you!

Sorry, I didn’t notice that CircularReparam uses ImproperUniform under the hood. To draw a trace with ImproperUniform, we need to use an initialization strategy like

trace = numpyro.handlers.trace(
        numpyro.handlers.seed(min_example, rng_seed=0),

I run your code and saw that we get NaN at the obs site

{name: site["fn"].log_prob(site["value"]) for name, site in trace.items()
 if site["type"] == "sample"}

Then printed out the parameters of the likelihood gives me valid values

d = trace["obs"]["fn"]
{name: getattr(d, name) for name in d.arg_constraints}
{'phi_loc': DeviceArray(0.29627275, dtype=float32),
 'psi_loc': DeviceArray(-0.03310037, dtype=float32),
 'phi_concentration': DeviceArray(19.143257, dtype=float32),
 'psi_concentration': DeviceArray(44.272667, dtype=float32),
 'correlation': DeviceArray(-0.6486334, dtype=float32)}

So I think this is a bug in SineBivariateVonMises. Could you please open an issue in github for this? For the reproducible code, we only need the data and the distribution parameters. Thanks!