c ++ 11 — C ++ Генерация случайных чисел: генерировать функцию cos квадрат
Распределение вероятностей интереса
double x; // range: -pi/2.0 to +pi/2.0
double y = std::pow(std::cos(x), 2.0);
Эта функция может быть интегрирована аналитически, но не может быть инвертирована. Поэтому обычная уловка сопоставления равномерного распределения с требуемым распределением вероятности не может быть выполнена.
Есть ли другой метод, который можно использовать для генерации случайной величины cos ^ 2 (тета) распределения?
Может быть возможно найти обратную функцию численно, однако я не знаю эффективного (запоминающего и вычислительного) способа сделать это.
Решение
От Выборка обратного преобразования: вы можете генерировать числа выборок в произвольном порядке из любого распределения вероятностей, учитывая его cdf.
Скажи, что ты хочешь, потому что2x распределение, от -pi / 2 до pi / 2. Поскольку интеграл от cos2x от -pi / 2 до pi / 2 равно pi / 2, вам нужно уменьшить масштаб так, чтобы интеграл равнялся 1. Таким образом, pdf P (x) = (2 / pi) cos2Икс
Следующим шагом является вычисление cdf из заданного pdf, являющегося интегралом pdf. Вы можете использовать любой численный метод, чтобы найти интеграл от P (x). Или вы можете перейти к Wolfram Alpha и получить ответ: cdf это F (x) = (2 / pi) (0.5x + 0.25sin2x) + 0.5
Далее вам нужно рассчитать F-1(Икс). Поскольку F (x) является монотонно возрастающей функцией, вы можете использовать метод деления пополам (бинарный поиск), чтобы найти F-1(х) легко. Вольфрам Альфа не имеет этого F-1(х) формула хотя.
Затем сгенерируйте равномерное действительное число u от 0 до 1. Ваше пользовательское распределение F-1(И).
#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 *
Другие решения
Других решений пока нет …