The Nature of Code閱讀心得與Python實作：0.5 A Custom Distribution...

這一節的標題是
A Custom Distribution of Random Numbers
因為方格子標題字數限制，所以沒完整顯現

先前的隨機漫步寫法，會有oversampling的問題，也就是同一個地方可能會走過很多次。老是在相同的地方走來走去，那可不是個好習慣，因為這樣會浪費不少時間和精力。那怎樣才可以避免這樣的問題呢？時不時的跨大步是個解決的辦法。原本的走法，步伐大小是固定的，利用時不時跨大步的方式，就可以減少oversampling的數量。這種時不時跨大步的走法，稱為Lévy flight。

根據維基百科上的說明，Lévy flight是Benoît Mandelbrot所創造的名詞，用來定義呈現特定分布方式的步伐大小。Mandelbrot這位大師似乎很喜歡創造名詞，fractal這個字，就是他的傑作。

在0.3節中有提到幾個方法，能夠用來產生符合特定分布的亂數。利用那些方法，就能讓walker具有Lévy flight的能力。例如

r = random.random()
if r < 0.01:
    x_step = random.uniform(-100, 100)
    y_step = random.uniform(-100, 100)
else:
    x_step = random.uniform(-1, 1)
    y_step = random.uniform(-1, 1)

原書針對這個寫法指出，這樣子寫，採取大步伐的機率會是1%。但這是不很精確的說法，因為從random.uniform(-100, 100)取得的亂數，也有可能會落在-1～1之間，而這是屬於小步伐的範圍，這會導致採取大步伐的機率，實際上會比1%還小一些。所以，比較精確的說法，應該是walker有1%的機率從比較大範圍的距離，也就是-100～100中，挑選出要用的步伐大小。不過不管怎麼說，這個機率是固定的，不管是100或99，被選中的機率都是一樣的。

雖然上述的寫法能夠讓walker具有Lévy flight的能力，但如果我們想要的，是像「越大的步伐，被使用的機率越小；越小的步伐，被使用的機率越高」這樣更複雜的能力時，又該怎麼辦呢？這時候，我們會需要用到具有特定分布方式的亂數，而這種亂數，是沒辦法使用Python提供的亂數產生器，如random.random()、random.uniform()、random.gauss()等，直接產生的。在遇到這種需要自己製作符合特定分布方式的亂數時，可以考慮使用accept-reject algorithm。

accept-reject algorithm也叫做rejection sampling或acceptance-rejection method，是蒙地卡羅方法(Monte Carlo method)的一種，可以用來產生符合特定機率分布的亂數。假設只考慮0～1間的亂數，其做法如下：

步驟1：取uniform distribution亂數r1。

步驟2：依照所需的機率分布，計算r1的機率p。

步驟3：取另一個uniform distribution亂數r2。

步驟4：若r2 < p，則r1符合資格，是我們所要的亂數。結束程式。

步驟5：若r2 ≥ p，則r1不符合資格，丟棄。回步驟1。

根據這個演算法，如果我們設定p = 1 - r1，那當r1 = 0.1時，就代表r1有90%的機率會符合資格而被取用；如果r1 = 0.83，則r1符合資格被取用的機率就是17%。很明顯的，這樣子的設定，會讓越大的數字，被取用的機率越低。反之，如果設定p = r1，則會讓越大的數字，被取用的機率越高。

雖然accept-reject algorithm可以產生符合我們所要的分布方式的亂數，但它的效率並不高，因為有可能會丟棄一大堆的r1之後，還是沒能找到符合資格可以取用的亂數。在取r1時，改用其他分布方式的亂數，可以改善執行效率，但有些額外的條件需要考慮。詳細的原理與說明，可以參考這篇文章：〈What is Rejection Sampling?〉

下面這個例子，就是應用accept-reject algorithm，在設定p = r1的情況下，符合資格的數字的分布狀況。

Example 0.5: An Accept-Reject Distribution

從圖可以看出，數字越大，符合資格的機率就越大。程式如下：

def accept_reject():
    while True:
        r1 = random.random()
        probability = r1
        r2 = random.random()
        if (r2 < probability):
            return r1


# python version 3.10.9
import random
import sys

import pygame  # version 2.3.0


pygame.init()

pygame.display.set_caption("Example 0.5: An Accept-Reject Distribution")

BLACK = (0, 0, 0)
WHITE = (255, 255, 255)

screen_size = WIDTH, HEIGHT = 640, 360
screen = pygame.display.set_mode(screen_size)

FPS = 60
frame_rate = pygame.time.Clock()

counter = [0]*20

rect_width = WIDTH/len(counter)

while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
            sys.exit()
    
    screen.fill(WHITE)

    index = int(len(counter)*accept_reject())
    counter[index] += 1
    
    for i in range(len(counter)):
        rect = pygame.Rect(i*rect_width, HEIGHT-counter[i], rect_width, counter[i])
        pygame.draw.rect(screen, BLACK, rect, 2)
    
    pygame.display.update()
    frame_rate.tick(FPS)

Exercise 0.6

把accept_reject()裡頭的probability設定成r1*r1：

def accept_reject():
    while True:
        r1 = random.random()
        probability = r1*r1
        r2 = random.random()
        if (r2 < probability):
            return r1

另外，把Walker類別的step()修改成

def step(self):
    step_size = 10*accept_reject()        
    
    step_x = random.uniform(-step_size, step_size)
    step_y = random.uniform(-step_size, step_size)
    
    self.x += int(step_x)
    self.y += int(step_y)

因為accept_reject()產生的是介於0～1間的數字，所以必須乘上10來轉換成實際上的步伐大小。另外要注意的是，step_x和step_y都是浮點數，若使用set_at()來顯示walker的軌跡，因為只接受整數座標，所以必須用int()將其轉為整數。