φυβλαςのβλογ
บล็อกของ phyblas



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

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

ในบทนี้เป็นเรื่องของการแจกแจงวิชาร์ต (威沙特分布, Wishart distribution) ซึ่งเป็นการแจกแจงความน่าจะเป็นของเมทริกซ์ความเที่ยงตรงของการแจกแจงแบบปกติหลายตัวแปร




เมทริกซ์ความเที่ยงตรง

ดังที่เคยเขียนถึงไปในบทที่ ๑๓ และบทที่ ๑๕ สมการการแจกแจงแบบปกติมักเขียนในลักษณะนี้


โดย Σ คือเมทริกซ์ความแปรปรวนร่วมเกี่ยว

แต่ในการคำนวณบางกรณีจะเป็นการสะดวกกว่าถ้าเขียนในรูปของเมทริกซ์ผกผันของ Σ


Λ นี้เรียกว่าเป็นเมทริกซ์ความเที่ยงตรง (精度矩阵, precision matrix) ซึ่งเป็นส่วนกลับของเมทริกซ์ความแปรปรวนร่วมเกี่ยว กล่าวคือ Σ ถ้ายิ่งมีตัวเลขมากยิ่งหมายถึงว่าค่ามีการแจกแจงมาก แต่ถ้า Λ มากจะหมายถึงค่ามีการแจกแจงน้อย ค่ามีความเที่ยงตรงมาก

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


เวลาเขียนแทนด้วย แบบนี้ปกติค่าตัวที่สองที่ใส่จะต้องเป็นเมทริกซ์ความแปรปรวนร่วมเกี่ยว ดังนั้นในกรณีที่ใช้เป็นเมทริกซ์ความเที่ยงตรงมักจะเขียนในรูปผกผัน Λ-1 แบบนี้ ไม่ใช่ Λ เฉยๆ




การแจกแจงวิชาร์ต

การแจกแจงวิชาร์ตมีรูปทั่วไปเป็นดังนี้


ที่เป็นตัวหนา ได้แก่ V และ Λ เป็นเมทริกซ์ขนาด m×m (โดย m เป็นจำนวนตัวแปร)

ถ้าเมทริกซ์ล้อมด้วย | | เช่น |V| และ |Λ| จะหมายถึงดีเทอร์มิแนนต์

ส่วน tr หมายถึงผลบวกค่าทั้งหมดในแนวทแยงของเมทริกซ์นั้น

พารามิเตอร์มี 2 ตัวคือ

- คือค่า องศาเสรี (自由度, degree of freedom) เทียบเท่าได้กับ ν/2 ในการแจกแจงแกมมา

- คือ เมทริกซ์สเกลปรับสเกล (尺度矩陣, scale matrix) มีขนาดเท่ากับเมทริกซ์ Λ

Γm เป็นฟังก์ชันแกมมาแบบหลายตัวแปร


จะเห็นว่าค่อนข้างยุ่ง แต่ก้อนนี้เป็นแค่ส่วนหนึ่งของค่าคงที่ซึ่งมีไว้หารเพื่อให้ผลรวมความน่าจะเป็นทั้งหมดเป็น 1 เท่านั้น ในที่นี้ไม่จำเป็นต้องสนใจก็ได้

ค่าต่างๆที่ไม่ได้เกี่ยวกับ Λ จะดึงมารวมไว้ใน แล้วเขียนใหม่ได้ดังนี้


ถ้า m=1 หมายถึงเป็นมิติเดียว V และ Λ ก็จะกลายเป็นเลขตัวเดียว แล้วหากเปลี่ยน เป็น และให้ค่าตัวเดียวใน V เป็น 1/β การแจกแจงนี้ก็จะกลายมาเป็นการแจกแจงแกมมา

ขอยกการแจกแจงแกมมามาเทียบให้ดูตรงนี้


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

ค่าคาดหมายของการแจกแจงของเมทริกซ์ความเที่ยงตรงของการแจกแจงวิชาร์ตมีค่าเป็น


ค่าของเมทริกซ์ความเที่ยงตรงที่มีความน่าจะเป็นในการแจกแจงสูงสุดคือ





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

ต่อไปจะพิจารณาการแจกแจงความน่าจะเป็นของเมทริกซ์ความเที่ยงตรง Λ เช่นเดียวกับที่ทำกับค่าพารามิเตอร์ความเที่ยงตรง α ในบทที่แล้ว

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

ความน่าจะเป็นก่อนหน้าเมื่อไม่มีข้อมูลก็สามารถให้เริ่มจากค่าคงที่ ซึ่งจะได้จากกรณีที่ และ Λ เป็นเมทริกซ์ที่ค่าข้างในมีค่าเล็กๆเข้าใกล้ 0


ให้ เป็นข้อมูล m มิติ (m ตัวแปร) ที่ได้จากการสุ่มแบบปกติหลายตัวแปร ฟังก์ชันควรจะเป็นก็จะได้เป็น


แล้วก็จะคำนวณความน่าจะเป็นภายหลังได้เป็นการแจกแจงวิชาร์ต


และถ้าหากมีข้อมูล n ตัวก็จะได้การแจกแจงวิชาร์ตดังนี้


กรณีที่มีการแจกแจงก่อนหน้าเป็น และ


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


ต่อไปเป็นตัวอย่างโค้ด ให้สร้างข้อมูลสองมิติที่ทำการสุ่มแบบปกติหลายตัวแปร แล้วลองหาการแจกแจงความน่าจะเป็น
import numpy as np

ν0 = 5
# เมทริกซ์ความแปรปรวนตั้งต้น
Σ0 = np.array([[0.5,0.1],
               [0.1,0.5]])

# เมทริกซ์ความแปรปรวนจริงๆที่ต้องการหา
Σ = np.array([[5.,1],
              [1,2.]])
μ = np.array([0,0])
x = np.random.multivariate_normal(μ,Σ,3600)
m = len(Σ) # จำนวนมิติ
for i in range(11):
    n = (6*i)**2 # จำนวนจุดข้อมูล
    # คำนวณค่า ν และเมทริกซ์ V ใหม่ในแต่ละรอบ
    ν = ν0+n
    V = np.linalg.inv(np.dot((x[:n]-μ).T,(x[:n]-μ))+Σ0)
    Σ_max = np.linalg.inv((ν-m-1)*V) # จุดที่ความน่าจะเป็นสูงที่สุด
    print(f'ν = {ν}\nV = {V}\nΣ = {Σ_max}')

ผลที่ได้
ν = 5
V = [[ 2.08333333 -0.41666667]
 [-0.41666667  2.08333333]]
Σ = [[0.25 0.05]
 [0.05 0.25]]

ν = 41
V = [[0.0055187  0.00104572]
 [0.00104572 0.01859767]]
Σ = [[ 4.81982819 -0.2710123 ]
 [-0.2710123   1.4302435 ]]

ν = 149
V = [[ 0.00176174 -0.00041632]
 [-0.00041632  0.00410707]]
Σ = [[3.98321883 0.40377066]
 [0.40377066 1.70861973]]

ν = 329
V = [[ 0.00068913 -0.00027996]
 [-0.00027996  0.00172312]]
Σ = [[4.76578006 0.77430338]
 [0.77430338 1.90599849]]

ν = 581
V = [[ 0.00036217 -0.00017873]
 [-0.00017873  0.00095164]]
Σ = [[5.26499485 0.98885953]
 [0.98885953 2.00375218]]

ν = 905
V = [[ 0.00024671 -0.00011678]
 [-0.00011678  0.00061443]]
Σ = [[4.93794367 0.93853582]
 [0.93853582 1.98272111]]

ν = 1301
V = [[ 1.72833326e-04 -7.97935605e-05]
 [-7.97935605e-05  4.19152283e-04]]
Σ = [[4.88708963 0.93034989]
 [0.93034989 2.0151434 ]]

ν = 1769
V = [[ 1.22722015e-04 -5.43669747e-05]
 [-5.43669747e-05  3.06883104e-04]]
Σ = [[5.00706708 0.88704489]
 [0.88704489 2.00231735]]

ν = 2309
V = [[ 9.57755868e-05 -4.67090851e-05]
 [-4.67090851e-05  2.38106922e-04]]
Σ = [[5.00678507 0.98217367]
 [0.98217367 2.01391784]]

ν = 2921
V = [[ 7.47142508e-05 -3.57389306e-05]
 [-3.57389306e-05  1.87296108e-04]]
Σ = [[5.04752753 0.9631446 ]
 [0.9631446  2.01350814]]

ν = 3605
V = [[ 6.11264567e-05 -2.97188289e-05]
 [-2.97188289e-05  1.53627692e-04]]
Σ = [[5.01329716 0.96980772]
 [0.96980772 1.99472562]]

ตัวแปรที่แจกแจงในการแจกแจงวิชาร์ตคือเมทริกซ์ความเที่ยงตรง แต่ในที่นี้เพื่อให้เห็นภาพง่ายจึงแสดงเป็นเมทริกซ์ความแปรปรวนเป็นหลัก เวลาคำนวณก็ใช้ np.linalg.inv เพื่อคำนวณเมทริกซ์ผกผัน

วาดภาพเคลื่อนไหวแสดงการแจกแจงที่เปลี่ยนไปตามเวลาได้แบบนี้



ภาพด้านบนคือค่าความหนาแน่นความน่าจะเป็นที่ค่า ต่างๆที่ สูงสุด โดยสีแดงคือค่าสูงสีม่วงคือค่าน้อย ค่าทั้ง 3 ตัวที่แสดงด้านบนคือค่าที่ตำแหน่งสูงสุด จุดรูปดาวแสดงตำแหน่งที่ค่าสูงสุด

ภาพด้านล่างแสดงการแจกแจงของจุดข้อมูลที่มี ทั้งหมด n จุด

ยิ่งใส่จุดข้อมูลมาก การแจกแจงค่า ที่คำนวณได้ก็จะยิ่งมากองอยู่ใกล้จุด ซึ่งเป็นค่าจริงๆของการแจกแจงนี้



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



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

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

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

หมวดหมู่

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

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

สารบัญ

รวมคำแปลวลีเด็ดจากญี่ปุ่น
มอดูลต่างๆ
-- numpy
-- matplotlib

-- pandas
-- manim
-- opencv
-- pyqt
-- pytorch
การเรียนรู้ของเครื่อง
-- โครงข่าย
     ประสาทเทียม
ภาษา javascript
ภาษา mongol
ภาษาศาสตร์
maya
ความน่าจะเป็น
บันทึกในญี่ปุ่น
บันทึกในจีน
-- บันทึกในปักกิ่ง
-- บันทึกในฮ่องกง
-- บันทึกในมาเก๊า
บันทึกในไต้หวัน
บันทึกในยุโรปเหนือ
บันทึกในประเทศอื่นๆ
qiita
บทความอื่นๆ

บทความแบ่งตามหมวด



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

  ค้นหาบทความ

  บทความแนะนำ

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

ไทย

日本語

中文