- int n = int(3000.0 / int_step_size + 0.5);
- double h = 3000.0 / double(n);
- double sum = evaluate_int_point(k, a, prodai_precompute, kfac_precompute, r1, mu2, sigma2, winfac, 0.0);
-
- for (int i = 1; i < n; i += 2) {
- sum += 4.0 * evaluate_int_point(k, a, prodai_precompute, kfac_precompute, r1, mu2, sigma2, winfac, i * h);
- }
- for (int i = 2; i < n; i += 2) {
- sum += 2.0 * evaluate_int_point(k, a, prodai_precompute, kfac_precompute, r1, mu2, sigma2, winfac, i * h);
- }
- sum += evaluate_int_point(k, a, prodai_precompute, kfac_precompute, r1, mu2, sigma2, winfac, 3000.0);
-
- return (h/3.0) * sum;
-}
-
-static inline double evaluate_int_point(int k, double a, double prodai_precompute, double kfac_precompute, double r1, double mu2, double sigma2, double winfac, double x)
-{
- double probscore = prob_score_real(k, a, prodai_precompute, kfac_precompute, (x - r1)*winfac);
- double z = (x - mu2)/sigma2;
- double gaussian = exp(-(z*z/2.0));
- return probscore * gaussian;
+ return simpson_integrate(ProbScoreEvaluator(k, a, prodai_precompute, kfac_precompute, r1, mu2, sigma2, winfac), 0.0, 3000.0, int_step_size);