φυβλαςのβλογ
phyblasのブログ



ความน่าจะเป็นเบื้องต้นสำหรับเขียนโปรแกรม บทที่ ๑๖: การแจกแจงแกมมากับพารามิเตอร์ของการแจกแจงความน่าจะเป็นแบบต่างๆ
เขียนเมื่อ 2020/09/07 19:08
แก้ไขล่าสุด 2021/09/28 16:42

ต่อจาก บทที่ ๑๕

ในบทนี้จะเป็นเรื่องของการแจกแจงแกมมา ซึ่งใช้เป็นการแจกแจงความน่าจะเป็นของค่าพารามิเตอร์ของการแจกแจงอื่นๆหลายชนิด

การแจกแจงแกมมา

การแจกแจงแกมมา (γ分布, gamma distribution) เป็นการแจกแจงที่อาจใช้เป็นความน่าจะเป็นก่อนหน้าสังยุคของการแจกแจงแบบปกติ, การแจกแจงปัวซง และ การแจกแจงเอ็กซ์โพเนนเชียล

ก่อนอื่นขอเริ่มจากแนะนำรูปทั่วไปของการแจกแจงแกมมา

ฟังก์ชันการแจกแจงความหนาแน่นความน่าจะเป็นของการแจกแจงแกมมาเป็นดังนี้


พารามิเตอร์มีอยู่ 2 ตัวคือ ν กับ β ซึ่งค่าอาจเป็นจำนวนจริงบวกใดๆ

ในบางครั้งอาจเขียนในรูปแบบนี้แทน


โดย θ=1/β เป็นส่วนกลับของ β

จะเขียนแบบไหนแล้วแต่ว่าจะใช้ในสมการไหน สำหรับกรณีที่กำลังจะใช้ต่อไปนี้จะเขียนในรูป β เนื่องจากสะดวกในการคำนวณมากกว่า

นอกจากนี้อักษรที่ใช้เป็นตัวแปรก็ต่างกันไปมากขึ้นอยู่กับว่าตำราไหน ที่จริงตัวแปรที่เขียนเป็น ν ในที่นี้มักถูกเขียนด้วย α หรือ λ มากกว่า แต่ทั้งสองตัวนี้ถูกใช้เป็นตัวแปรอย่างอื่นไปแล้วจึงขอใช้ ν (นิว)

ส่วน Γ คือฟังก์ชันแกมมา ซึ่งในกรณีทั่วไปคำนวณได้จาก


ถ้า ν มีค่าเป็นจำนวนเต็ม ก็จะเท่ากับแฟกทอเรียลลบหนึ่ง ซึ่งสามารถคำนวณได้ง่ายเป็น


หรือกรณีที่เป็นจำนวนเต็มบวกด้วย 1/2 จะคำนวณได้เป็น


กรณีที่ใช้กับพารามิเตอร์ของการแจกแจงเอ็กซ์โพเนนเชียลหรือการแจกแจงปัวซง ปกติแล้ว ν มักจะเป็นจำนวนเต็ม และเมื่อใช้กับการแจกแจงแบบปกติมักเป็นจำนวนเต็มหรือมีเศษ 1/2

กราฟแสดงตัวอย่างค่าความหนาแน่นของความน่าจะเป็นของการแจกแจงแกมมาในกรณีต่างๆ



ค่าคาดหมายของการแจกแจงแกมมาจะอยู่ที่


ค่าสูงสุดของการแจกแจงแกมมาจะอยู่ที่


เพียงแต่ถ้า ν<1 จะไม่มีจุดสูงสุด

ความแปรปรวนเท่ากับ


ในที่นี้จะใช้ แทนการแจกแจงแกมมาโดยเขียนแบบนี้





การแจกแจงความน่าจะเป็นของพารามิเตอร์ความเที่ยงตรงในการแจกแจงแบบปกติ

ต่อไปเริ่มยกตัวอย่างการใช้การแจกแจงแกมมา โดยเริ่มจากใช้เป็นความน่าจะเป็นก่อนหน้าสังยุคของการแจกแจงแบบปกติ

รูปทั่วไปของการแจกแจงแบบปกติคือ


ในบทที่ ๑๕ ได้เขียนถึงการแจกแจงพารามิเตอร์ μ โดยถือว่า σ เป็นค่าที่รู้อยู่แน่นอนแล้ว ส่วนในบทนี้จะทำในกรณีตรงกันข้ามกัน คือมีค่า μ อยู่ตายตัวแล้วต้องการหา σ

เพียงแต่เพื่อความสะดวก แทนที่จะพิจารณาการแจกแจงของส่วนเบี่ยงเบนมาตรฐาน σ หรือความแปรปรวน σ2 โดยทั่วไปมักจะคำนวณในรูปของส่วนกลับของความแปรปรวน เรียกว่า พารามิเตอร์ความเที่ยงตรง (精度参数, precision parameter)


ที่จริงถ้าจะหาการแจกแจงในรูปของความแปรปรวน σ2 ก็ทำได้ แต่จะกลายเป็นการแจกแจงแกมมาผกผัน (逆γ分布, inverse gamma distribution) ซึ่งก็คล้ายการแจกแจงแกมมาแต่ไม่นิยมใช้มากเท่า ในที่นี้ก็จะใช้ α เพื่อจะได้เป็นการแจกแจงแกมมาธรรมดา

เมื่อแทนด้วย α การแจกแจงแบบปกติจะเขียนใหม่ได้เป็นดังนี้


ปกติเมื่อเขียนแทนด้วย ค่าตัวขวาสุดจะเขียนเป็นความแปรปรวน σ2 ในที่นี้แม้จะเขียนในรูปของ α ก็ให้เขียนเป็นส่วนกลับ α-1 ซึ่งเท่ากับ σ2

ในที่นี้จะพิจารณาค่าพารามิเตอร์ α โดยใช้ทฤษฎีบทของเบย์ เช่นเดียวกับที่ทำกับ μ ในบทที่แล้ว แต่ครั้งนี้ตัวแปรที่พิจารณาคือ α

เมื่อมีข้อมูล x1 ที่เกิดจากการสุ่มโดยมีการแจกแจงความน่าจะเป็นแบบปกติ สามารถหา α ตามทฤษฎีบทของเบย์ได้โดย


ความน่าจะเป็นก่อนหน้าเมื่อไม่มีข้อมูลอะไรอาจให้เป็นค่าคงตัวไป ซึ่งการแจกแจงแบบค่าคงตัวก็ถือเป็นฟังก์ชันแกมมารูปแบบหนึ่ง คือในกรณีที่ ν=1 และ β เป็นค่าเล็กๆเข้าใกล้ 0


ส่วนฟังก์ชันควรจะเป็นคือ


เมื่อแทนลงไปแล้วก็จะคำนวณการแจกแจงความน่าจะเป็นภายหลังได้เป็น


เมื่อพิจารณาค่า α จะพบว่า α สามารถเขียนในรูปการแจกแจงแกมมาเป็น


ก็จะได้ว่า


และถ้าหากมีข้อมูล x2 ที่เป็นผลจากการแจกเดิมเพิ่มเข้ามาอีกก็จะได้ว่า


และเมื่อมีข้อมูลเข้ามาทั้งหมด n ตัวก็จะได้ว่า


จะเห็นว่าเมื่อมีข้อมูลใส่เข้ามา ค่า α จะบวกเพิ่มไปตัวละ 1/2 ส่วน β จะเป็นครึ่งหนึ่งของผลรวมของระยะห่างกำลังสองจาก μ

ในกรณีที่ความน่าจะเป็นตั้งต้นมีค่าอยู่แล้วเป็นฟังก์ชันแกมมาที่มีพารามิเตอร์เริ่มต้นเป็น ν=ν0, β=β0


การแจกแจงความน่าจะเป็นภายหลังก็จะได้เป็น


ต่อไปเป็นโค้ดตัวอย่าง โดยจะสร้างข้อมูลที่สุ่มโดยแจกแจงแบบปกติ โดยที่รู้ค่า μ จริงอยู่แล้วแต่ไม่รู้ σ แล้วพยายามหาค่า σ จากการแจกแจงนี้ ดูลักษณะการแจกแจงและจุดสูงสุดที่ได้เมื่อเพิ่มข้อมูลเข้ามาเรื่อยๆ
import random

ν0 = 1 # ν เริ่มต้น
β0 = 0.01 # β เริ่มต้น
μ = 8 # μ จริงๆของการแจกแจง
σ = 2 # σ จริงๆของการแจกแจง
α = 1/σ**2 # แปลงจาก σ เป็นค่าความเที่ยงตรง
x = [] # ลิสต์เก็บค่าที่สุ่มได้ทั้งหมด
for n in range(21):
    # คำนวณค่า ν และ β ใหม่ในแต่ละรอบ
    ν = ν0+n/2
    β = β0+sum((xi-μ)**2/2 for xi in x)
    α_max = (ν-1)/β # หาจุดสูงสุด
    print('n=%d, ν=%.2f, β=%.2f, α_max=%.2f'%(n,ν,β,α_max))
    # สุ่มค่าใหม่ในแต่ละรอบ
    x += [random.gauss(μ,σ)]

ผลที่ได้
n=0, ν=1.00, β=0.01, α_max=0.00
n=1, ν=1.50, β=0.55, α_max=0.91
n=2, ν=2.00, β=0.55, α_max=1.80
n=3, ν=2.50, β=0.56, α_max=2.70
n=4, ν=3.00, β=1.32, α_max=1.51
n=5, ν=3.50, β=1.36, α_max=1.84
n=6, ν=4.00, β=1.45, α_max=2.06
n=7, ν=4.50, β=9.46, α_max=0.37
n=8, ν=5.00, β=11.80, α_max=0.34
n=9, ν=5.50, β=12.79, α_max=0.35
n=10, ν=6.00, β=13.37, α_max=0.37
n=11, ν=6.50, β=17.01, α_max=0.32
n=12, ν=7.00, β=17.28, α_max=0.35
n=13, ν=7.50, β=17.91, α_max=0.36
n=14, ν=8.00, β=17.93, α_max=0.39
n=15, ν=8.50, β=18.33, α_max=0.41
n=16, ν=9.00, β=24.77, α_max=0.32
n=17, ν=9.50, β=25.38, α_max=0.33
n=18, ν=10.00, β=25.44, α_max=0.35
n=19, ν=10.50, β=35.43, α_max=0.27
n=20, ν=11.00, β=39.57, α_max=0.25

เมื่อนำมาแสดงเป็นภาพเคลื่อนไหวจะได้แบบนี้



ภาพด้านบนแสดงการแจกแจงของ α จะเห็นว่ายิ่งเพิ่มข้อมูล ตำแหน่ง α ที่ค่าสูงสุดก็ยิ่งลู่เข้าตำแหน่งเส้นประสีชมพู คือ α=0.25 ซึ่งเป็นค่าจริง และกราฟการแจกแจงก็ค่อยๆบีบแคบลงมา

ส่วนภาพด้านล่างแสดงจุดตำแหน่งของข้อมูลที่สุ่มได้ และเส้นกราฟแสดงการแจกแจงของค่า x โดยเส้นประสีฟ้าคือค่าจริง ส่วนเส้นเขียนคือค่าที่ใช้ α สูงสุดจากการแจกแจงในแต่ละรอบซึ่งเปลี่ยนไปเรื่อยๆเมื่อได้ข้อมูลเข้ามา ยิ่งผ่านไปกราฟทั้งสองก็ยิ่งใกล้เคียงกัน นั่นคือเราสามารถหาค่า α ได้ใกล้เคียงค่าจริง

สำหรับการแจกแจงแบบปกติหลายตัวแปรจะทำคล้ายกันแต่จะใช้การแจกแจงวิชาร์ตเป็นความน่าจะเป็นก่อนหน้าสังยคแทน ซึ่งเขียนถึงในบทที่ ๑๗


การแจกแจงความน่าจะเป็นของพารามิเตอร์ของการแจกแจงเอ็กซ์โพเนนเชียล

ต่อไปมาดูกรณีของการแจกแจงเอ็กซ์โพเนนเชียล ซึ่งใช้การแจกแจงแกมมาเป็นความน่าจะเป็นก่อนหน้าสังยุคเช่นกัน

ในที่นี้ขอเขียนแทนการแจกแจงเอ็กซ์โพเนนเชียลด้วย แบบนี้


ค่าพารามิเตอร์ที่พิจารณาคราวนี้คือ λ

เนื่องจากขั้นตอนต่างๆก็จะคล้ายกับที่ทำไปในกรณีการแจกแจงแบบปกติ ต่อไปนี้จะอธิบายแค่คร่าวๆอย่างรวบรัด

รายละเอียดเกี่ยวกับการแจกแจงเอ็กซ์โพเนนเชียลอ่านได้ในบทที่ ๑๐

เมื่อมีค่า x1 ซึ่งได้จากการสุ่มด้วยการแจกแจงแบบเอ็กซ์โพเนนเชียล ฟังก์ชันควรจะเป็นคือ


ให้การแจกแจงความน่าจะเป็นก่อนหน้าเป็นค่าคงตัว


แล้วก็จะคำนวณความน่าจะเป็นภายหลังได้ในรูปของการแจกแจงแกมมาตามนี้


และถ้ามีข้อมูลจากการสุ่มอยู่ n ตัว ก็จะได้ผลตามนี้


โดย ν จะเพิ่มขึ้นทีละ 1 ตามจำนวนข้อมูล และ β จะเป็นผลรวมของค่า

หากเตรียมความน่าจะเป็นก่อนหน้าเป็นฟังก์ชันแกมมาซึ่งมีค่าเริ่มต้น νเป็นν0 และ β เป็น 0


ก็จะได้การแจกแจงความน่าจะเป็นภายหลังของ λ เป็น




ตัวอย่าง ลองทำการสุ่มจริงแล้วหาการแจกแจงความน่าจะเป็นของพารามิเตอร์ λ จากข้อมูลที่ได้ซึ่งเพิ่มขึ้นเรื่อยๆในแต่ละขั้น
ν0 = 1 # ν เริ่มต้น
β0 = 0.01 # β เริ่มต้น
λ = 20 # λ จริงๆของการแจกแจง
# สุ่มค่า t
t = []
for i in range(100*λ):
    t += [random.uniform(0,100)]
t = sorted(t) # จัดเรียงใหม่
n = 0
x = [] # ลิสต์เก็บค่าระยะห่างระหว่าง t
for i in range(13):
    # คำนวณค่า ν และ β ใหม่ในแต่ละรอบ
    ν = ν0+n
    β = β0+sum(x)
    λ_max = (ν-1)/β # λ ที่ค่าสูงที่สุด
    print('n=%d, ν=%d, β=%.2f, λ_max=%.2f'%(n,ν,β,λ_max))
    
    # เอาค่า t ที่สุ่มได้มาเทียบระยะห่างทีละคู่
    for j in range(24):
        x += [t[n+1]-t[n]]
        n += 1

ผลที่ได้
n=0, ν=1, β=0.01, λ_max=0.00
n=24, ν=25, β=1.17, λ_max=20.54
n=48, ν=49, β=2.68, λ_max=17.91
n=72, ν=73, β=4.38, λ_max=16.44
n=96, ν=97, β=5.30, λ_max=18.10
n=120, ν=121, β=6.38, λ_max=18.80
n=144, ν=145, β=8.29, λ_max=17.37
n=168, ν=169, β=9.57, λ_max=17.55
n=192, ν=193, β=10.44, λ_max=18.40
n=216, ν=217, β=11.75, λ_max=18.38
n=240, ν=241, β=12.70, λ_max=18.89
n=264, ν=265, β=13.73, λ_max=19.23
n=288, ν=289, β=14.70, λ_max=19.59

ทำเป็นภาพเคลื่อนไหวได้แบบนี้



เส้นสีฟ้าคือจุดสูงสุดที่หาได้ขณะนั้น ส่วนเส้นประสีเขียวคือที่ λ=20 ซึ่งเป็นค่าจริง

ยิ่งใช้ผลการสุ่มหลายค่ามาก การแจกแจงก็ยิ่งบีบแคลลงโดยไปกองรวมที่ λ=20




การแจกแจงความน่าจะเป็นของพารามิเตอร์ของการแจกแจงปัวซง

สุดท้ายมาดูกรณีของการแจกแจงปัวซง ซึ่งใช้การแจกแจงแกมมาเป็นความน่าจะเป็นก่อนหน้าสังยุคเช่นกัน

รายละเอียดเกี่ยวกับการแจกแจงปัวซงอ่านได้ในบทที่ ๗

ในที่นี้จะเขียนแทนการแจกแจงปัวซงด้วย ดังนี้


พารามิเตอร์ของการแจกแจงปัวซงก็คือค่า λ

ให้ข้อมูล k1 ซึ่งมีการแจกแจงแบบปัวซง ฟังก์ชันควรจะเป็นคือ


เมื่อให้การแจกแจงความน่าจะเป็นก่อนหน้าเป็นค่าคงตัว


ก็จะได้การแจกแจงภายหลังดังนี้


เมื่อมีข้อมูล n ตัวก็จะได้เป็น


จะเห็นว่าตรงกันข้ามกับกรณีของการแจกแจงเอ็กซ์โพเนนเชียล โดย ν จะเพิ่มขึ้นตามผลรวมของค่าของข้อมูล ส่วน β จะเพิ่มตามจำนวนข้อมูล

กรณีที่ให้ความน่าจะเป็นก่อนหน้าเป็นฟังก์ชันแกมมาที่มีพารามิเตอร์เริ่มต้นเป็น ν0 และ β0


ก็จะได้การแจกแจงความน่าจะเป็นภายหลังดังนี้


ตัวอย่าง ทำการสุ่มข้อมูลให้มีการแจกแจงแบบปัวซง ค่อยๆเพิ่มข้อมูลเข้ามาทีละนิดแล้วดูความเปลี่ยนแปลงของการแจกแจงความน่าจะเป็นของค่า λ ที่ทำนายได้
ν0 = 1 # ค่า ν เริ่มต้น
β0 = 0.5 # ค่า β เริ่มต้น
λ = 4 # ค่า λ จริงๆของการแจกแจง
k = [0]*300
for i in range(300*λ):
    k[random.randint(0,299)] += 1

for i in range(11):
    n = i*30 # จำนวนข้อมูลที่จะใช้ในรอบนั้นๆ
    ν = ν0+sum(k[:n])
    β = β0+n
    λ_max = (ν-1)/β
    print('n=%d, ν=%d, β=%.1f, λ_max=%.3f'%(n,ν,β,λ_max))

ผลที่ได้
n=0, ν=1, β=0.5, λ_max=0.000
n=30, ν=125, β=30.5, λ_max=4.066
n=60, ν=233, β=60.5, λ_max=3.835
n=90, ν=346, β=90.5, λ_max=3.812
n=120, ν=468, β=120.5, λ_max=3.876
n=150, ν=610, β=150.5, λ_max=4.047
n=180, ν=737, β=180.5, λ_max=4.078
n=210, ν=850, β=210.5, λ_max=4.033
n=240, ν=969, β=240.5, λ_max=4.025
n=270, ν=1090, β=270.5, λ_max=4.026
n=300, ν=1201, β=300.5, λ_max=3.993

นำมาวาดเป็นภาพเคลื่อนไหวแสดงแต่ละขั้นตอนได้ดังนี้



ในที่นี้ภาพด้านบนคือการแจกแจงของค่า λ ส่วนภาพด้านล่างคือจำนวนรวมของค่า k แต่ละค่า ซึ่งรวมกันทั้งหมดเท่ากับ n แต่ละขั้นให้เพิ่ม n ทีละ 30

เมื่อข้อมูลเพิ่มมากขึ้นก็จะหาจุดสูงสุดได้ใกล้เคียงค่าจริงคือที่ λ=4 และการแจกแจงก็จะบีบแคบลง



บทถัดไป >> บทที่ ๑๗



-----------------------------------------

囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧

ดูสถิติของหน้านี้

หมวดหมู่

-- คณิตศาสตร์ >> ความน่าจะเป็น
-- คอมพิวเตอร์ >> เขียนโปรแกรม >> python

ไม่อนุญาตให้นำเนื้อหาของบทความไปลงที่อื่นโดยไม่ได้ขออนุญาตโดยเด็ดขาด หากต้องการนำบางส่วนไปลงสามารถทำได้โดยต้องไม่ใช่การก๊อปแปะแต่ให้เปลี่ยนคำพูดเป็นของตัวเอง หรือไม่ก็เขียนในลักษณะการยกข้อความอ้างอิง และไม่ว่ากรณีไหนก็ตาม ต้องให้เครดิตพร้อมใส่ลิงก์ของทุกบทความที่มีการใช้เนื้อหาเสมอ

目次

日本による名言集
モジュール
-- numpy
-- matplotlib

-- pandas
-- manim
-- opencv
-- pyqt
-- pytorch
機械学習
-- ニューラル
     ネットワーク
maya
javascript
確率論
日本での日記
中国での日記
-- 北京での日記
-- 香港での日記
-- 澳門での日記
台灣での日記
北欧での日記
他の国での日記
qiita
その他の記事

記事の類別



ติดตามอัปเดตของบล็อกได้ที่แฟนเพจ

  記事を検索

  おすすめの記事

การสร้างแบบจำลองสามมิติเป็นไฟล์ .obj วิธีการอย่างง่ายที่ไม่ว่าใครก็ลองทำได้ทันที
รวมรายชื่อนักร้องเพลงกวางตุ้ง
ภาษาจีนแบ่งเป็นสำเนียงอะไรบ้าง มีความแตกต่างกันมากแค่ไหน
ทำความเข้าใจระบอบประชาธิปไตยจากประวัติศาสตร์ความเป็นมา
เรียนรู้วิธีการใช้ regular expression (regex)
หลักการเขียนทับศัพท์ภาษาจีนกวางตุ้ง
การใช้ unix shell เบื้องต้น ใน linux และ mac
หลักการเขียนทับศัพท์ภาษาจีนกลาง
g ในภาษาญี่ปุ่นออกเสียง "ก" หรือ "ง" กันแน่
ทำความรู้จักกับปัญญาประดิษฐ์และการเรียนรู้ของเครื่อง
ค้นพบระบบดาวเคราะห์ ๘ ดวง เบื้องหลังความสำเร็จคือปัญญาประดิษฐ์ (AI)
หอดูดาวโบราณปักกิ่ง ตอนที่ ๑: แท่นสังเกตการณ์และสวนดอกไม้
พิพิธภัณฑ์สถาปัตยกรรมโบราณปักกิ่ง
เที่ยวเมืองตานตง ล่องเรือในน่านน้ำเกาหลีเหนือ
บันทึกการเที่ยวสวีเดน 1-12 พ.ค. 2014
แนะนำองค์การวิจัยและพัฒนาการสำรวจอวกาศญี่ปุ่น (JAXA)
เล่าประสบการณ์ค่ายอบรมวิชาการทางดาราศาสตร์โดยโซวเคนได 10 - 16 พ.ย. 2013
ตระเวนเที่ยวตามรอยฉากของอนิเมะในญี่ปุ่น
เที่ยวชมหอดูดาวที่ฐานสังเกตการณ์ซิงหลง
บันทึกการเที่ยวญี่ปุ่นครั้งแรกในชีวิต - ทุกอย่างเริ่มต้นที่สนามบินนานาชาติคันไซ
หลักการเขียนทับศัพท์ภาษาญี่ปุ่น
ทำไมจึงไม่ควรเขียนวรรณยุกต์เวลาทับศัพท์ภาษาต่างประเทศ
ทำไมถึงอยากมาเรียนต่อนอก
เหตุผลอะไรที่ต้องใช้ภาษาวิบัติ?

ไทย

日本語

中文