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