作者PanJC ((#‵Д′)f〒﹌﹌﹌﹌﹌ꄩ
看板NCTU-STAT100
标题[闲聊] Bivariate normal (using box-muller)
时间Thu Jul 28 09:17:05 2011
#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