Gasdev(), Ran1() 함수 설명

gasdev(), ran1() 함수는 어떤때 쓰이는 어떻게 이용하는 것이가?
(참고 :  함수에 있느는 값이 어떤 것을 의미하고 어떻게 작동하는지에 대한 원리는 모른다. 다만 사용법을 설명한다.)

통계에서 정규분포 N( , ) 형식을 보았을 것이다. 우리는 gasdev()와 ran1() 함수를 이용하여 정규분포 N(10, 5^2)을 따르는 모집단을 만들 것이다.
(10, 5^2)에서 ^2는 2의 제곱을 듯한다. 즉, (10, 25)라는 것이다. 2^3 = 2*2*2 //설마 이정도는 알겠지;;;;;;

쓰는 방법은 간단하다.

먼저 ran1()함수에 쓰이는 #define IM 2147483647 값을 메인 함수에서 아래와 같이 선언해 준다.
long idum=IM;  //long 타입 선언, idum는 자신이 마음에 드는 것으로 바꿔도 무방.

이제 gasdev()함수에 IM값을 선언한 idum을 넣어준다. (만약 idum을 다른것으로 바꿨다면 그 값을 넣어준다. 물론 앞에 주소를 알려주는 &기호를  붙이고....)

// 정규분포 N(10, 5^2)을 따르는 1개의 값을 생성

5*gasdev(&idum)+10; 

위처럼 gasdev에 5를 곱하면 5^2를 따르는 분포가 형성되고 뒤에 10을 더하면 평균 10을 가지는 정규분포 N(10, 5^2) 1개의 값을 가지게 된다.
이제 모집단 50개를 만든다면 이 작업을 50번 반복해서 배열에 넣으면 되는 것이다.

정규분포 N(10, 5^2)을 따르는 모집단 10000개를 생성하고 평균과 분포를 구하여라. (대학교 통계 과목 과제중 일부분의 문제입니다.)

작업환경 : C 언어, VS2008


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

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

//아래 define는 내가 직접 설정
#define POP 10000 //모집단의 크기 설정
#define COUNT 5000 //표본추출 횟수 설정


float ran1(long *idum);
float gasdev(long *idum);




int main()
{
float population[POP] = {0};//모집단 만개
float sum=0; //합
float avg=0; //평균
float dispersion=0; //분산
int i=0;
float tmp=0; //분산용 (임시)
long idum=IM;




////////////모집단(POP) 크기에 따른 값 생성////////////
while (i<POP)
{
population[i]=5*gasdev(&idum)+10; // 정규분포 N(10, 5^2)을 따르는 10000개의 값 생성
i++;
}

i=0;
/////////////모집단이 가진 값들의 평균값 구하기//////////
while(i<POP)
{
sum+=population[i];
i++;
}

i=0;
avg=sum/POP; //모집단의 평균값

/////////////모집단이 가진 분산값 구하기//////////////////
while(i<POP)
{
tmp+=(population[i]-avg)*(population[i]-avg); //(Xi(모집단이가진 값-u(평균))^2 값
i++;
}
dispersion=tmp/POP; //분산값

/////////모집단의 크기, 평균, 분산 화면에 출력////////////
printf("-----------모집단 정보--------------\n");
printf("크기 : %d\n평균 : %f\n분산 : %f\n", POP, avg, dispersion);
printf("------------------------------------\n");


return 0;

}


float ran1(long *idum)
//"Minimal" random number generator of Park and Miller with Bayes-Durham shuffle and added safeguards. Returns a uniform random deviate between 0.0 and 1.0(excluxive of the endpoint values). Call with idum a negative integer to initialize; thereafter, do not alter idum between successive deviates in a sequence. RNMX should approximate the largest floating value that is less than 1.
{
int j;
long k;
static long iy=0;
static long iv[NTAB];
float temp;

if (*idum <= 0 || !iy) { // Initialize.
if (-(*idum) < 1) *idum=1; // Be sure to prevent idum=0
else *idum = -(*idum);
for (j=NTAB+7;j>=0;j--){ // Load the shuffle table (after 8 warm-ups).
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
if (j < NTAB) iv[j] = *idum;
}
iy=iv[0];
}
k=(*idum)/IQ; // Start here when not initializing.
*idum=IA*(*idum-k*IQ)-IR*k; // Compute idum=(IA*idum) % IM without overflow
if (*idum < 0) *idum += IM; // by Schrage's method.
j=iy/NDIV; // Will be in the range 0..NTAB-1.
iy=iv[j]; // Output previously stored value and refill the
iv[j] = *idum; // shuffle table.
if ((temp=AM*iy) > RNMX) return RNMX; // Beacause users don't expect endpoint values.
else return temp;
}

float gasdev(long *idum)
//Returns a normally distributed deviate with zero mean and unit variance ran1(idum) as the source of uniform deviates.
{
float ran1(long *idum);
static int iset=0;
static float gset;
float fac, rsq, v1, v2;

if (*idum < 0) iset=0; //Reinitialize.
if (iset == 0) { //We don't have an extra deviate handy, so
do {
v1=2.0*ran1(idum)-1.0; //pick two uniform numbers in the square
v2=2.0*ran1(idum)-1.0; // extending from-1 to +1 in each dirction,
rsq=v1*v1+v2*v2; //see if they are in the unit circle,
} while (rsq >= 1.0 || rsq == 0.0); //and if they are not, try again.
fac=sqrt(-2.0*log(rsq)/rsq);
//Now make the Box-Muller transformation to get two normal deviates.
//Return one and save the other for next time.
gset=v1*fac;
iset=1; //Set flag.
return v2*fac;
} else { //We have an extra deviate handy,
iset=0; //so unset the flag,
return gset; //and return it.
}
}