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$
3
$\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$
  • $\begingroup$ Thank you also for the explanation. $\endgroup$ – mattiav27 Sep 24 '18 at 14:43
  • $\begingroup$ I still get negative values: could you show how to use controlled precision? $\endgroup$ – mattiav27 Sep 25 '18 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 '18 at 21:56
  • $\begingroup$ thank you @John Doty $\endgroup$ – mattiav27 Sep 26 '18 at 8:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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