cgamma-distributionmean-square-error

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


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?


Solution

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