RejectionSampling-拒绝性采样 ripple rejection

我们需要随机生成一个随机变量 X 满足概率分布f 时,我们可以选择一个更容易生成的概率分布g 。当然,不是任意的 g 都可以胜任,g 必须满足:

这种方法就是拒绝性采样。

理论支持

采用一组相互独立并符合均匀分布的随机变量 (Un) ,和另一组相互独立的随机变量(Yn) 满足概率分布g 并且与 Un 相互独立. 那么,我们定义:

随机变量 n0满足几何分布,其参数是 (所以期望是 a)。 X 是一个满足概率分布 f 的随机变量.

需要特别强调的是,当我们选择的 a 越接近1, 模拟程序的运行就越快。 极限情况是当 a = 1 时,f = g.

拒绝性算法

  1. 模拟
  2. 接受或者拒绝

应用

我们这里介绍几个JAVA模拟的具体例子。我们把被接受的点的数目在所有点(被接受的 +被拒绝的)中所占的比率成为效率。所测试的点越多,效率就越接近 。 在以下的两个例子中,设 .

在单位圆上模拟均匀分布

我们是要模拟满足概率分布 f 的随机变量,其中f 满足 .
我们把 g 选择为均匀分布在一个中心在原点,边长为2 的正方形上。然后我们选择 。 我们可以得出:
RejectionSampling-拒绝性采样 ripple rejection


对标准正态分布的模拟

我们是要模拟满足概率分布 f 的随机变量,并有f 满足 .
我们可以很合理地假设方程 f 在区间 [ − 5,5] 外为零. 我们把 g 选择为一个定义域在区间 [ −5,5] 的均匀分布,并选择 。 那么,我们有:


编码

Cercle.java

  1. import java.awt.*;
  2. import java.applet.*;
  3. public class Cercle extends Applet {
  4.     static final long serialVersionUID = 1L;
  5.     int n = 10000;      // number of points
  6.     int count = 0;      // number of tested points
  7.     Point tab[];   // point array
  8.     // --------------------------------------------------
  9.     //  The applet is initialised
  10.     public void init() {
  11.    // resizing
  12.    this.setSize(500, 500);
  13.    // the array is initialised
  14.    tab = new Point[n];
  15.    int k = 0;
  16.    // the array is filled with n samplings
  17.    while (k < n) {
  18.        double x, y;
  19.        x = 2. * Math.random() - 1.;
  20.        y = 2. * Math.random() - 1.;
  21.        // acceptatance or rejection
  22.        if (x * x + y * y <= 1.) {
  23.            tab[k] = new Point(x, y);
  24.            k++;
  25.        }
  26.        count++;
  27.    }
  28.     } // init()
  29.     // --------------------------------------------------
  30.     public void start() {
  31.     } // start()
  32.     //  --------------------------------------------------
  33.     //  The applet is displayed
  34.     public void paint(Graphics g) {
  35.    int k;
  36.    // displaying the circle
  37.    g.drawOval(0, 0, 500, 500);
  38.    // displaying the points
  39.    for (k = 0; k < n; k++) {
  40.        int x, y;
  41.        x = (int) (250. * (1. + tab[k].x));
  42.        y = (int) (250. * (1. - tab[k].y));
  43.        g.drawLine(x, y, x, y);
  44.    }
  45.    // displaying the algorithm's efficiency (#accepted divided by #tested)
  46.    g.setColor(Color.BLACK);
  47.    g.drawString("eff = " + ((double) n) / count, 30, 30);
  48.     } // paint()
  49.     // --------------------------------------------------
  50. } // class Cercle

Gauss.java

  1. import java.awt.*;
  2. import java.applet.*;
  3. public class Gauss extends Applet {
  4.     static final long serialVersionUID = 1L;
  5.     int n = 5000;       // number of points
  6.     int count = 0;      // number of tested points
  7.     int histo[];   // histogram of the distribution
  8.     int m;               // maximum of the distribution values
  9.     // --------------------------------------------------
  10.     // Gaussian distribution function we want to simulate 
  11.     double f(double t) {
  12.    return Math.exp(-t * t / 2.) / Math.sqrt(2. * Math.PI);
  13.     } // f()
  14.     // --------------------------------------------------
  15.     // uuniform probability distribution across the interval [-5,5] (easily simulated)
  16.     double g(double t) {
  17.    return 1. / 10;
  18.     } // g()
  19.     // --------------------------------------------------
  20.     // The applet is initialised
  21.     public void init() {
  22.    // resizing
  23.    this.setSize(500, 500);
  24.    // the histogram is initialised
  25.    histo = new int[100];
  26.    int k;
  27.    for (k = 0; k < 100; k++) {
  28.        histo[k] = 0;
  29.    }
  30.    // the histogram is filled with n samplings
  31.    k = 0;
  32.    while (k < n) {
  33.        double a = 10. / Math.sqrt(2. * Math.PI);
  34.        double x, u;
  35.        // x is sampled from the g-distribution
  36.        x = 10. * (Math.random() - 0.5);
  37.        // a uniform variable is simulated
  38.        u = Math.random();
  39.        // acceptance or rejection
  40.        if (u < f(x) / (a * g(x))) {
  41.            int i;
  42.            i = (int) (10. * (x + 5.));
  43.            histo[i]++;
  44.            k++;
  45.        }
  46.        count++;
  47.    }
  48.    // the maximum of the histogram values is computed
  49.    m = 0;
  50.    for (k = 0; k < 100; k++) {
  51.        m = Math.max(m, histo[k]);
  52.    }
  53.     } // init()
  54.     // --------------------------------------------------
  55.     public void start() {
  56.     } // start()
  57.     // --------------------------------------------------
  58.     // The applet is displayed
  59.     public void paint(Graphics g) {
  60.    int i;
  61.    // displaying the histogram
  62.    g.setColor(Color.RED);
  63.    for (i = 0; i < 100; i++) {
  64.        // the histogram values are converted for the viewport (size = 500*500)
  65.        int y;
  66.        y = (int) (500. * (1. - ((double) histo[i]) / m));
  67.        g.drawRect(5 * i, y, 5, 500 - y);
  68.    }
  69.    // displaying the Gaussian curve
  70.    g.setColor(Color.BLUE);
  71.    i = 0;
  72.    for (i = 0; i < 500; i++) {
  73.        double b = 500. * Math.sqrt(2 * Math.PI);
  74.        // the curve values are converted for the viewport (size = 500*500)
  75.        double x1, x2;
  76.        x1 = (i - 250.) / 50;
  77.        x2 = (i + 1 - 250.) / 50;
  78.        int y1, y2;
  79.        y1 = 500 - (int) (b * f(x1));
  80.        y2 = 500 - (int) (b * f(x2));
  81.        g.drawLine(i, y1, i + 1, y2);
  82.    }
  83.    // displaying the algorithm's efficiency (#acceptated divided by #tested)
  84.    g.setColor(Color.BLACK);
  85.    g.drawString("eff = " + ((double) n) / count, 30, 30);
  86.     } // paint()
  87.     // --------------------------------------------------
  88. } // class Gauss

  

爱华网本文地址 » http://www.413yy.cn/a/25101011/102580.html

更多阅读

窦性心律不齐不是病? 呼吸性窦性心律不齐

窦性心律不齐不是病?——简介由于心脏有精致而完整的传导系统,在正常情况下,窦房结发出的搏动频率主宰整个心脏的跳动,为心脏跳动的主导节律,称之为窦性心律。而实际上窘窦性心律不齐本身并不是一种病,下面就来详细解析下!窦性心律不齐不

如何预防月经性哮喘的发作 哮喘发作

?  有些女性哮喘患者,每到月经前或月经期就会哮喘发作,真叫人感到莫名其妙,医学上把这些哮喘称为“月经性哮喘”。月经性哮喘可发生于月经前,也可发生于月经期。那么,如何预防月经性哮喘的发作?  如何预防月经性哮喘的发作?这类哮

怎么样延长性时间

怎么样延长性时间——简介怎么样延长性时间成了大家最关注的话题,有些人认为肾虚早*泄是一件不光彩的事情,选择了拖延给治疗带来了很大的困难。怎么样延长性时间——方法/步骤怎么样延长性时间 1、体质性因素

热性水果有哪些? 哪些水果是热性的

热性水果有哪些?——简介一般情况,身体不舒服的人,会有一些饮食的禁忌,例如,不能吃寒性的水果,或者是不能吃热性的水果。但,好些人不知道哪些水果是寒性的,哪些水果是热性的,今天,我们一起来看一下,热性的水果有哪些。热性水果有哪些?——方法/

声明:《RejectionSampling-拒绝性采样 ripple rejection》为网友樱花分享!如侵犯到您的合法权益请联系我们删除