# Negative values with RandomVariate

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

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.

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