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.
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.
lambda 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 51 5 223 5 223 52 6 224 6 224 53 7 225 7 225 54 7 226 7 226 55 8 227 8 227 56 8 228 8 228 57 9 229 9 229 58 10 230 10 230 59 10 231 10 231 60 11 232 11 232 61 11 233 11 233 62 12 234 12 234 63 13 235 13 235 64 13 236 13 236 65 14 237 14 237 66 15 238 15 238 67 15 239 15 239 68 16 240 16 240 69 16 241 16 241 70 17 242 17 242 71 18 243 18 243 72 18 244 18 244 73 19 245 19 245 74 20 246 20 246 75 20 247 20 247
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.
mrmc-users@lists.rwth-aachen.de
-
Elvio Amparore