[閒聊] Bivariate normal (using box-muller)

[PanJC ((#‵Д′)f〒﹌﹌﹌﹌﹌ꄩ]

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.141592653589793

void gen2Dnorm(double *nrm1, double *nrm2, int n);
double mean(double *arr, int n);
double variance(double *arr, int n, double sum);
double covariance(double *arr1, double *arr2, int n, double m1, double m2);

int main(){
        srand(10);
        double *norm1, *norm2;
        int n;
        printf("Input the sample size: ");
        scanf("%d",&n);
        norm1 = (double *) calloc(n,sizeof(double));
        norm2 = (double *) calloc(n,sizeof(double));
        gen2Dnorm(norm1, norm2, n);

        double m1, m2, v1, v2, cov;
        m1 = mean(norm1,n);
        m2 = mean(norm2,n);
        v1 = variance(norm1,n,m1);
        v2 = variance(norm2,n,m2);
        cov= covariance(norm1, norm2, n, m1, m2);
        printf("means: %lf\t%.4lf\n",m1, m2);
        printf("cov_matrix: %.4lf\t%.4lf\n",v1, cov);
        printf("            %.4lf\t%.4lf\n",cov, v2);
        free(norm1);
        free(norm2);
}

void gen2Dnorm(double *nrm1, double *nrm2, int n){
// input
//  nrm1: double 陣列指標名稱,
//  nrm2: double 陣列指標名稱,
//  n   : 陣列長度
//
// return: 無 return
// object: 產生一組二元標準常態分配
        int i;
        double u1, u2;
        for( i = 0 ; i < n ; i++ ){
                u1 =  (double) rand() / RAND_MAX;
                u2 =  (double) rand() / RAND_MAX;
                nrm1[i] = sqrt(-2.0*log(u1)) * cos(2.0*PI*u2);
                nrm2[i] = sqrt(-2.0*log(u1)) * sin(2.0*PI*u2);
        }
}

double mean(double *arr, int n){
        int i;
        double sum=0.0;
        for( i = 0 ; i < n ; i++ ){
                sum += arr[i];
        }
        return(sum/n);
}

double variance(double *arr, int n, double sum){
        int i;
        double tmp=0.0;
        for( i = 0 ; i < n ; i++ ){
                tmp += pow(arr[i],2);
        }
        return(tmp/n - pow(sum,2));
}

double covariance(double *arr1, double *arr2, int n, double m1, double m2){
        double tmp=0.0;
        int i;
        for( i = 0 ; i < n ; i++ ){
                tmp += arr1[i]*arr2[i];
        }
        return(tmp/n-m1*m2);
}

--
※ 發信站: 批踢踢實業坊(ptt.cc) 
◆ From: 140.113.7.248