Why is Mean Square error of a gamma variable not being calculated?

Question

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?

Answer 1

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