Fix bug in SineBivariateVonMises sampler (#1628)

* sbvm sampler fix by comparison with R code

* use rsqrt instead of 1/sqrt

numpyro/distributions/directional.py

@@ -477,7 +477,7 @@ def update_fn(curr):
             phi_key, key = random.split(key)
             accept_key, acg_key, phi_key = random.split(phi_key, 3)
 
-            x = jnp.sqrt(1 + 2 * eig / b0) * random.normal(acg_key, shape)
+            x = lax.rsqrt(1 + 2 * eig / b0) * random.normal(acg_key, shape)
             x /= jnp.linalg.norm(
                 x, axis=1, keepdims=True
             )  # Angular Central Gaussian distribution