# Coding complex equations into C++

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");
}

• What have you already done to debug your problem? – Wolfgang Bangerth Mar 17 '14 at 1:09
• @markovchain: Welcome to SciComp! People volunteer their time to answer questions here, and debugging questions tend either to involve a lot of effort to answer in detail, or to devolve into stock answers (use good style, make your code modular, unit test it against inputs with known outputs, step through it with a debugger, etc.). In the interest of time, please tell us what you've done (as far as testing, etc.), and the information you've gotten out of it. If you haven't done any testing, try it, see how far you get, and then add what you've learned from the process to your question. – Geoff Oxberry Mar 17 '14 at 7:59
• To add to @GeoffOxberry comment, there are many unit test frameworks with some of the popular ones being Google Test and Boost Test. Both are cross platform and easy to use. – Damien Mar 17 '14 at 10:58

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.