# How to prevent Mathematica from rounding some PDF values to 0?

Dear Mathematica community,
I'm having a problem with the PDF function for the hypoexponential distribution. Specifically, Mathematica returns 0 as soon as the value is very small. Below is an example:

b = .08;
L = 775;
B[q_] := b*(1-q/L)*q;
H0[q1_, q2_, t1_, t2_] := PDF[HypoexponentialDistribution[Table[B[k], {k, q1, q2}]], t2-t1]/B[q2] // N;


Then I get: $$H0[247, 300, 0, 1]=1.47678\times10^{-15}$$ and $$H0[246, 300, 0, 1]=0.0$$. As you can see, the output jumps from 10^(-15) to 0 as the input decreases from 247 to 246.
I read similar topics (e.g., here), but I couldn't make their solution work for my problem.
Any suggestions are more than welcome! Please let me know if you need me to clarify anything. Thanks in advance for your help!

Edit n°1: I did what @Nasser and @user293787 advised me to do (thank you both for your quick reply), but below is what I get:

Edit n°2:
I used the method suggested by Bob Hanlon, which helped me a lot. However, there are still some input values for which the function H0 returns 0. For example, if I have b=13/200 and L=926, I get H0[47, 74, 0, 1]=2.77146*10^-15, but H0[46, 74, 0, 1]=0. I thought it would be enough to increase the calcPrep and/or MaxExtraPrecision, but apparently not. Any suggestions? It would also be great if one of you could explain to me or share with me some webpages/documents that explain how to deal with Mathematica precision.

• does using 8/100 instead of 0.08 fixes your problem? Screen shot !Mathematica graphics Commented Oct 19, 2022 at 14:10
• Thank you @Nasser for your quick reply! I just tried what you advised me to do. Surprisingly, the things get worse... I replaced b=.08 with b=8/100 and I got H0[247, 300, 0, 1]=-1.35129*10^72 and H0[246, 300, 0, 1]=1.87379*10^73... Commented Oct 19, 2022 at 14:38
• Drop the final //N in the definition of H0. Commented Oct 19, 2022 at 14:59
• Thank you very much, @user293787! It works! Commented Oct 19, 2022 at 15:03

The same as suggested by user293787 except most of the work is moved into the function definition.

$Version (* "13.1.0 for Mac OS X x86 (64-bit) (June 16, 2022)" *) Clear["Global*"] b = .08 // Rationalize (* req'd to maintain precision in B *); L = 775; B[q_] := b*(1 - q/L)*q; H0[q1v_, q2v_, t1v_, t2v_, outPrec_ : MachinePrecision, calcPrec_ : 200] := Module[{ q1 = Rationalize[q1v, 0], q2 = Rationalize[q2v, 0], t1 = Rationalize[t1v, 0], t2 = Rationalize[t2v, 0]}, Block[{$MaxExtraPrecision = 200},
PDF[
HypoexponentialDistribution[
Table[B[k], {k, q1, q2}]], t2 - t1]/B[q2] //
N[#, calcPrec] & //
N[#, outPrec] &]];


Examples:

H0[247, 300, 0, 1]

(* 1.47678*10^-15 *)

H0[247, 300, 0, 1, 10]

(* 1.476776375*10^-15 *)

H0[246, 300, 0, 1]

(* 3.72087*10^-16 *)

H0[246, 300, 0, 1, 10]

(* 3.720866046*10^-16 *)


Using 8/100 as pointed out by @Nasser and some extra precision:

b=8/100;
L=775;
B[q_]:=b*(1-q/L)*q;
H0[q1_,q2_,t1_,t2_]:=PDF[HypoexponentialDistribution[Table[B[k],{k,q1,q2}]],t2-t1]/B[q2];

Block[{\$MaxExtraPrecision=1000},
N[N[H0[246,300,0,1],1000]]]
(* 3.72087*10^-16 *)
`