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



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

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

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

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

การแจกแจงแกมมา (γ分布, 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
-- opencv
-- pytorch
機械学習
-- ニューラル
     ネットワーク
maya
javascript
確率論
日本での日記
中国での日記
-- 北京での日記
-- 香港での日記
-- 澳門での日記
台灣での日記
北欧での日記
他の国での日記
qiita
その他の記事

記事の類別



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

  記事を検索

  おすすめの記事

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

月別記事

2021年

1月 2月 3月 4月
5月 6月 7月 8月
9月 10月 11月 12月

2020年

1月 2月 3月 4月
5月 6月 7月 8月
9月 10月 11月 12月

2019年

1月 2月 3月 4月
5月 6月 7月 8月
9月 10月 11月 12月

2018年

1月 2月 3月 4月
5月 6月 7月 8月
9月 10月 11月 12月

2017年

1月 2月 3月 4月
5月 6月 7月 8月
9月 10月 11月 12月

もっと前の記事

ไทย

日本語

中文