爱华网这种方法就是拒绝性采样。
理论支持
采用一组相互独立并符合均匀分布的随机变量 (Un) ，和另一组相互独立的随机变量(Yn) 满足概率分布g 并且与 Un 相互独立. 那么，我们定义:
和
随机变量 n0满足几何分布，其参数是 (所以期望是 a)。 X 是一个满足概率分布 f 的随机变量.
需要特别强调的是，当我们选择的 a 越接近1， 模拟程序的运行就越快。 极限情况是当 a = 1 时，f = g.
拒绝性算法
- 模拟
- 接受或者拒绝
应用
我们这里介绍几个JAVA模拟的具体例子。我们把被接受的点的数目在所有点（被接受的 +被拒绝的）中所占的比率成为效率。所测试的点越多，效率就越接近 。 在以下的两个例子中，设 .
在单位圆上模拟均匀分布
我们是要模拟满足概率分布 f 的随机变量，其中f 满足 .
我们把 g 选择为均匀分布在一个中心在原点，边长为2 的正方形上。然后我们选择 。 我们可以得出:
对标准正态分布的模拟
我们是要模拟满足概率分布 f 的随机变量，并有f 满足 .
我们可以很合理地假设方程 f 在区间 [ − 5,5] 外为零. 我们把 g 选择为一个定义域在区间 [ −5,5] 的均匀分布，并选择 。 那么，我们有:
编码
Cercle.java
import java.awt.*;
import java.applet.*;
public class Cercle extends Applet {static final long serialVersionUID = 1L;
int n = 10000; // number of points
int count = 0; // number of tested points
Point tab[]; // point array
// --------------------------------------------------
// The applet is initialised
public void init() {// resizing
this.setSize(500, 500);
// the array is initialised
tab = new Point[n];
int k = 0;
// the array is filled with n samplings
while (k < n) {double x, y;
x = 2. * Math.random() - 1.;
y = 2. * Math.random() - 1.;
// acceptatance or rejection
if (x * x + y * y <= 1.) {tab[k] = new Point(x, y);
k++;
}
count++;
}
} // init()
// --------------------------------------------------
public void start() {} // start()
// --------------------------------------------------
// The applet is displayed
public void paint(Graphics g) {int k;
// displaying the circle
g.drawOval(0, 0, 500, 500);
// displaying the points
for (k = 0; k < n; k++) {int x, y;
x = (int) (250. * (1. + tab[k].x));
y = (int) (250. * (1. - tab[k].y));
g.drawLine(x, y, x, y);
}
// displaying the algorithm's efficiency (#accepted divided by #tested)
g.setColor(Color.BLACK);
g.drawString("eff = " + ((double) n) / count, 30, 30);} // paint()
// --------------------------------------------------
} // class Cercle
Gauss.java
import java.awt.*;
import java.applet.*;
public class Gauss extends Applet {static final long serialVersionUID = 1L;
int n = 5000; // number of points
int count = 0; // number of tested points
int histo[]; // histogram of the distribution
int m; // maximum of the distribution values
// --------------------------------------------------
// Gaussian distribution function we want to simulate
double f(double t) {return Math.exp(-t * t / 2.) / Math.sqrt(2. * Math.PI);
} // f()
// --------------------------------------------------
// uuniform probability distribution across the interval [-5,5] (easily simulated)
double g(double t) {return 1. / 10;
} // g()
// --------------------------------------------------
// The applet is initialised
public void init() {// resizing
this.setSize(500, 500);
// the histogram is initialised
histo = new int[100];
int k;
for (k = 0; k < 100; k++) {histo[k] = 0;
}
// the histogram is filled with n samplings
k = 0;
while (k < n) {double a = 10. / Math.sqrt(2. * Math.PI);
double x, u;
// x is sampled from the g-distribution
x = 10. * (Math.random() - 0.5);
// a uniform variable is simulated
u = Math.random();
// acceptance or rejection
if (u < f(x) / (a * g(x))) {int i;
i = (int) (10. * (x + 5.));
histo[i]++;
k++;
}
count++;
}
// the maximum of the histogram values is computed
m = 0;
for (k = 0; k < 100; k++) {m = Math.max(m, histo[k]);
}
} // init()
// --------------------------------------------------
public void start() {} // start()
// --------------------------------------------------
// The applet is displayed
public void paint(Graphics g) {int i;
// displaying the histogram
g.setColor(Color.RED);
for (i = 0; i < 100; i++) {// the histogram values are converted for the viewport (size = 500*500)
int y;
y = (int) (500. * (1. - ((double) histo[i]) / m));
g.drawRect(5 * i, y, 5, 500 - y);
}
// displaying the Gaussian curve
g.setColor(Color.BLUE);
i = 0;
for (i = 0; i < 500; i++) {double b = 500. * Math.sqrt(2 * Math.PI);
// the curve values are converted for the viewport (size = 500*500)
double x1, x2;
x1 = (i - 250.) / 50;
x2 = (i + 1 - 250.) / 50;
int y1, y2;
y1 = 500 - (int) (b * f(x1));
y2 = 500 - (int) (b * f(x2));
g.drawLine(i, y1, i + 1, y2);
}
// displaying the algorithm's efficiency (#acceptated divided by #tested)
g.setColor(Color.BLACK);
g.drawString("eff = " + ((double) n) / count, 30, 30);} // paint()
// --------------------------------------------------
} // class Gauss
