包裹正态分布 C++

Question

有谁知道在 C++ 中生成一个包装的正态分布？我有一个平均和循环标准偏差，想知道 C++ 是否有一些函数可以从一个包裹的正态分布中生成值，类似于它们的标准正态分布。

Answer 1

简而言之，C++ 中没有生成包裹正态分布的原生函数，但有几种简单的计算方法。

1。使用 Boost

包裹正态分布 (WND) 与 Jacobi theta 函数 (see here) 相关。具体来说，它可以用 Boost 中可用的“第三个”Jacobi theta 函数来表述：

    #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。使用迭代公式

不幸的是，Boost 没有提供一种简单的方法来实现累积 WND。然而，使用第三个 Jacobi theta 函数的简单迭代来实现这一点（以及 WND）相当容易。详细信息在the Wikipedia entry 中给出。 可以使用相同的技巧来计算不定积分，该不定积分可以通过每次替换的积分从迭代公式导出，see here。

对于较大的 sigma 值，两次迭代都将快速收敛到机器精度（对于 20° 以上的 sigma，大约需要 30 次迭代），较小的 sigma 值需要大约 800 次迭代才能达到机器精度。如果这太多了，您可以在给定代码中显式设置更大的 eps。

    /// 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);
    }

当您使用循环来迭代地评估 acc 的迭代而不是昂贵的 pow 函数时，可以节省大量运行时间。