下面的帖子已经解决了,这个问题是因为对http://www.cplusplus.com/reference/random/piecewise_constant_distribution/上的公式解释错误导致的强烈建议读者考虑页面:http://en.cppreference.com/w/cpp/numeric/random/piecewise_constant_distribution
我有以下奇怪的现象让我很困惑!:
我有一个分段常数概率密度,如下所示
using RandomGenType = std::mt19937_64;
RandomGenType gen(51651651651);
using PREC = long double;
std::array<PREC,5> intervals {0.59, 0.7, 0.85, 1, 1.18};
std::array<PREC,4> weights {1.36814, 1.99139, 0.29116, 0.039562};
// integral over the pdf to normalize:
PREC normalization =0;
for(unsigned int i=0;i<4;i++){
normalization += weights[i]*(intervals[i+1]-intervals[i]);
}
std::cout << std::setprecision(30) << "Normalization: " << normalization << std::endl;
// normalize all weights (such that the integral gives 1)!
for(auto & w : weights){
w /= normalization;
}
std::piecewise_constant_distribution<PREC>
distribution (intervals.begin(),intervals.end(),weights.begin());
当我从这个分布中抽取 n
个随机数(以毫米为单位的球体半径)并计算球体的质量并将它们相加时:
unsigned int n = 1000000;
double density = 2400;
double mass = 0;
for(int i=0;i<n;i++){
auto d = 2* distribution(gen) * 1e-3;
mass += d*d*d/3.0*M_PI_2*density;
}
我得到mass = 4.3283 kg(见LIVE here)
在 Mathematica 中做完全相同的事情,例如:
给出假定正确的值 4.5287 kg。 (参见 mathematica)
这不一样,而且种子也不同,C++ 和 Mathematica 永远不匹配! ?是数字不准确吗,我怀疑是……? 问题:C++ 中的采样有什么问题?
简单的 Mathematica 代码:
pdf[r_] = 2*Piecewise[{{0, r < 0.59}, {1.36814, 0.59 <= r <= 0.7},
{1.99139, Inequality[0.7, Less, r, LessEqual, 0.85]},
{0.29116, Inequality[0.85, Less, r, LessEqual, 1]},
{0.039562, Inequality[1, Less, r, LessEqual, 1.18]},
{0, r > 1.18}}];
pdfr[r_] = pdf[r] / Integrate[pdf[r], {r, 0, 3}];(*normalize*)
Plot[pdf[r], {r, 0.4, 1.3}, Filling -> Axis]
PDFr = ProbabilityDistribution[pdfr[r], {r, 0, 1.18}];
(*if you put 1.18=2 then we dont get 4.52??*)
SeedRandom[100, Method -> "MersenneTwister"]
dataR = RandomVariate[PDFr, 1000000, WorkingPrecision -> MachinePrecision];
Fold[#1 + (2*#2*10^-3)^3 Pi/6 2400 &, 0, dataR]
(*Analytical Solution*)
PDFr = ProbabilityDistribution[pdfr[r], {r, 0, 3}];
1000000 Integrate[ 2400 (2 InverseCDF[PDFr, p] 10^-3)^3 Pi/6, {p, 0, 1}]
更新: 我做了一些分析:
读入从 Mathematica 生成的数字(64 位 double ) C++ -> 计算总和,结果与 Mathematica 相同
通过还原计算的质量:4.52528010260687096888432279229将 C++(64 位 double )生成的数字读入 Mathematica -> 计算总和,得到相同的 4.32402
我几乎得出结论,
std::piecewise_constant_distribution
的采样不准确(或与 64 位 float 一样准确)或有错误...或者我的有问题权重?在 http://coliru.stacked-crooked.com/a/ca171bf600b5148f 中密度计算错误
std::piecewise_constant_distribution
===> 好像是个bug!
CPP 生成值与所需分布的直方图:
file = NotebookDirectory[] <> "numbersCpp.bin";
dataCPP = BinaryReadList[file, "Real64"];
Hpdf = HistogramDistribution[dataCPP];
h = DiscretePlot[ PDF[ Hpdf, x], {x, 0.4, 1.2, 0.001},
PlotStyle -> Red];
Show[h, p, PlotRange -> All]
文件在这里生成:Number generation CPP
最佳答案
看来 std::piecewise_constant_distribution
的概率公式写错了
http://www.cplusplus.com/reference/random/piecewise_constant_distribution/
权重的求和是在不乘以区间长度的情况下完成的!
正确的公式是: http://en.cppreference.com/w/cpp/numeric/random/piecewise_constant_distribution
这解决了以前发现的每一个愚蠢的怪癖,如错误/浮点错误等等!
关于c++ - 与 Mathematica 相比,C++ 中的 float 学舍入很奇怪,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/28862895/