Build custom random method with reject sampling

Questions
1

I am trying to define sine-skewed von mises distribution(SSVM).
SSVM was not implemented in scipy.stats so I have to write it by my self.

Most of the random methods in pymc3 seem to be using scipy.stats.rvs or inverse transform sampling.
Is it ok to implement custom _random method with reject sampling?
(the cumulative density function of SSVM looks too complicated for inverse transform sampling…)

If it’s ok, how can I include it properly in custom density distribution?

Sine-skewed von mises distribution(SSVM)
Reference : Characterization of Sine- Skewed von Mises Distribution

1.Probability density function
f(\theta)=(1+\lambda\sin(\theta-\mu)) \frac{1}{2\pi I_{0}(\kappa)}\exp(\kappa \cos(\theta - \mu))
where \lambda is skewness(|\lambda| \leq 1), \mu is the mean(|\mu| \leq \pi),
\kappa is the concentration(0 \leq \kappa), I_{0} is the modified Bessel function of order zero.

2.Cumulative density function
F(\theta) = \frac{1}{2\pi} \left[ (\pi + \theta -\mu) + \ 2\sum_{j=1}^{\infty} \frac{I_j(\kappa)}{j} \sin(j(\theta-\mu)) \right] + \ \frac{\lambda}{2\pi I_0(\kappa)}(\exp^{-\kappa} - \exp^{\kappa\cos(\theta-\mu)})
where
I_j(\kappa) = \left( \frac{\kappa}{2} \right)^j \sum_{i=0}^{\infty} \ \left[ \left( \frac{\kappa}{2} \right)^{2i} \left( \frac{1}{i!(j+1)} \right) \right]

Simple rejection sampling works well.

rejection
rejection.png6238×2122 270 KB

2

This is interesting! Adding the custom random method will allow prior predictive sampling, or posterior predictive sampling, whenever the RV is observed.

I see no reason you could not include it - the hardest part will be dealing with shapes. You will need to allow for vectorized inputs, and return similarly vectorized outputs. I would try copying a bunch of the code from the existing Von Mises distribution as a starting point, and implementing a custom generate_samples.

3

Thank you so much for your kind reply!

I have two another questions.
1.Is there any sample of custom generate_samples?
2.Is there any other simple way to implement random of complex distribution?

To be honest, I almost have no statistics nor programming background so I’ll start with studying generate_samples.
Using pymc3 inference is not so hard but understanding the inside is difficult for me now…