Gamma变量均方误差计算异常：代码问题排查求助

以下是你的代码中存在的关键问题，以及对应的修复方法：

1. 未初始化随机数生成器 #

rand()函数在调用前需要用srand()初始化，否则每次运行程序都会生成完全相同的随机数序列，导致模拟结果失去统计意义。

• 修复：在main函数开头添加srand(time(NULL));，同时需要包含<time.h>头文件。

2. 未初始化变量g #

变量float g;声明后未初始化就执行g += ...，会读取内存中的垃圾值，引发未定义行为。

• 修复：在每次生成Gamma变量前，将g初始化为0，即改为float g = 0.0;。

3. 循环条件使用浮点数比较 #

alpha代表指数分布的个数（固定为3），却被定义为long double类型，循环for(k=0;k<alpha;k++)中用整数和浮点数比较，写法不规范且存在潜在风险。

• 修复：将alpha的类型改为int，即int alpha = 3;。

4. MSE计算的除数错误 #

suma和sumb是N次模拟的平方和，计算均方误差时应该除以模拟次数N，而非样本量n。

• 修复：将msea = suma/n;改为msea = suma / N;，mseb = sumb/n;改为mseb = sumb / N;。

5. 冗余且冲突的变量定义 #

在i循环内部定义的float msea=0.0,mseb=0.0完全未被使用，还会屏蔽外层的long double类型同名变量，属于无效代码。

• 修复：直接删除这两个变量的定义。

6. 整数除法导致的计算错误 #

-1/beta中，1是整数类型，会先执行整数除法（结果为0）再转换为浮点型，导致生成的指数分布值完全错误。

• 修复：将-1改为浮点型常量，比如-1.0L，即(-1.0L / beta)。

7. 可能出现的log(0)异常 #

当rand()返回RAND_MAX时，1-((double)rand()/RAND_MAX)会等于0，此时log(0)会产生无穷大值，破坏计算结果。

• 修复：将随机数生成改为rand()/(RAND_MAX + 1.0)，确保得到[0,1)区间的数值，1减去后得到(0,1]，避免log(0)的情况。

修复后的完整代码 #

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>

int main()
{
    int n, N; // Sample size and number of simulations
    int alpha = 3;
    long double beta, suma = 0.0, sumb = 0.0, msea, mseb;
    int i, j, k;

    srand(time(NULL)); // 初始化随机数生成器

    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++)
    {
        double sum = 0.0, sumsq = 0.0; // 使用double提升精度
        for(j = 0; j < n; j++) // each sample
        {
            double g = 0.0;
            for(k = 0; k < alpha; k++)
                g += (-1.0L / beta) * log(1 - ((double)rand() / (RAND_MAX + 1.0))); // 避免log(0)
            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 = %Lf\n", msea);
    printf("MSE of beta is = %Lf\n", mseb);

    return 0;
}

内容的提问来源于stack exchange，提问作者Anweshan Goswami