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.
 A: 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;
}


// https://g...content-available-to-author-only...b.com/tbrown122387/pf/blob/23d0e94ff8bd987997efafad66c120c07aaeddee/include/rv_eval.h#L120
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.
A: 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()))
[1] 0.9945046
> mean(1/propz(simox()))
[1] 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)
}

