通过C + GSL与通过R的平行模拟

Question

我在研究工作中使用R进行模拟。最近，由于速度原因，我转向C + GSL。为了说明，我首先通过R写一些简单的代码，如下所示。

    n = 10
    nsim = 20
    out = array(0,c(nsim,1))
    set.seed(123) # for reproducibility
    for (i in 1:nsim){
        d = rnorm(n, mean=1, sd=0.1) # generate data
        m = mean(d)
        out[i] = m
    }
    print(out)

然后我尝试使用C + GSL做同样的事情，其中​​读取

    #include <stdio.h>
    #include <stddef.h>
    #include <stdlib.h>
    #include <math.h>
    #include <time.h>
    #include <gsl/gsl_vector.h>
    #include <gsl/gsl_matrix.h>
    #include <gsl/gsl_blas.h>
    #include <gsl/gsl_linalg.h>
    #include <gsl/gsl_rng.h>
    #include <gsl/gsl_randist.h>
    #include <gsl/gsl_statistics.h>

    int main(void){
        const gsl_rng_type * T;
        gsl_rng * r;

        int i, j, n=10, nsim=20;
        double d[n], out[nsim];
        double mean = 1, sigma = 0.1;

        gsl_rng_env_setup();
        T = gsl_rng_default;
        r = gsl_rng_alloc (T);

        for (int i = 0; i < nsim; i++) {
            for (j = 0; j < n; j++) {
                double d[j] = mean + gsl_ran_gaussian(r, sigma);
                // printf (" %f", x);
            }
            double out[i] = gsl_stats_mean(d, 1, n);    
        }   
    printf ("The output is %g\n", out);
    gsl_rng_free (r);
    return 0;
    }

但是关于我的C + GSL代码是错误的。我是C + GSL的新手。有帮助吗？

Answer 1

for (j = 0; j < n; j++) {
    double d[j] = mean ...

这个关键字double现在宣布一个大小为d的双倍的新数组j。那不是你想要的;你已经声明了你的数组：在这里删除关键字double。

同时在double中删除关键字double out[i] = ...

请注意，在printf ("The output is %g\n", out);中out是一个数组，但你不能一次打印一个数组;你必须逐个打印，所以像这样：

for (int i = 0; i < nsim; i++)
    printf ("%g\n", out[i]);

---

So your code could look like (only relevant parts):

    double d[n], out[nsim];
    ...
    for (int i = 0; i < nsim; i++) {
        for (j = 0; j < n; j++) {
            d[j] = mean + gsl_ran_gaussian(r, sigma);
        }
        out[i] = gsl_stats_mean(d, 1, n);    
    }   
    for (int i = 0; i < nsim; i++)
        printf ("%g\n", out[i]);