Dear MRMC developers,
I think I have found another (more severe) problem in the MRMC FoxGlynn implementation. The algorithm generates wrong Left Poisson tail values for lambdas in interval [25; 40]. I send you a table of the expected values (L', R'), and the values generated by MRMC (L, R) with epsilon = 1.0e-7.
L' R' L R 1 0 173 0 173 2 0 174 0 174 3 0 175 0 175 4 0 176 0 176 5 0 177 0 177 6 0 178 0 178 7 0 179 0 179 8 0 180 0 180 9 0 181 0 181 10 0 182 0 182 11 0 183 0 183 12 0 184 0 184 13 0 185 0 185 14 0 186 0 186 15 0 187 0 187 16 0 188 0 188 17 0 189 0 189 18 0 190 0 190 19 0 191 0 191 20 0 192 0 192 21 0 193 0 193 22 0 194 0 194 23 0 195 0 195 24 0 196 0 196 25 0 197 23 197 26 0 198 24 198 27 0 199 25 199 28 0 200 26 200 29 0 201 27 201 30 0 202 28 202 31 0 203 29 203 32 0 204 30 204 33 0 205 31 205 34 0 206 32 206 35 0 207 33 207 36 0 208 34 208 37 0 209 35 209 38 0 210 36 210 39 0 211 37 211 40 0 212 38 212 41 0 213 0 213 42 0 214 0 214 43 1 215 1 215 44 1 216 1 216 45 2 217 2 217 46 3 218 3 218 47 3 219 3 219 48 4 220 4 220 49 4 221 4 221 50 5 222 5 222
I send you also a PDF with an Excel plot of these values. As you can see, the left tail is computed wrongly with lambda=[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.