1

In a computer program (written in C++), given $x\in[0,1)$ and $\sigma>0$, I need to sample $y$ from the wrapped normal distribution $\mathcal W_{x,\:\sigma^2}$ with mean $x$ and variance $\sigma^2$ and evaluate the density $q$ of $\mathcal W_{x,\:\sigma^2}$ at $y$ at the same time.

If $\varphi_{\sigma^2}$ denotes the density of the normal distribution with mean $0$ and variance $\sigma^2$ and $$\psi_{\sigma^2}(y):=\sum_{k\in\mathbb Z}\varphi_{\sigma^2}(k+y)\;\;\;\text{for }y\in\mathbb R,$$ then $$q(y):=\psi_{\sigma^2}(y-x)\;\;\;\text{for }y\in[0,1)$$ is the density of $\mathcal W_{x,\:\sigma^2}$ with respect to the Lebesgue measure.

We can sample $y$ by drawing $\xi\sim\mathcal N_{0,\:1}$, setting $z:=x+\sigma\xi$ and $y:=z-\lfloor z\rfloor$. But now I would need to evaluate $q(y)$ separately.

So, my question is: Can we somehow sample $y$ in a clever way such that we obtain $q(y)$ as a byproduct? If not, how should we evaluate $q(y)$ (the problematic thing being that it's an infinite sum)?

0xbadf00d
  • 271
  • It's an infinite sum, but for $\sigma$ not much larger than $1$, it will converge pretty quickly. Consider $\sigma = 2$ and $y = 1/2$; the 10th term is $\sim 2.5\times 10^{-6}$, and subsequent terms decrease by more than a factor of $10$. With $\sigma = 1$, the sixth term is $\sim 10^{-7}$ and the seventh $\sim 3\times 10^{-10}$, for example. – jbowman Jan 08 '20 at 16:56
  • @jbowman So, should I sum from $k=0,1,2,\ldots$ adding up $a_k:=\varphi_{\sigma^2}(y+k)$ and $b_k:=\varphi_{\sigma^2}(y-k)$ and stop once $a_k<\epsilon$ for some suitable $\epsilon$ (maybe $\epsilon=10^{-6}$?)? Or should I check $b_k$ as well? – 0xbadf00d Jan 08 '20 at 20:00
  • When $\sigma$ is only a little larger than $1,$ for all practical computational purposes this density is uniform: see the last method at https://stats.stackexchange.com/a/117711/919. For smaller $\sigma,$ as @jbowman indicates, the sum converges very rapidly. It is a theta function; some software includes procedures to evaluate it (buried, for instance, in the KS test). – whuber Jan 08 '20 at 20:13
  • @whuber Thank you for your comment. The $\sigma$ I've got in mind is smaller than $1$, e.g. $\sigma=0.01$. – 0xbadf00d Jan 08 '20 at 20:18
  • 1
    Then you can evaluate the density anywhere with a single term of the series! (At many points--more than about $0.05$ from the mean--you don't need any terms at all, because the density is practically zero there.) – whuber Jan 08 '20 at 20:20
  • @whuber Well, I guess you're right. Didn't thought about that (+1)! However, I guess I can't obtain the density simultaneously with sampling or is there a possibility? – 0xbadf00d Jan 08 '20 at 20:28
  • I don't see why sampling wouldn't work as it always does. Just estimate the density using a circular method so that values near the endpoints ${0,1}$ are treated appropriately. Indeed, you may safely shift all values (modulo $1$) to make any single observation the new origin and treat the problem as estimating a density from a Normal sample. The chance of erring is far less than $10^{-100}.$ But why is there any need? You already know what distribution you are sampling from! – whuber Jan 08 '20 at 20:33
  • @whuber I need the density for an importance sampling correction. I don't understand your suggestion. Do you say that I should sample $y$ as usual and then compute density as described in my first comment (where it is most probably sufficient to only take the first sumand $k=0$ for my choice of $\sigma$)? Or what do you mean by a "circular method"? How would that be different from what I've described in the question? – 0xbadf00d Jan 09 '20 at 05:28
  • In R, just apply this function to a vector of the data: function(x) {y <- (x-x[1]+1/2) %% 1; c((mean(y)+x[1]-1/2) %% 1, sd(y))}) Assuming $\sigma \ll 1$ ($\sigma \lt 1/10$ should work well), it will provide an accurate estimate of the mean and sd. – whuber Jan 09 '20 at 14:06
  • @whuber Thanks for your input, but I need it in C++. – 0xbadf00d Jan 09 '20 at 15:11
  • I'm sure you can write the C++ code that implements this very simple calculation! If not, SE is the place to ask about it. – whuber Jan 09 '20 at 15:15
  • @whuber Of course, I could, but since I never used R I'm not able to fully understand your code ;) – 0xbadf00d Jan 09 '20 at 15:17
  • About the only explanation it needs for someone familiar with C syntax is that %% reduces a value modulo $1.$ – whuber Jan 09 '20 at 15:18
  • How about https://www.researchgate.net/publication/268443292_Efficient_Evaluation_of_the_Probability_Density_Function_of_a_Wrapped_Normal_Distribution ? – Severin Pappadeux Jan 09 '20 at 15:41

0 Answers0