0
$\begingroup$

I am using the following code to generate some random number from a custom distribution:

me =511000;
M = 938*10^6;
z = 1;
\[Rho] = 19.32;
\[Tau][T_?NumericQ] := T/M ;
\[Beta][T_?NumericQ] := Sqrt[1 - (1/(\[Tau][T] + 1))^2];
wm[T_?NumericQ] := (2 me \[Beta][T]^2)/(1 - \[Beta][T]^2);
Zs[T_?NumericQ] = z (1 - Exp[-((125 \[Beta][T])/z^(2/3))]);
(*Numero delta prodotti da protone da 100 MeV*)

pr[en_] := 
 ProbabilityDistribution[(0.307075*79/196.96655 (10^6 \[Rho])/2 (
     10^-4 Zs[en]^2)/(\[Beta][en]^2 (w)^2) (1 - (\[Beta][en]^2 (w))/
       wm[en] + (Pi \[Beta][en] Zs[en]^2)/
        137 Sqrt[(w)/wm[en]] (1 - (w)/wm[en])))/
   Integrate[
    0.307075*79/196.96655 (10^6 \[Rho])/2 (
     10^-4 Zs[en]^2)/(\[Beta][en]^2 (w)^2) (1 - (\[Beta][en]^2 (w))/
       wm[en] + (Pi \[Beta][en] Zs[en]^2)/
        137 Sqrt[(w)/wm[en]] (1 - (w)/wm[en])), {w, 10, wm[en]}], {w, 
   10, wm[en]}]

Table[RandomVariate@pr[10^11], 10]

but I get some negative values as a result:

{23.6665, -4.36639*10^9, 18.1011, 12.4939, 10.6993, 
 2.23713*10^8, 12.4214, 12.8859, 24.9426, 67.1775}

Why does it happen? If I understand correctly the documentation {w, 10, wm[en]} sets the possible range of results between 10 and wm[en] and both are positive numbers so I should get no negative numbers.

I use version 11.3 for windows

$\endgroup$

1 Answer 1

4
$\begingroup$

Your definition for your probability distribution contains very small and very large machine numbers. Results of arithmetic on such quantities are approximate, with untracked precision. The wide range of magnitudes apparently causes substantial precision loss.

Converting all machine numbers to exact rationals fixes the problem:

p = pr[10^11] // Rationalize[#, 0] &
Table[RandomVariate[p], 100]

yielding:

{12.0533, 29.3423, 11.0546, 29.1306, 21.2312, 16.6761, 27.9732, 
17.4662, 33.744, 48719.1, 28.8642, 142.309, 12.0813, 953690., 
18.6152, 10.8712, 453819., 14.9188, 31.6657, 16.0243, 41.7741, 
151.724, 14.4706, 13.3611, 285.572, 32.0111, 54.6824, 52.284, 
65.0357, 25.6183, 77.4526, 435.646, 23.935, 10.0398, 18.4747, 
1.07676*10^6, 18.9915, 13.3078, 
3.04228*10^6, 21.6419, 11.6097, 269.668, 16.2735, 19.3627, 12.7885, 
19.2447, 11.5729, 41.7535, 11.8968, 57.5484, 22.0166, 15.2088, 
46.0215, 12.0185, 44.7919, 354108., 206349., 23.9537, 222125., 
10.5913, 13.9834, 279.787, 48.0365, 69.8381, 59.0337, 16.789, 
76.9336, 11.8326, 10.9745, 25.226, 65.1776, 34.8532, 23.2601, 
60.5109, 11.2216, 15.0242, 63.258, 15.5377, 14.8823, 25.6393, 
17.6715, 12.3322, 10.5237, 16.0636, 109.666, 
5.05526*10^6, 11.7306, 14.2062, 15.7391, 431.196, 10.0336, 122058., 
23.2851, 42.7913, 374302., 17.1834, 11.9845, 956747., 59.316, 16.5455}

Suitably chosen controlled precision will presumably also work.

$\endgroup$
4
  • $\begingroup$ Thank you also for the explanation. $\endgroup$
    – mattiav27
    Sep 24, 2018 at 14:43
  • $\begingroup$ I still get negative values: could you show how to use controlled precision? $\endgroup$
    – mattiav27
    Sep 25, 2018 at 9:01
  • $\begingroup$ Using RandomVariate[p, WorkingPrecision -> 30] seems to get rid of the negative numbers as well as the extreme positive outliers. $\endgroup$
    – John Doty
    Sep 25, 2018 at 21:56
  • $\begingroup$ thank you @John Doty $\endgroup$
    – mattiav27
    Sep 26, 2018 at 8:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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