门特卡罗方法计算渗透模型中的p阈值
问题：
研究人员对以下问题感兴趣：如果site为open的概率为p（因此为blocked的概率为1- p），那系统是渗透系统的概率是多少？
框架
为了实现蒙特卡洛方法计算得到渗透模型的p值，我们不妨首先搭好框架如下：
public static void main(String[] args) {
double result = 0;
System.out.println("Net_Size == "+NET_SIZE);
for (int i = 0; i < RUNTIMES; i++) {
double temp = motokaro_p();
result += temp;
//System.out.println("第"+i+"次的p值："+temp);
}
System.out.println((double) result / RUNTIMES);
}
其中，NET_SIZE为网格的大小，RUNTIMES为进行蒙特卡洛模拟的次数，将每次运行得到的p值保存在temp中，再累加到cnt上，最好除以总次数即为p的数学期望。
蒙特卡洛实现
在蒙特卡洛方法的具体实现中，由于每次模拟开始前，都会预先对网格的渗透进行判断，因此cnt的初始值置为1。网格数组预先置为CLOSED，由于当所有网格都被打开时，网格必定渗透，因此得到for循环的结束条件。
在每次循环中，应用了如下算法：
- 每次循环前对网格的渗透情况进行判断，若已经渗透，则结束循环，计算得到本次的p值作为返回值；
- 如果当前获得的随机方块之前已经打开，则跳过本次循环，同时计数减一来还原循环条件；
- 激活该方块。
其具体实现如下：
public static double motokaro_p() {
int cnt = 1;
Arrays.fill(motokaro_NET, CLOSED);
Random r = new Random();
for (int i = 0; i < NET_SIZE * NET_SIZE; i++) {
if (!is_unblocked()) {
int x = r.nextInt(NET_SIZE);
int y = r.nextInt(NET_SIZE);
if (motokaro_NET[x * NET_SIZE + y] != CLOSED) {
i--;
continue;
}
cnt++;
pour_into(x, y);
//print_net();
} else break;
}
//System.out.println(cnt);
return (double) cnt / (NET_SIZE * NET_SIZE);
}
方块激活方法
对于获取到的方块坐标：
- 若其为第一排方块，直接置为开放；
- 否则对方块的情况使用
check_all方法进行检查。
private static void pour_into(int x, int y) {
motokaro_NET[x * NET_SIZE + y] = OPENED;
if (x == 0) {
motokaro_NET[x * NET_SIZE + y] = UNBLOCKED;
}
check_all(x, y);
return;
}
在check_all方法中，依据如下步骤来激活该方块和其周围的方块：
- 若该方块为打开状态，当其周围存在方块为已渗透状态，则该方块也置为已渗透状态；
- 若该方块为已渗透状态，若其周围存在方块为打开状态，则置该邻居方块为已渗透状态，再利用
check_all方法进行递归检查。
这样的算法即为类比水中放入一枚石子，其激发的涟漪会改变周围水域的状态，直到遇到障碍停住；其具体实现如下：
private static void check_all(int x, int y) {
if (motokaro_NET[x * NET_SIZE + y] == OPENED && ((x > 0 && motokaro_NET[(x - 1) * NET_SIZE + y] == UNBLOCKED) || (x < NET_SIZE - 1 && motokaro_NET[(x + 1) * NET_SIZE + y] == UNBLOCKED)
|| (y > 0 && motokaro_NET[x * NET_SIZE + y - 1] == UNBLOCKED) || (y < NET_SIZE - 1 && motokaro_NET[x * NET_SIZE + y + 1] == UNBLOCKED)))
motokaro_NET[x * NET_SIZE + y] = UNBLOCKED;
if (motokaro_NET[x * NET_SIZE + y] == UNBLOCKED) {
if (x > 0 && motokaro_NET[(x - 1) * NET_SIZE + y] == OPENED) {
motokaro_NET[(x - 1) * NET_SIZE + y] = UNBLOCKED;
check_all(x - 1, y);
}
if (x < NET_SIZE - 1 && motokaro_NET[(x + 1) * NET_SIZE + y] == OPENED) {
motokaro_NET[(x + 1) * NET_SIZE + y] = UNBLOCKED;
check_all(x + 1, y);
}
if (y > 0 && motokaro_NET[x * NET_SIZE + y - 1] == OPENED) {
motokaro_NET[x * NET_SIZE + y - 1] = UNBLOCKED;
check_all(x, y - 1);
}
if (y < NET_SIZE - 1 && motokaro_NET[x * NET_SIZE + y + 1] == OPENED) {
motokaro_NET[x * NET_SIZE + y + 1] = UNBLOCKED;
check_all(x, y + 1);
}
}
return;
}
网格渗透判断方法
网格渗透判断方法相对简单，直接判断最底层的方块是否存在已渗透状态即可。
private static boolean is_unblocked() {
for (int i = 0; i < NET_SIZE; i++)
if (motokaro_NET[(NET_SIZE - 1) * NET_SIZE + i] == UNBLOCKED)
return true;
return false;
}
其他辅助方法的说明
除此之外还调用了一个方法来打印网格的情况以便debug：
private static void print_net() {
System.out.println("NET:");
for (int i = 0; i < NET_SIZE; i++) {
for (int j = 0; j < NET_SIZE; j++) {
System.out.print(String.format(" %2d", motokaro_NET[i * NET_SIZE + j]));
}
System.out.print("\n");
}
}
总体代码
全部代码如下：
package com.company;
import java.util.Random;
import java.util.Arrays;
public class Main {
public final static int CLOSED = -1;
public final static int OPENED = 0;
public final static int UNBLOCKED = 1;
public final static int RUNTIMES = 1000;
public final static int NET_SIZE = 100;
private static int[] motokaro_NET = new int[NET_SIZE * NET_SIZE + 1];
private static int[] NET_visited = new int[NET_SIZE * NET_SIZE + 1];
private static void print_net() {
System.out.println("NET:");
for (int i = 0; i < NET_SIZE; i++) {
for (int j = 0; j < NET_SIZE; j++) {
System.out.print(String.format(" %2d", motokaro_NET[i * NET_SIZE + j]));
}
System.out.print("\n");
}
}
private static void check_all(int x, int y) {
if (motokaro_NET[x * NET_SIZE + y] == OPENED && ((x > 0 && motokaro_NET[(x - 1) * NET_SIZE + y] == UNBLOCKED) || (x < NET_SIZE - 1 && motokaro_NET[(x + 1) * NET_SIZE + y] == UNBLOCKED)
|| (y > 0 && motokaro_NET[x * NET_SIZE + y - 1] == UNBLOCKED) || (y < NET_SIZE - 1 && motokaro_NET[x * NET_SIZE + y + 1] == UNBLOCKED)))
motokaro_NET[x * NET_SIZE + y] = UNBLOCKED;
if (motokaro_NET[x * NET_SIZE + y] == UNBLOCKED) {
if (x > 0 && motokaro_NET[(x - 1) * NET_SIZE + y] == OPENED) {
motokaro_NET[(x - 1) * NET_SIZE + y] = UNBLOCKED;
check_all(x - 1, y);
}
if (x < NET_SIZE - 1 && motokaro_NET[(x + 1) * NET_SIZE + y] == OPENED) {
motokaro_NET[(x + 1) * NET_SIZE + y] = UNBLOCKED;
check_all(x + 1, y);
}
if (y > 0 && motokaro_NET[x * NET_SIZE + y - 1] == OPENED) {
motokaro_NET[x * NET_SIZE + y - 1] = UNBLOCKED;
check_all(x, y - 1);
}
if (y < NET_SIZE - 1 && motokaro_NET[x * NET_SIZE + y + 1] == OPENED) {
motokaro_NET[x * NET_SIZE + y + 1] = UNBLOCKED;
check_all(x, y + 1);
}
}
return;
}
private static void pour_into(int x, int y) {
motokaro_NET[x * NET_SIZE + y] = OPENED;
if (x == 0) {
motokaro_NET[x * NET_SIZE + y] = UNBLOCKED;
}
check_all(x, y);
return;
}
private static boolean is_unblocked() {
for (int i = 0; i < NET_SIZE; i++)
if (motokaro_NET[(NET_SIZE - 1) * NET_SIZE + i] == UNBLOCKED)
return true;
return false;
}
public static double motokaro_p() {
int cnt = 1;
Arrays.fill(motokaro_NET, CLOSED);
Random r = new Random();
for (int i = 0; i < NET_SIZE * NET_SIZE; i++) {
if (!is_unblocked()) {
int x = r.nextInt(NET_SIZE);
int y = r.nextInt(NET_SIZE);
if (motokaro_NET[x * NET_SIZE + y] != CLOSED) {
i--;
continue;
}
cnt++;
pour_into(x, y);
//print_net();
} else break;
}
//System.out.println(cnt);
return (double) cnt / (NET_SIZE * NET_SIZE);
}
public static void main(String[] args) {
double result = 0;
System.out.println("Net_Size == "+NET_SIZE);
for (int i = 0; i < RUNTIMES; i++) {
double temp = motokaro_p();
result += temp;
//System.out.println("第"+i+"次的p值："+temp);
}
System.out.println((double) result / RUNTIMES);
}
}
版权所有©️spio 抄作业不可取
浙公网安备 33010602011771号