# Where is the error in my computation of the wrapped normal distribution density?

Let

• $$\sigma\in(0,1)$$
• $$\phi(x):=\frac1{\sqrt{2\pi\sigma^2}}e^{-\frac{x^2}{2\sigma^2}}\;\;\;\text{for }x\in\mathbb R$$ and $$\psi(x):=\sum_{k\in\mathbb Z}\phi(x+k)\;\;\;\text{for }x\in\mathbb R$$
• $$\beta\in[0,1]$$
• $$d\in\mathbb N$$ and $$\lambda'$$ denote the Lebesgue measure on $$[0,1)^d$$
• $$u'(x',y'):=\beta+(1-\beta)\prod_{i=1}^d\psi(y'_i-x'_i)\;\;\;\text{for }x',y'\in[0,1)^d$$ and $$Q'(x',B'):=\int_{B'}\lambda'({\rm d}y')u'(x',y')\;\;\;\text{for }(x',B')\in[0,1)^d\times\mathcal B([0,1)^d)$$

I'm using $$Q'$$ as the proposal kernel for the Metropolis-Hastings algorithm. Since I've encountered a huge error in my estimates, I presume that something is wrong with my computation of the density $$u'$$.

I'm unsure what the best way to verify my implementation is, since I've never thought about this before. In the following code you find a complete (simplified) implementation (you can run the code online here: https://coliru.stacked-crooked.com/a/435e82ba145c450c). As an example, I've build a uniformly distributed $$x'\in[0,1)^d$$ and then created independent samples $$y'_0,\ldots,y'_{n-1}\sim Q'(x',\;\cdot\;)$$. The only sensible test which came to my mind was $$\frac1n\sum_{i=0}^{n-1}\frac1{u'(x',y'_i)}\xrightarrow{n\to\infty}1\tag1.$$ I've printed the result of this estimate to the command line. Note that the function sampler::density returns $$u'(x',y'_i)$$ for the $$i$$ of the current iteration:

#include <algorithm>
#include <cassert>
#include <iostream>
#include <random>

template<typename T = double>
T const pi = std::acos(-T(1));

template<typename T = double>
auto normal_distribution_density(T x, T mu = 0, T sigma = 1)
{
auto const y = (x - mu) / sigma;
return 1 / (sigma * std::sqrt(2 * pi<T>)) * std::exp(-y * y / 2);
}

template<typename T = double>
auto wrapped_normal_distribution_density(T x, T mu = 0, T sigma = 1, T epsilon = 0)
{
T s = normal_distribution_density(x, mu, sigma);
for (int k = 1;; ++k)
{
T a = normal_distribution_density(x + k, mu, sigma),
b = normal_distribution_density(x - k, mu, sigma);
s += a + b;
if (a <= epsilon && b <= epsilon)
break;
}
return s;
}

template<typename RealType = double>
class sampler
{
public:
using real_type = RealType;

sampler(real_type beta, real_type sigma, std::vector<real_type> const& x)
: m_beta(beta),
m_sigma(sigma),
m_x(x)
{}

template<class Generator>
void begin_iteration(Generator& g)
{
m_large_step = m_uniform_distribution(g) < m_beta;
m_sample_index = 0;
m_second_density = 1;
}

template<class Generator>
real_type generate(Generator& g)
{
assert(m_sample_index < m_x.size());

real_type sample;
if (!m_large_step)
{
std::normal_distribution<real_type> normal_distribution(
m_x[m_sample_index], m_sigma);
auto const normal_sample = normal_distribution(g);
sample = normal_sample - std::floor(normal_sample);
}
else
sample = m_uniform_distribution(g);

m_second_density *= wrapped_normal_distribution_density(
sample, m_x[m_sample_index], m_sigma);

++m_sample_index;
return sample;
}

real_type density() const
{
assert(m_sample_index == m_x.size());
return m_beta + (1 - m_beta) * m_second_density;
}

private:
real_type m_beta,
m_sigma;
std::uniform_real_distribution<real_type> m_uniform_distribution;

std::size_t m_sample_index;
std::vector<real_type> const& m_x;

bool m_large_step;
real_type m_second_density;
};

int main()
{
std::size_t d = 1;
std::mt19937 g{ std::random_device{}() };

std::vector<double> x;
x.reserve(d);
{// initialize x ~ U_{[0, 1)^d}
std::uniform_real_distribution<> u;
std::generate_n(std::back_inserter(x), d, [&]() { return u(g); });
}

double const beta = .3,
sigma = .01;
sampler<> s{ beta, sigma, x };

std::size_t const n = 1e6;
double acc{};
for (std::size_t i = 0; i < n; ++i)
{
s.begin_iteration(g);

std::vector<double> y;
y.reserve(d);
std::generate_n(std::back_inserter(y), d, [&]() { return s.generate(g); });

acc += 1 / s.density();
}

acc /= n;
std::cout << acc << std::endl;

return 0;
}


While the error is not as huge as in my Metropolis-Hastings estimates, the computed result is significantly off from $$1$$.

Maybe there error is due to floating point imprecision. Should I compute something in a different way? And please feel free to tell me if there are other simple tests which I might want to consider.

While using the harmonic mean of the simulation densities as an estimator of one is "the worst Monte Carlo method ever", I checked its convergence by coding the simulation from $$u(x,\cdot)$$ on my own and I did not spot any discrepancy:

> mean(1/propz(simox()))
 0.9945046
> mean(1/propz(simox()))
 1.001786


Here is my R code for completion.

wrap<-function(x, mu, sigma=.1){
termini = trunc(5*sigma + 2)
s = sum(dnorm(x + (-termini):termini, mu, sigma))
return(s)}

simox = function(N=1e4,beta=.5,mu=.5,sigma=.1){
unz = (runif(N)<beta)
termini = trunc(5*sigma + 2)
prbz = pnorm(-mu + (-termini):termini, sd=sigma)
qrbz = diff(prbz)
ndx = sample((-termini+1):termini,N,rep=TRUE,pr=qrbz)+termini
z = sigma*qnorm(prbz[ndx]+runif(N)*qrbz[ndx])-ndx+mu+termini+1
return(c(runif(sum(unz)),z[!unz]))
}

propz<-function(y,beta=.5,mu=.5,sigma=.1){
beta+(1-beta)*apply(as.matrix(y),1,wrap,mu=mu,sigma=sigma)
}

• +1 Thank you for your answer. Well, the discrepancy is that $0.9945046$ and $1.001786$ are significantly off from $1$, but I guess this is not due to a computational error but the worseness of the estimator. Mar 23 '20 at 14:21
• As you may know, I've asked this question in the hope that a bad computation is causing the error I've described in my other question. However, I've noticed that even in the case $\beta=1$ (in which $u'\equiv 1$) my summation of $\frac1n\sum_{i=1}^n\rho'(X'_{i-1},Y'_i)=\frac1n\sum_{i=1}^n\left(\frac pq\circ\varphi\right)(Y'_i)$ is dramatically off from the value of $c$ which I've computed by ordinary Monte Carlo integration. I've no clue what's going wrong. Mar 23 '20 at 14:21
• As I cannot invest time to check alien code, I have no clue either. The only generic reason I can think of would be that $q$ has too thin a tail compared with $p$, which would imply the variance of the IS estimator is infinity. Mar 23 '20 at 14:38

Observe that \begin{align*} \log \psi(x) &= \log \sum_{k\in\mathbb Z}\phi(x+k) \\ &= \log \sum_k \exp\left[ \log \phi(x+k) \right] \\ &= m + \log \sum_k \exp\left[ \log \phi(x+k) - m \right] \end{align*}

The goal will be to perform the sum in log-space, and then exponentiate at the very end to get $$\psi(x)$$.

First, if we take $$m := \max_k\{ \log \phi(x+k) \}$$, then the exponentiations inside the sum will have a better chance of not underflowing to $$0$$.

Second, we can calculate $$\log \phi(x+k)$$ intelligently. We won't exponentiate and normalize, only to take the log again. Instead we will use $$\log \phi(x) = -\log(\sigma) - \frac{1}{2}\log(2\pi) - \frac{ x^2}{2\sigma^2}.$$ NB: this is nonstandard notation..typically $$\phi$$ refers to the normal density with mean $$0$$ and variance $$1$$.

Here's the code:

#include <iostream>
#include <cmath>
#include <limits>

const double log_pi (1.1447298858494);
const double log_two_pi (1.83787706640935);
const double inv_sqrt_2pi(0.3989422804014327);

template<typename T = double>
auto normal_distribution_density(T x, T mu = 0, T sigma = 1)
{
auto const y = (x - mu) / sigma;
return 1 / (sigma * std::sqrt(2 * std::exp(log_pi))) * std::exp(-y * y / 2);
}

// https://g...content-available-to-author-only...b.com/tbrown122387/pf/blob/master/include/rv_eval.h#L140template<typename T = double>
template<typename float_t>
float_t evalUnivNorm(float_t x, float_t mu, float_t sigma, bool log)
{
float_t exponent = -.5*(x - mu)*(x-mu)/(sigma*sigma);
if( sigma > 0.0){
if(log){
return -std::log(sigma) - .5*log_two_pi + exponent;
}else{
return inv_sqrt_2pi * std::exp(exponent) / sigma;
}
}else{
if(log){
return -std::numeric_limits<float_t>::infinity();
}else{
return 0.0;
}
}
}

template<typename T = double>
auto wrapped_normal_distribution_density(T x, T mu = 0, T sigma = 1, T epsilon = 0)
{
T s = normal_distribution_density(x, mu, sigma);
for (int k = 1;; ++k)
{
T a = normal_distribution_density(x + k, mu, sigma),
b = normal_distribution_density(x - k, mu, sigma);
s += a + b;
if (a <= epsilon && b <= epsilon)
break;
}
return s;
}

template<typename float_t>
float_t log_sum_exp(float_t a, float_t b)
{
float_t m = std::max(a,b);
return m + std::log(std::exp(a-m) + std::exp(b-m));
}

template<typename T = double>
auto wrapped_normal_distribution_density_taylor(T x, T mu = 0, T sigma = 1, T epsilon = 0)
{
// k=0
T log_result = evalUnivNorm(x, mu, sigma, true);
T last_iter_log_r;
for (int absk = 1; absk < 1000; absk++)
{
last_iter_log_r = log_result;
log_result = log_sum_exp<T>(log_result, evalUnivNorm<T>(x + absk, mu, sigma, true));
log_result = log_sum_exp<T>(log_result, evalUnivNorm<T>(x - absk, mu, sigma, true));
if (last_iter_log_r == log_result)
return std::exp(log_result);
}
}

int main() {

std::cout << "your eval: " << wrapped_normal_distribution_density<double>(10.0) << "\n";
std::cout << "my eval: " << wrapped_normal_distribution_density_taylor<double>(10.0) << "\n";
return 0;
}


When I run it here, I get the following output:

your eval: 1
mine eval: 1


Note that I'm not using $$m$$ as the max of all the log densities. I am taking pairwise maxima.

Edit: woops I was missing a few iterations. They look equivalent to me.

• Thank you for your answer. Shouldn't it be $\log \phi(x) = -\log(\sigma) - \frac{1}{2}\log(\color{red}2\pi) - \frac{ x^2}{2\sigma^2}.$? Mar 22 '20 at 19:54
• @0xbadf00d yes thanks. Code remains unchanged, though. Mar 22 '20 at 19:58
• (a) I didn't get the thing with the pairwise maxima. Could you explain the idea behind it? Actually, while it's clear to me that your identity $\ln\psi(x)=m+\ln\sum_{k\in\mathbb Z}\exp\left(\ln\phi(x+k)-m\right)$ holds for all $m\in\mathbb R$, I didn't catch the idea behind this either and why we would like to choose $m= \max_{k\in\mathbb Z}\ln\phi(x+k)$. (b) I'm unsure whether you've intended this, but please note that your first assignment to log_result inside the for-loop is never read, since you've immediately overwrite the assigned value. Mar 22 '20 at 20:18
• @0xbadf00d a.) perhaps mc-stan.org/docs/2_18/stan-users-guide/… alternatively, try using both formulas side-by-side on numbers that are either really negative or really large; b.) those lines replace the += operator Mar 22 '20 at 22:04
• (b) Sorry, since I didn't fully understand the theory behind it, I didn't read the code carefully enough. (a) I will need to take a closer look to it. However, I've noticed that for x = 0.22143192861467834, mu = 0.21656738867920472 and sigma = 0.01 the computation of wrapped_normal_distribution_density_taylor(x, mu, sigma) is extremely slow (I've seen over 569389625 iterations and the exit criterion is still not satisfied). Is there any better exit criterion (like the epsilon in my implementation)? Mar 23 '20 at 6:14