C++生成正态分布样本、M-H采样

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef vector<double> vec;
typedef vector<vec> mat;

const int INF=0x3f3f3f3f;
const int MAX_N=100;
const double PI=acos(-1.0);

double uniform(){//[0,1]均匀分布
    return (double)rand() / RAND_MAX;
}

double gaussrand(double E,double V){//均值E,标准差V
    static double V1, V2, S;
    static int phase = 0;
    double X;

    if ( phase == 0 ) {
        do {
            double U1 = (double)rand() / RAND_MAX;
            double U2 = (double)rand() / RAND_MAX;

            V1 = 2 * U1 - 1;
            V2 = 2 * U2 - 1;
            S = V1 * V1 + V2 * V2;
        } while(S >= 1 || S == 0);
        X = V1 * sqrt(-2 * log(S) / S);
    }
    else{
        X = V2 * sqrt(-2 * log(S) / S);
    }
    phase = 1 - phase;

    X = X * V + E;
    return X;
}

double get_gauss(double x,double E,double V){//正态分布函数,f(x),均值E,标准差V
    return (1/(sqrt(2*PI)*V))*exp(-1*(x-E)*(x-E)/(2*V*V));
}

int main()
{
    srand(time(NULL));
    double pi[10000];
    for(int i=1;i<10000;i++){
        double pi_star=gaussrand(pi[i-1],1);
        double alpha=min(1.0,get_gauss(pi_star,3,2)/get_gauss(pi[i-1],3,2));
        double u=uniform();
        if(u<alpha)pi[i]=pi_star;
        else pi[i]=pi[i-1];
    }
    for(int i=0;i<10000;i++){
        printf("%.12f,",pi[i]);
    }
    return 0;
}