I can send you also a PDF with a plot of these values, if you need. As you can
see, the left tail is computed wrongly with lambda in interval [25 to 40].
I suspect that the problem lies in the following lines (132-157) of foxglynn.c:
if( lambda >= lambda_25 )
{
/*The starting point for looking for the left truncation point*/
const double start_k = 1.0 / ( sqrt( 2.0 * lambda ) );
/*Here we choose the max possible value of k such that
(m - k * sqrt_lambda - 3/2) >= 0*/
const double stop_k = ( m - 3.0/2.0 ) / sqrt_lambda;
/*Start looking for the left truncation point*/
double k, k_ltp = 0;
/*For lambda >= 25 we can use the upper bound for b_lambda,
here lambda is always at least 400:
b_lambda_sup = 1.05;
b_lambda = ( 1.0 + 1.0 / lambda ) * exp( 1.0 / ( 8.0*lambda ) ) <=
a_lambda_sup*/
const double b_lambda_sup = 1.05;
double k_prime, c_m_inf, result;
for( k = start_k; k <= stop_k; k = k + 1 )
if( b_lambda_sup * exp(- 1.0 / 2.0 * k*k) / ( k * sqrt_2_pi ) <
epsilon / 2.0 )
{
k_ltp = k;
break;
}
/*Finally the left truncation point is found*/
pFG->left = (int) floor( m - k_ltp * sqrt_lambda - 3.0 / 2.0 );
I dont know why stop_k is defined as:
const double stop_k = ( m - 3.0/2.0 ) / sqrt_lambda;
because I couldn't find it in the paper, however I believe that this
upper limit estimate is too small. Maybe a value of lambda would be
more appropriate, but I'm not sure. I think that this code should be
modified in this way:
if( lambda >= lambda_25 )
{
/*The starting point for looking for the left truncation point*/
const double start_k = 1.0 / ( sqrt( 2.0 * lambda ) );
/*Here we choose the max possible value of k such that
(m - k * sqrt_lambda - 3/2) >= 0*/
const double stop_k = lambda; //( m - 3.0/2.0 ) / sqrt_lambda;
<<<<<<<<<<<MODIFIED
/*Start looking for the left truncation point*/
double k, k_ltp = 0;
/*For lambda >= 25 we can use the upper bound for b_lambda,
here lambda is always at least 400:
b_lambda_sup = 1.05;
b_lambda = ( 1.0 + 1.0 / lambda ) * exp( 1.0 / ( 8.0*lambda ) ) <=
a_lambda_sup*/
const double b_lambda_sup = 1.05;
double k_prime, c_m_inf, result;
for( k = start_k; k <= stop_k; k = k + 1 )
if( b_lambda_sup * exp(- 1.0 / 2.0 * k*k) / ( k * sqrt_2_pi ) <
epsilon / 2.0 )
{
k_ltp = k;
break;
}
/*Finally the left truncation point is found*/
pFG->left = (int) floor( m - k_ltp * sqrt_lambda - 3.0 / 2.0 );
pFG->left = max(0, pFG->left); <<<<<<<<<<MODIFIED
With these changes the code seems to behave correctly, but I'm not
completely sure of the upper limit stop_k. Note also the last line,
that checks that pFG->left is never < 0.
P.S.: I don't know how to reply in the mailing list, so send me
e-mails to my email address elvio.amparore AT
gmail.com if you have
any questions.