C++ 随机数生成：生成 cos 平方函数答案

【问题标题】：C++ Random Number Generation: Generate cos squared functionC++ 随机数生成：生成 cos 平方函数
【发布时间】：2017-02-26 13:01:39
【问题描述】：

感兴趣的概率分布是

double x; // range: -pi/2.0 to +pi/2.0
double y = std::pow(std::cos(x), 2.0);

这个函数可以解析积分，但不能倒置。因此，无法执行将均匀分布映射到所需概率分布的通常技巧。

还有其他方法可以用来生成随机变量 cos^2(theta) 分布吗？

也许可以通过数值方法找到反函数，但是我不知道有一种有效的（记忆和计算）方法可以做到这一点。

【问题讨论】：

拒绝抽样是否可以接受？
这不可能是分布 - 它在该范围内的积分不等于 1。如果您按 2/pi 缩放它可以挽救它。为了省事，这个案例的累积分布函数是F(x) = (x + sin(2*x)/2 + pi/2) / pi
@IgorTandetnik 好点——我没有对其进行归一化，因为归一化取决于 x 的范围——因为我写了这篇文章，所以我认为使用 0 到 pi/2.0 的范围可能更容易，而不是问题中指定的原始范围。

标签： c++ c++11 random probability probability-density

【解决方案1】：

来自Inverse transform sampling：你可以从任何概率分布中随机生成样本数，给定它的 cdf。

假设你想要 cos²x 分布，从 -pi/2 到 pi/2。由于 cos²x 从 -pi/2 到 pi/2 的积分是 pi/2，因此您需要按比例缩小以使积分为 1。因此，pdf P(x) = (2 /pi)cos²x

下一步是根据给定的 pdf 计算 cdf，这是 pdf 的积分。您可以使用任何数值方法来求 P(x) 的积分。或者你可以去 Wolfram Alpha 得到答案：cdf is F(x) = (2/pi)(0.5x + 0.25sin2x) + 0.5

接下来您需要计算 F^-1(x)。由于 F(x) 是单调递增函数，因此可以使用二分法（二分查找）轻松找到 F^-1(x)。 Wolfram Alpha 没有这个 F^-1(x) 公式。

然后生成一个从0到1的统一实数u。你的自定义分布是F^-1(u)。

#include <iostream>
#include <cmath>
#include <random>
#include <boost/random/random_device.hpp>
#include <vector>
#include <iomanip>

const double pi = 3.14159265358979323846;

const double LOW = -pi/2;
const double HIGH = pi/2;
double pdf(double x)
{
    return cos(x) * cos(x);
}

double cdf(double x) //integral of pdf
{
    return (2/pi)*(x/2 + sin(2*x)/4) + 0.5; //from Wolfram Alpha
}

double inverse_cdf(double u)
{   //bisection, not 100% accurate
    double low  = LOW;
    double high = HIGH;
    double epsilon = 1e-10; //any small number, e.g. 1e-15
    while (high - low > epsilon)
    {
        double mid = (low + high) / 2;
        if (cdf(mid) == u) return mid;
        if (cdf(mid) < u) low = mid; else high = mid;
    }
    return (low + high) / 2;
}

double custom_distribution(std::mt19937& rng)
{
    double u = std::uniform_real_distribution<double>(0,1)(rng);
    return inverse_cdf(u);
}

int main()
{
    std::mt19937 rng{boost::random::random_device{}()};

    std::vector<double> xCount(15);
    int nSamples = 10000;
    double gap = (HIGH-LOW) / xCount.size();
    while (nSamples--) xCount[(int)( (custom_distribution(rng) - LOW) / gap )]++;
    for (int i = 0; i < xCount.size(); ++i)
    {
        std::cout << std::setw(2) << i << ":" << xCount[i] << "\t";
        for (int bar = xCount[i]/15; bar--; std::cout << '*');
        std::cout << "\n";
    }
}

样本输出：

 0:17   *
 1:135  *********
 2:305  ********************
 3:604  ****************************************
 4:859  *********************************************************
 5:1106 *************************************************************************
 6:1256 ***********************************************************************************
 7:1353 ******************************************************************************************
 8:1271 ************************************************************************************
 9:1102 *************************************************************************
10:876  **********************************************************
11:614  ****************************************
12:334  **********************
13:143  *********
14:25   *

【讨论】：