난수(변환법, 거부법, 가우시안 분포, 몬테카를로 적분)

학교수업, CS/수치해석

난수(변환법, 거부법, 가우시안 분포, 몬테카를로 적분)

빨대도둑 2022. 11. 30. 14:47

난수 발생

알고리즘

적당한 크기의 배열을 잡아서 그 안의 내용을 균일한 편차를 갖는 난수로 채운다. srand()로 초기화하고 rand()로 난수를 발생시킨다
난수 발생을 요구할 때마다 배열의 인덱스를 난수로 발생시키고 그 내부값을 리턴하기 전에 한번 사용한 값은 버리고 rand()를 이용하여 새로운 값으로 채운다
2번 과정을 필요한 횟수만큼 반복한다

<수식 1>

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

#define seed 31
#define NUM 100
#define MX 32767

static double table[100], index;

double ran_gen() {
	unsigned int i;
	i = (int)(1 + 99.0 * index / MX);
	if (i > 99 || i < 0) {
		printf("error! \n");
		exit(0);
	}
	index = table[i];
	table[i] = rand();
	return((double)index / MX);
}

int main(void) {
	int i, dummy;
	double a;
	srand(seed);
	for (i = 0; i < 101; i++)
		dummy = rand();
	for (i = 0;i < 101; i++)
		table[i] = rand();
	for (i = 0; i < NUM; i++) {
		a = ran_gen();
		printf("%8.6lf", a);
		if ((i + 1) % 5 == 0)
			printf("\n");
	}
}

<수식 2>

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

#define NUM 100
#define seed -13
#define MX1 167823
#define R1 1.0/167823.0
#define A1 2167
#define C1 20471

#define MX2 298331
#define R2 1.0/298331.0
#define A2 741
#define C2 25897

#define MX3 109283
#define R3 1.0/109283.0
#define A3 4351
#define C3 18761

static long x1, x2, x3;
static double table[100];

double ran_gen() {
	double temp;
	unsigned int i;
	x1 = (A1 * x1 + C1) % MX1;
	x2 = (A2 * x2 + C2) % MX2;
	x3 = (A3 * x3 + C3) % MX3;

	i = 1 + ((99 * x3) * R3);

	if (i > 99 || i < 0) {
		printf("Error \n");
			exit(0);
	}
	temp = table[i];
	table[i] = (x1 + x2 * R2) * R1;
	return temp;
}

int main(void) {
	int i;
	double a;

	x1 = (C1 - seed) % MX1;
	x1 = (A1 * x1 + C1) % MX1;
	x2 = x1 % MX2;
	x1 = (A1 * x1 + C1) % MX1;
	x3 = x1 % MX3;

	for (i = 0; i < 100; i++) {
		x1 = (A1 * x1 + C1) % MX1;
		x1 = (A2 * x2 + C2) % MX2;
		table[i] = (double)(x1 + (double)x2 / MX2) / MX1;
	}
	
	for (i = 0; i < NUM; i++) {
		a = ran_gen();
		printf("%8.6lf", a);
		if ((i + 1) % 5 == 0)
			printf("\n");
	}
	
}

변환법

변환 방법은 균일 확률 분포를 갖는 수열을 원하는 확률분포로 변환시키는 방법이다. 원하는 확률밀도의 적분인 확률 분포 함수를 구하고, 함수의 역함수를 구하여 균일 확률 분포 함수를 변환시키면 원하는 확률분포를 갖는 난수열을 만들 수 있다.

<수식>

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <conio.h>
#include <errno.h>
#pragma warning(disable:4996)//fopen, fclose를 해결하기 위해서 이 방법을 사용함

#pragma warning(suppress : 4996)

#define NUM 100000
#define SEED -13
#define MX1 167823
#define R1 1.0/167823.0
#define A1 2167
#define C1 20471

#define MX2 298331
#define R2 1.0/298331.0
#define A2 741
#define C2 25897

#define MX3 109283
#define R3 1.0/109283.0
#define A3 4351
#define C3 18761

static long x1, x2, x3;
static double table[100];


double ran_gen() {
	double temp;
	unsigned int i;
	static ini = 0;

	if (ini == 0) {
		x1 = (C1 - SEED) % MX1;
		x1 = (A1 * x1 + C1) % MX1;
		x2 = x1 % MX2;
		x1 = (A1 * x1 + C1) % MX1;
		x3 = x1 % MX3;

		for (i = 0; i < 100; i++) {
			x1 = (A1 * x1 + C1) % MX1;
			x2 = (A2 * x2 + C2) % MX2;
			table[i] = (double)(x1 + (double)x2 / MX2) / MX1;
		}
		ini = 1;
	}
	x1 = (A1 * x1 + C1) % MX1;
	x2 = (A2 * x2 + C2) % MX2;
	x3 = (A3 * x3 + C3) % MX3;

	i = 1 + ((99 * x3) * R3);

	if (i > 99 || i < 0) {
		printf("Error \n");
		exit(0);
	}
	temp = table[i];
	table[i] = (x1 + x2 * R2) * R1;
	return temp;
}

int main(void) {
	int i;
	double a;
	FILE* fp1, * fp2;

	fp1 = fopen("exp.dat", "wb");
	if (fp1 == NULL) { exit(0); }
	fp2 = fopen("test.dat", "w");
	if (fp2 == NULL) { exit(0); }

	for (i = 0; i < NUM; i++) {
		a = ran_gen();
		a = 1.0 - a;
		if (a == 0)continue;
		a = -1.0 * log(a);
		fprintf(fp2, "%8.6lf", a);
		if ((i + 1) % 7 == 0)fprintf(fp2, "\n");

		fwrite(&a, 4, 1, fp1);
	}
	fcloseall();

}

거부법

변환법은 구하고자 하는 확률 분포 함수의 역함수를 알아야 한다는 단점을 가지고 있다.

거부법은 역함수를 구하지 않고도 원하는 분포의 난수열을 발생시킬 수 있는 방법이다. 역함수 대신 구하고자 하는 확률 밀도 함수의 위쪽에 있는, 역함수를 아는 임의의 함수를 사용한다. 이 방법도 우선 균일 확률 분포 난수를 발생시키고 이를 적절히 조작하여 주어진 조건을 만족하는 대강의 분포를 갖도록 한 다음에 원하는 분포 내에 들어가는지 확인하여 거부, 선택을 하는 방법이다.

<수식>

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <conio.h>
#include <errno.h>
#pragma warning(disable:4996)//fopen, fclose를 해결하기 위해서 이 방법을 사용함

#pragma warning(suppress : 4996)

#define NUM 100000
#define SEED -13
#define MX1 167823
#define R1 1.0/167823.0
#define A1 2167
#define C1 20471

#define MX2 298331
#define R2 1.0/298331.0
#define A2 741
#define C2 25897

#define MX3 109283
#define R3 1.0/109283.0
#define A3 4351
#define C3 18761

static long x1, x2, x3;
static double table[100];


double ran_gen() {
	double temp;
	unsigned int i;
	static ini = 0;

	if (ini == 0) {
		x1 = (C1 - SEED) % MX1;
		x1 = (A1 * x1 + C1) % MX1;
		x2 = x1 % MX2;
		x1 = (A1 * x1 + C1) % MX1;
		x3 = x1 % MX3;

		for (i = 0; i < 100; i++) {
			x1 = (A1 * x1 + C1) % MX1;
			x2 = (A2 * x2 + C2) % MX2;
			table[i] = (double)(x1 + (double)x2 / MX2) / MX1;
		}
		ini = 1;
	}
	x1 = (A1 * x1 + C1) % MX1;
	x2 = (A2 * x2 + C2) % MX2;
	x3 = (A3 * x3 + C3) % MX3;

	i = 1 + ((99 * x3) * R3);

	if (i > 99 || i < 0) {
		printf("Error \n");
		exit(0);
	}
	temp = table[i];
	table[i] = (x1 + x2 * R2) * R1;
	return temp;
}

double gam(double a) {
	int j;
	double am, e, s, v1, v2, x, y;

	if (a < 1.0) {
		printf("a should be greater than 1 \n");
		exit(0);
	}
	if (a < 6.0)/*direct methd, adding wating times*/
	{
		x = 1.0;
		for (j = 1; j < a + 1; j++)
			x *= ran_gen();
		x = -1.0 * log(x);
	}
	else {
		do {
			do {
				do {
					v1 = 2.0 * ran_gen() - 1.0;
					v2 = 2.0 * ran_gen() - 1.0;
				} while (v1 * v1 + v2 * v2 > 1.0);
				y = v2 / v1;
				am = a - 1;
				s = sqrt(2.0 * am - 1.0);
				x = s * y + am;
			} while (x <= 0.0);
			e = (1.0 + y * y) * exp(am * log(x / am) - x * y);
		} while (ran_gen() > e);
	}return x;
}

int main(void) {
	int i;
	double a;
	FILE* fp1;

	fp1 = fopen("gam.dat", "wb");
	if (fp1 == NULL) { exit(0); }

	for (i = 0; i < NUM; i++) {
		a = gam(3.0);
		fwrite(&a, 4, 1, fp1);
	}
	fcloseall();

}

가우시안 분포

가우시안 분포는 자연계에서 가증 흔하게 볼 수 있는 분포이다. 가우시안 확률 밀도 함수는 수학적으로도 매우 우수한 성질을 가질 뿐만 아니라 실험적으로도 많이 사용한다.

가우시안 분포를 갖는 난수를 발생시키기 위해서 중간 극한 정리를 이용한 방법을 사용했다. 중간 극한 정리의 엄밀한 정의와 증명은 통계학이나 불규칙 신호론을 참고하면된다.

<수식>

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <conio.h>
#include <errno.h>
#pragma warning(disable:4996)//fopen, fclose를 해결하기 위해서 이 방법을 사용함

#pragma warning(suppress : 4996)

#define NUM 100000
#define SEED -13
#define MX1 167823
#define R1 1.0/167823.0
#define A1 2167
#define C1 20471

#define MX2 298331
#define R2 1.0/298331.0
#define A2 741
#define C2 25897

#define MX3 109283
#define R3 1.0/109283.0
#define A3 4351
#define C3 18761

#define NS 20
#define MEAN 0.5

static long x1, x2, x3;
static double table[100];


double ran_gen() {
	double temp;
	unsigned int i;
	static ini = 0;

	if (ini == 0) {
		x1 = (C1 - SEED) % MX1;
		x1 = (A1 * x1 + C1) % MX1;
		x2 = x1 % MX2;
		x1 = (A1 * x1 + C1) % MX1;
		x3 = x1 % MX3;

		for (i = 0; i < 100; i++) {
			x1 = (A1 * x1 + C1) % MX1;
			x2 = (A2 * x2 + C2) % MX2;
			table[i] = (double)(x1 + (double)x2 / MX2) / MX1;
		}
		ini = 1;
	}
	x1 = (A1 * x1 + C1) % MX1;
	x2 = (A2 * x2 + C2) % MX2;
	x3 = (A3 * x3 + C3) % MX3;

	i = 1 + ((99 * x3) * R3);

	if (i > 99 || i < 0) {
		printf("Error \n");
		exit(0);
	}
	temp = table[i];
	table[i] = (x1 + x2 * R2) * R1;
	return temp;
}


int main(void) {
	int i, j;
	double a, ss, s, sig, sum_a;
	FILE* fp1;

	fp1 = fopen("gau.dat", "wb");
	if (fp1 == NULL) { exit(0); }

	ss = s = 0.0;
	for (i = 0; i < NUM; i++) {
		sum_a = 0.0;
		for (j = 0; j < NS; j++) {
			a = ran_gen();
			sum_a += a;
		}
		sig = sqrt(NS);
		a = sum_a;

		ss += a * a;
		s += a;
		fwrite(&a, 4, 1, fp1);
	}
	printf("MEAN : %lf \n VAR : %lf\n", s/NUM, ss/NUM - s/NUM * s/NUM);
	fcloseall();
}

몬테카를로 적분

컴퓨터를 이용하여 어떤 함수를 적분할 때 주로 쓰이는 방법이다

<수식>

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <conio.h>
#include <errno.h>
#include <time.h>
#pragma warning(disable:4996)//fopen, fclose를 해결하기 위해서 이 방법을 사용함
#define randomize() srand((unsigned)time(NULL))

#pragma warning(suppress : 4996)

#define NUM 100
#define SEED -13
#define MX1 167823

#define R1 1.0/167823.0
#define A1 2167
#define C1 20471

#define MX2 298331
#define R2 1.0/298331.0
#define A2 741
#define C2 25897

#define MX3 109283
#define R3 1.0/109283.0
#define A3 4351
#define C3 18761

#define NS 20
#define MEAN 0.5

static long x1, x2, x3;
static double table[100];

double fun(double xx) {
	double v;
	v = xx * xx * xx;
	return v;
}


double ran_gen() {
	double temp;
	unsigned int i;
	static ini = 0;

	if (ini == 0) {
		x1 = (C1 - SEED) % MX1;
		x1 = (A1 * x1 + C1) % MX1;
		x2 = x1 % MX2;
		x1 = (A1 * x1 + C1) % MX1;
		x3 = x1 % MX3;

		for (i = 0; i < 100; i++) {
			x1 = (A1 * x1 + C1) % MX1;
			x2 = (A2 * x2 + C2) % MX2;
			table[i] = (double)(x1 + (double)x2 / MX2) / MX1;
		}
		ini = 1;
	}
	x1 = (A1 * x1 + C1) % MX1;
	x2 = (A2 * x2 + C2) % MX2;
	x3 = (A3 * x3 + C3) % MX3;

	i = 1 + ((99 * x3) * R3);

	if (i > 99 || i < 0) {
		printf("Error \n");
		exit(0);
	}
	temp = table[i];
	table[i] = (x1 + x2 * R2) * R1;
	return temp;
}


int main(void) {
	int i, num;
	double lb, ub, xx, fx, intv, ran, est, s2, sum, ssum;

	printf("\n Type the number of samples : ");
	scanf_s("%d", &num);

	printf("\n Type the lower limit of the integral : ");
	scanf_s("%lf", &lb);

	printf("\n Type the upper limit of the integral : ");
	scanf_s("%lf", &ub);

	intv = ub - lb;
	sum = ssum = 0.0;

	randomize();

	for (i = 0; i < NUM; i++) {
		xx = lb + ran_gen() * intv;
		fx = fun(xx);
		sum += fx;
		ssum += (fx * fx);
	}
	est = intv * sum / (double)num;
	s2 = (ssum * intv * intv - est * est * (double)num) / (double)(num - 1.0);
	s2 = sqrt(s2 / (double)num);

	printf("\n the estimated intrgral value if  %lf \n", est);
	printf("it can be stated with a confidence of 95 percent \n");
	printf("that the true integral value exists in this boundary. \n");
	printf("\n %lf <I < %lf", est - 2.0 * s2, est + 2.0 * s2);

}

'학교수업, CS > 수치해석' 카테고리의 다른 글

수치미분(테일러시리즈, 리차드슨, 직사각형 수치미분, 사다리꼴 적분, 심프슨의 수치적분) (0)	2022.11.30
푸리에 변환(DFT, FFT, 파워스펙트럼 추정, 자기상관도) (0)	2022.11.30
보간법(선형, 랑그라제, 네빌레, 뉴튼 다항식, 전/후향 차분법, 스플라인) (0)	2022.11.30
최소 자승법(촐레스키법을 이용한 행렬의 제곱근) (0)	2022.11.30
고유치(파데브-레브리어의 알고리즘) (0)	2022.11.30

현재글난수(변환법, 거부법, 가우시안 분포, 몬테카를로 적분)

brilliant-star