Coding complex equations into C++

Question

I have here equations from a paper by E. Bradlow. They're for counting the events in a Weibull-distributed data set.

$$\begin{align} \Pr(N(t)=n) &= \sum_{j=n}^\infty{\frac{(-1)^{j+n}(\lambda t^c)^j \alpha_j^n}{\Gamma(cj+1)}} \\ E(N(t)) &= \sum_{n=1}^\infty{\sum_{j=n}^\infty{\frac{n(-1)^{j+n}(\lambda t^c)^j \alpha_j^n}{\Gamma(cj+1)}}} \end{align}$$

Where $Pr(N(t) = n)$ is the probability that the count will be $n$, and $E(N(t))$ is the expected value. Below are the other vectors.

$$\begin{align} \alpha_m^l &= \sum_{m=0}^{l-1}{\frac{\Gamma(cm+1)\Gamma(cl-cm+1)}{\Gamma(m+1)\Gamma(l-m+1)}} \\ \alpha_j^0 &= \frac{\Gamma(cj+1)}{\Gamma(j+1)}, \quad j=0,1,2,\ldots \\ \alpha_j^{n+1} &= \sum_{m=n}^{j-1}{\alpha_m^n \frac{\Gamma(cj-cm+1)}{\Gamma(j-m+1)}} \end{align}$$

I've programmed these into C++, but I'm not getting the right results. I know because when I get the expected value, I get an extremely negative number (too many zeroes to count). I'd like some input on where I got my code wrong?

    #include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

using namespace std;

double alpha_m(int n, double shape){
       double sum = 0;
       for(int m = 0; m < n-1; m += 1){
             //printf("%i < %i\n", m, n-1);
             double num = tgamma((shape*m) + 1) * tgamma((shape*n) - (shape*m) + 1);
//             printf("num %f , shape %f, m %i\n", num, shape, m);
//             system("pause");
             double denom = tgamma(m + 1) * tgamma(n - m + 1);
             //printf("denum %f \n", denom);
             sum += num/denom;
             //printf("sum%f\n", sum);
       }
//       printf("sum%f\n", sum);
//         printf("done");
       return sum;
}

double alpha_j(int n, int j, double shape){
       //printf("n: %i", n);
       //system("pause");
       double sum = 0;
       if(n == 0){
            double result = tgamma((shape*j) + 1)/tgamma(j + 1);
            return result;
       } else if(n >= 1) {
            for(int m = n-1; m < j-1; m += 1){
                  //printf("%i < %i\n", m, j-1);
                  double result = tgamma((shape*j) - (shape*m) + 1)/tgamma(j - m + 1);
//                  printf("result %f\n", result);
                  sum += result;
//                  printf("sum %f\n", sum);
            }
            double alpha = alpha_m(n,shape);
//            printf("sum %f, alpha %f\n", sum, alpha);
            sum = sum*alpha;
//            printf("sum %f\n", sum);
            return sum;
       }
}

int main()
{
/*  double x;


  printf("libm");
  for(x=1.0; x <= 10.0; x+=1.0) {
    printf("%15.8lf\n", tgamma(x/3.0));
  }
  system("pause");
  return 0;
*/
  double scale, shape;
  int freq;

  cout << "Enter scale param: ";
  cin >> scale;
  cout << "Enter shape param: ";
  cin >> shape;
  cout << "Enter Frequency of PM: ";
  cin >> freq;

  printf("Scale: %f\nc: %f\nFrequency: %i\n\n", scale, shape, freq);
  system("pause");


  double sum = 0;

  if(freq > 0){
      for(int n=1; n <= 99; n++){
            //printf("%i\n", n);
            for(int j=n; j<=99; j++){
                  double t = 365/freq;
                  double power = pow (t, shape);
                  double num = n*((-1)^(j+n))* pow (scale*power, j);
                  double denom = tgamma((shape*j) + 1);
                  //printf("%f", alpha_j(n, j, shape));
                  sum += num*alpha_j(n, j, shape)/denom;
//                  printf("t %f, power %f, num %f, denom %f, sum %f\n, alpha %f", t, power, num, denom, sum, alpha_j(n, j, shape));
//                  system("pause");
                  //printf("n = %i ;; j = %i\n", n, j);
            }
      }
  } else if(freq == 0){
         sum = (365*24)/(scale*tgamma((1/shape) + 1));
         //printf("Mean TTF gamma %f", tgamma((1/shape) + 1));
  }

  printf("The final sum is %f\n", sum);
  system("pause");
}

Answer 1

There are several issues that you will need to resolve.

1) In your function alpha_m the loop should be

for (int m = 0; m < n; m++) {
}

because the sum goes from 0 to n-1.

2) I'm not sure what you tried in the function alpha_j. But I believe the loop should not end at m = j-2. Try to match your variable names and starting indices with those of the equation you are trying to code.

3) Also, function alpha_j does not return anything if n < 0.

4) The main routine does not deal with the freq < 0 case.