Wrapped Normal Distribution C++

Question

Does anyone know of anyway to generate a wrapped normal distribution in C++? I have a mean and circular standard deviation and was wondering whether C++ had some function which would generate values from a wrapped normal distribution, similar to their standard normal distribution.

Answer 1

In short, there is no native function in C++ to generate a wrapped normal distribution, but there are several easy ways to calculate it.

1. Use Boost

The wrapped normal distribution (WND) is related to the Jacobi theta function (see here). Specifically, it can be formulated in terms of the "third" Jacobi theta function that is available in Boost:

    #include <boost/math/special_functions/jacobi_theta.hpp>

    /// Calculate the density function at angle given sigma and mean my.
    double wnd_density(double angle,double my,double sigma)
    {
        double a=(angle-my)/180.0*M_PI;
        double s=sigma/180.0*M_PI;
        return 1.0/(2.0*M_PI)*boost::math::jacobi_theta3tau(0.5*a,(0.5*s*s)/M_PI);
    }

2. Use iterative formula

Unfortunately Boost does not provide an easy way to implement the cumulative WND. Nevertheless, it is rather easy to implement this (and also the WND) using a simple iteration for the third Jacobi theta function. The details are given in the Wikipedia entry. The same trick can be used to calculate the indefinite integral that can be derived from the iterative formula by integration per substitution, see here.

For larger values of sigma both iteration will converge quickly to machine precision (roughly in under 30 iterations for sigma above 20°), smaller values of sigma take around 800 iterations before reaching machine precision. If this is too much you can explicitly set a bigger eps in the given code.

    /// Calculate the indefinite integral of the third Jacobi Theta function as described in
    /// https://en.wikipedia.org/wiki/Theta_function
    /// Set eps<0 to iterate until machine precision.
    double jacobi_theta3_iterate(double z,double q,double eps=-1.0)
    {
        double th1=1.0,acc=1.0;
        int k=0;
        // This will converge depending on the value of q.
        if(eps<0.0)
        {
            while(th1+acc>th1)
            {
                k++;
                acc=2.0*pow(q,k*k); // an estimate for the accuracy
                th1+=acc*cos(2.0*z*k);
            }
        }
        else
        {
            while(fabs(acc)>=eps)
            {
                k++;
                acc=2.0*pow(q,k*k); // an estimate for the accuracy
                th1+=acc*cos(2.0*z*k);
            }
        }
        return th1;
    }

    /// Calculate the density function at angle given sigma and mean my.
    double wnd_density(double angle,double my,double sigma)
    {
        double a=(angle-my)/180.0*M_PI;
        double s=sigma/180.0*M_PI;
        return 1.0/(2.0*M_PI)*jacobi_theta3_iterate(0.5*a,exp(-0.5*s*s));
    }

    /// Calculate the indefinite integral of the third Jacobi Theta function as described in
    /// https://functions.wolfram.com/EllipticFunctions/EllipticTheta3/21/01/01/
    /// Set eps<0 to iterate until machine precision.
    double jacobi_theta3_integral(double z,double q,double eps=-1.0)
    {
        double th1=z,acc=1.0;
        int k=0;
        // This will converge depending on the value of q.
        if(eps<0.0)
        {
            while(th1+acc>th1)
            {
                k++;
                acc=pow(q,k*k)/k; // an estimate for the accuracy
                th1+=acc*sin(2.0*z*k);
            }
        }
        else
        {
            while(fabs(acc)>eps)
            {
                k++;
                acc=pow(q,k*k)/k; // an estimate for the accuracy
                th1+=acc*sin(2.0*z*k);
            }
        }
        return th1;
    }

    /// Calculate the cumulated density function at angle given sigma and mean my.
    double wnd_cumulated_density(double angle,double my,double sigma)
    {
        double a=(angle-my)/180.0*M_PI;
        double za=(0.0-my)/180.0*M_PI;
        double s=sigma/180.0*M_PI;
        double zero_v=jacobi_theta3_integral(0.5*za,exp(-0.5*s*s));
        return 1.0/M_PI*(jacobi_theta3_integral(0.5*a,exp(-0.5*s*s))-zero_v);
    }

You can save a lot of runtime when you use the loop to also iteratively evaluate the iterations of acc instead of the costly pow function.