0
$\begingroup$

I have the following strange phenomenon which puzzles me!:

I have a piecewise constant probability density given as

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());

When I draw n random numbers (radius of sphere in millimeters) from this distribution and compute the mass of the sphere and sum them up like:

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

I get mass = 4.3283 kg (see LIVE and the c++ post)

Doing the EXACT identical thing in Mathematica like:

Graphic

Gives the assumably correct value of 4.5287 kg.

Which is not the same, also with different seeds , C++ and Mathematica never match! ? Is that numeric inaccuracy, which I doubt it is...? Question : What the hack is wrong with the sampling in C++?

Simple Mathematica Code:

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}]
$\endgroup$

closed as off-topic by Mr.Wizard Mar 10 '15 at 2:49

  • The question does not concern the technical computing software Mathematica by Wolfram Research. Please see the help center to find out about the topics that can be asked here.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 2
    $\begingroup$ the pdf you supply to ProbabilityDistribution must be a proper probability density function in the sense that it integrates to unity over the interval. ( yours integrates to ~1/2 ) $\endgroup$ – george2079 Mar 4 '15 at 19:35
  • $\begingroup$ ups jeah, sorry, should be 2 times this, (edit) $\endgroup$ – Gabriel Mar 4 '15 at 19:38
  • 1
    $\begingroup$ Ok, Thanks. I see that with n=2000000 your C-code gives 8.65219 and Mathematica 9.0509. (The difference increases with n). I suppose this should not be correct, isn't it? $\endgroup$ – Dargor Mar 4 '15 at 19:58
  • 2
    $\begingroup$ Isn't this more of a C++ question than a Mathematica one? $\endgroup$ – Simon Woods Mar 4 '15 at 21:29
  • 1
    $\begingroup$ The c++ question is here: stackoverflow.com/questions/28862895/… $\endgroup$ – Gabriel Mar 4 '15 at 23:59
6
$\begingroup$

a bit of an extended comment,

Note there is no need for a delayed defintion of your pdf:

 pdf[r_] = 
   Simplify[2 (Piecewise[{{0, r <= 0.59}, {1.36814, 
     Inequality[0.59, Less, r, LessEqual, 0.7]}, {0, r > 0.7}}, 
   Indeterminate] + 
  Piecewise[{{0, r <= 0.7}, {1.99139, 
     Inequality[0.7, Less, r, LessEqual, 0.85]}, {0, r > 0.85}}, 
   Indeterminate] + 
  Piecewise[{{0., r <= 0.85}, {0.29116, 
     Inequality[0.85, Less, r, LessEqual, 1]}, {0., r > 1}}, 
   Indeterminate] + 
  Piecewise[{{0, r <= 1}, {0.039562, 
     Inequality[1, Less, r, LessEqual, 1.18]}, {0, r > 1.18}}, 
   Indeterminate])];

This integrates to nearly 1, but lets normalize it so its exact: (the factor is ~0.999998 )

 pdf[r_] =  pdf[r]/Integrate[ pdf[r], {r, .59, 1.18}] // Simplify

enter image description here

Now you can get your result analytically: (this fails if you don't do the normalization)

 PDFr = ProbabilityDistribution[pdf[r], {r, 0.59, 1.18}];
 1000000 Integrate[ 2400 ( 2 InverseCDF[PDFr, p] 10^-3)^3 Pi/6  , {p, 0, 1}]

4.52594

This is true because: Math

This seems fairly convincing that your mathematica monte carlo is correct and the trouble lies in the c version.

$\endgroup$
  • $\begingroup$ Jeap, I tested that too! Thanks for the comment! it seems that the c++ is buggy! $\endgroup$ – Gabriel Mar 4 '15 at 20:32
  • $\begingroup$ Sorry, I was mistaken, I did not test that, can you quickly give me a hint, why you integrate the inverse cdf? The term Integrate[ 2400 ( 2 InverseCDF[PDFr, p] 10^-3)^3 Pi/6 , {p, 0, 1}]is the expected mass of one sphere I assume, but why :-)? $\endgroup$ – Gabriel Mar 4 '15 at 23:15

Not the answer you're looking for? Browse other questions tagged or ask your own question.