WYL's Homepage
0%

C语言蒙特卡罗

发表于 更新于 分类于
本文字数： 1.4k 阅读时长 ≈ 5 分钟

2023.9.19更新：学完统计回头来看，这只不过是一点非常fundamental的技术，似乎不值得专门挂一篇blog，但毕竟这小东西还挺别致，当时也让我着迷了一晚上，所以这篇就保留着吧。

今天看了《Statistical Data Analysis》一书中关于Monte Carlo方法的介绍，觉得非常有趣，遂动手用C将其实现，源代码在本篇blog的结尾处。

生成的Maxwell速率分布

程序整体由三个部分组成，分别是生成单位随机数、Transformation method和Acceptance Rejection method[1]。

生成单位随机数所使用的是Linear congruential generator算法，由于取的modulus比较大，所以能够实现比rand()函数更加精细的生成结果。最初的种子来自系统时间，这就保证了每次生成的结果是不同的。生成的随机数必然是整数，所以要将其单位化，之后还要根据实际需要将其进行平移和伸缩。

Transformation方法的概述如下：设$F(x)$为随机变量$x$的累积分布函数，如果我们取$y=F(x),F(x)\in [0,1]$作为一个新的随机变量，那么$\frac{\text{d}y}{\text{d}x}=F’(x)=f(x)$，所以随机变量$y$所满足的pdf(probability density function)函数为

$$g(y)=f(x)|\frac{\text{d}x}{\text{d}y}|=1$$

也就是$y$实际上服从$[0,1]$上的均匀分布。根据这一点，我们可以反向生成目标分布，因为均匀分布是非常容易实现的。现在我们将$x$视为$y$的函数，那么我们有

$$F(x(y))=\int_{-\infty}^{x(y)}f(t)\text{d}t=y$$

将这个式子积分就得到了$x(y)$和$y$的函数，从中可以（对于比较简单的函数）将$x(y)$反解出来。利用这个关系和已经生成好的$y$，就可以生成服从目标分布的$x$。

Transformation的概念非常清晰，但缺点也是显而易见的：当函数比较复杂的时候反解$x(y)$的工作几乎是不可能的。这时我们可以利用另一种方法——Acceptance Rejection method。

这种方法的步骤如下：（1）生成位于$[x_{min},x_{max}]$的均匀分布的$x$序列；（2）生成位于$[0,f_{max}]$的均匀分布的$u$序列；（3）如果$u<f(x)$，则接受$x$，否则就拒绝$x$。由于在$x$处，$x$被接受的概率正比于$f(x)$，所以得到的分布也是服从$f(x)$的。

然而这种方法也有其缺点：当$f(x)$有尖峰和长尾的时候，效率很低下，往往发生的情况是生成了20个$x$，也不一定有一个能被接受。所以我们可以对这种方法进行进一步的改进：将初始的矩形域修改为小一点的域。也就是首先寻找一个能将$f(x)$包括进去的曲线$g(x)$，然后（1）生成一个服从$\frac{g(x)}{\int g(t)\text{d}t}$的$x$；（2）生成一个服从$[0,g(x)]$均匀分布的$u$；（3）如果$u<f(x)$，则接受$x$，重复以上步骤。但是由于目前缺少关于这种方法的实例，我也没有迫切的改进效率的需求，我暂时没有在程序中包括这种改进后的方法。

以下是程序的源代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
#include<time.h>

/*定义multiplier和modulus*/
#define M 4294967296 // 2^32
#define A 1664525
#define C 1013904223
#define PI 3.14159265358979323846

/*定义整数和浮点数两种顺序表*/
struct intseq {
	int MAXNUM;
	unsigned int* element;
}; typedef struct intseq* Pintseq;
struct floatseq {
	int MAXNUM;
	float* element;
}; typedef struct floatseq* Pfloatseq;

/*建立空表*/
Pintseq createintseq(int m) {
	Pintseq pis = (Pintseq)malloc(sizeof(struct intseq));
	if (pis != NULL) {
		pis->element = (unsigned int*)malloc(sizeof(unsigned int) * m);
		if (pis->element) {
			pis->MAXNUM = m;
			return pis;
		}
		else free(pis);
	}
	printf("Out of space!"); return NULL;
}
Pfloatseq createfloatseq(int m) {
	Pfloatseq pfs = (Pfloatseq)malloc(sizeof(struct floatseq));
	if (pfs != NULL) {
		pfs->element = (float*)malloc(sizeof(float) * m);
		if (pfs->element) {
			pfs->MAXNUM = m;
			return pfs;
		}
		else free(pfs);
	}
	printf("Out of space!"); return NULL;
}

/*新建随机数（原始整数）顺序表，并返回指针*/
Pintseq Random(int size, unsigned see) {
	int i;
	unsigned int seed;
	if (see == 0) { seed = (unsigned)time(NULL); }
	else { seed = (unsigned)time(NULL); seed = seed * see; }
	Pintseq pis = createintseq(size);
	pis->element[0] = (A * seed + C) % M;
	for (i = 1; i < pis->MAXNUM; i++) {
		pis->element[i] = (A * pis->element[i - 1] + C) % M;
	}
	return pis;
}

/*输出（打印至屏幕）*/
void printintseq(Pintseq pis) {
	int i;
	for (i = 0; i < pis->MAXNUM; i++) {
		printf("%d\n", pis->element[i]);
	}
	return;
}
void printfloatseq(Pfloatseq pfs) {
	int i;
	for (i = 0; i < pfs->MAXNUM; i++) {
		printf("%f\n", pfs->element[i]);
	}
	return;
}

/*单位化原始随机数*/
Pfloatseq unifying(Pintseq pis) {
	Pfloatseq pfs = createfloatseq(pis->MAXNUM);
	int i;
	for (i = 0; i < pfs->MAXNUM; i++) {
		pfs->element[i] = (float)pis->element[i] / M;
	}
	return pfs;
}

/*对随机数进行transform*/
void transform(Pfloatseq pfs) {
	int i;
	for (i = 0; i < pfs->MAXNUM; i++) {
		pfs->element[i] = -log(pfs->element[i]);
	}
	return;
}

/*对随机数进行区间的平移和伸缩变换*/
void moveandscale(Pfloatseq pfs, float xmin, float xmax) {
	int i;
	for (i = 0; i < pfs->MAXNUM; i++) {
		pfs->element[i] = xmin + pfs->element[i] * (xmax - xmin);
	}
	return;
}

/*进行acceptance-rejection筛选，并输出结果*/
void acceptance_rejection(Pfloatseq pxs, Pfloatseq pus) {
	int i; float x, fxi;
	for (i = 0; i < pxs->MAXNUM; i++) {
		/*在这里指定f(x)的具体形式*/
		x = pxs->element[i];
		fxi = 4 / sqrt(PI) * exp(-pow(x, 2)) * pow(x, 2);

		if (pus->element[i] < fxi) {
			printf("%f\n", pxs->element[i]);
		}
	}
	return;
}

int main() {
	/*数据规模*/
	printf("Please type the numbers of the random numbers you want:\n");
	int size; scanf("%d", &size);

	/*生成标准化的x和u数列*/
	Pintseq pis = Random(size, 2);
	Pfloatseq pxs = unifying(pis);
	pis = Random(pxs->MAXNUM, 7);
	Pfloatseq pus = unifying(pis);

	/*伸缩变换*/
	moveandscale(pxs, 0, 5);
	moveandscale(pus, 0, 4 / sqrt(PI) * exp(-1));

	/*acceptance-rejection*/
	acceptance_rejection(pxs, pus);
}