So I have written this code that calculates mean square estimates of alpha and beta for different values of samples and betas. The value of alpha is 3, since we add three exponential distributions to create the gamma variable.
However, this code isn't giving me any correct output.
# include <stdio.h>
# include <math.h>
# include <stdlib.h> // Don't forget to include this
main()
{
int n,N; // Sample size and number of simulations
long double alpha=3,beta,suma=0.0,sumb=0.0,msea,mseb;
int i,j,k;
printf("Enter the number of simulations");
scanf("%d", &N);
printf("Enter the sample size");
scanf("%d", &n);
printf("Enter the value of beta");
scanf("%Lf", &beta);
//Simulation starts to calculate MSE
for(i=0;i<N;i++)
{
float msea=0.0,mseb=0.0,sum=0.0,sumsq=0.0; //Vlaue is reset at every iteration, so declared inside i loop
for(j=0;j<n;j++)//each sample
{
float g;
for(k=0;k<alpha;k++)
g += (double)(-1/beta)*log(1-((double)rand()/RAND_MAX));//sum of exp is gamma
sum += g;
sumsq += g*g;
}
float xbar = sum/n;
float var = sumsq/n - xbar*xbar;
suma += pow ((xbar*xbar/var - alpha),2);
sumb += pow ((xbar/var - beta),2);
}
msea = suma/n;
mseb = sumb/n;
printf("MSE of alpha is = %Lf", msea);
printf("MSE of beta is = %Lf", mseb);
return 0;
}
Can you help me find my error?
Code has at least these problems:
Uninitialized g
float g;
for (k = 0; k < alpha; k++)
g += (double) (-1 / beta) * log(1 - ((double) rand() / RAND_MAX));
log(0)
log(1-((double)rand()/RAND_MAX))
might deliver log(0)
. I suspect log(1 - (((double) rand() + 0.5)/ (RAND_MAX + 1u)))
will provide better distribution.
3 FP types
There is a lot of up converting, down converting going on with float, double, long double
.
Use double
throughout until it is determined other widths are needed.
Unused objects
float msea=0.0,mseb=0.0
are not used. Tip: Enable all warnings and save time.
Use %g
"%g"
is more informative.
// printf("MSE of alpha is = %Lf", msea);
printf("MSE of alpha is = %Lg", msea);
OP reports "I made the changes, it's still giving nan output ". I get
MSE of alpha is = 105.474
MSE of beta is = 31.4536
I suspect OP has not made all the changes.
# include <stdio.h>
# include <math.h>
# include <stdlib.h> // Don't forget to include this
int main() {
int n, NN; // Sample size and number of simulations
double alpha = 3, beta, suma = 0.0, sumb = 0.0, msea, mseb;
int i, j, k;
printf("Enter the number of simulations");
// scanf("%d", &NN);
printf("Enter the sample size");
// scanf("%d", &n);
printf("Enter the value of beta");
// scanf("%f", &beta);
NN = 1000;
n = 20;
beta = 1.5;
//Simulation starts to calculate MSE
for (i = 0; i < NN; i++) {
double sum = 0.0, sumsq = 0.0; //Vlaue is reset at every iteration, so declared inside i loop
for (j = 0; j < n; j++) //each sample
{
double g = 0;
for (k = 0; k < alpha; k++)
g += (double) (-1 / beta)
* log(1 - (((double) rand() + 0.5) / (RAND_MAX + 1u)));
sum += g;
sumsq += g * g;
}
double xbar = sum / n;
double var = sumsq / n - xbar * xbar;
suma += pow((xbar * xbar / var - alpha), 2);
sumb += pow((xbar / var - beta), 2);
}
msea = suma / n;
mseb = sumb / n;
printf("MSE of alpha is = %g\n", msea);
printf("MSE of beta is = %g\n", mseb);
return 0;
}