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



[python] การสร้างค่าสุ่มด้วยการแจกแจงแบบปกติหลายตัวแปร
เขียนเมื่อ 2018/05/25 19:23
ในบทความก่อนหน้านี้ได้อธิบายเรื่องการอธิบายความสัมพันธ์ระหว่างตัวแปรด้วยค่าความแปรปรวนร่วมเกี่ยว (协方差, covariance) และสหสัมพันธ์ (相关系数, correlation) ไป
https://phyblas.hinaboshi.com/20180517

สำหรับในบทความนี้จะพูดถึงเรื่องหนึ่งที่นำมาประยุกต์ใช้ได้ คือใช้สร้างข้อมูลที่มีการแจกแจงแบบปกติหลายตัวแปร (多元正态分布, multivariate normal distribution)

โดยพื้นฐานแล้ว numpy มีฟังก์ชันสำหรับสร้างข้อมูลสุ่มอยู่มากมาย ส่วนหนึ่งได้เขียนถึงไปแล้วใน https://phyblas.hinaboshi.com/numa15

แต่ในนี้จะขอแนะนำฟังก์ชัน np.random.multivariate_normal() ซึ่งไว้ใช้สร้างการแจกแจงแบบปกติหลายตัวแปร

การจะใช้ได้จำเป็นต้องอาศัยความเข้าใจเรื่องความแปรปรวนร่วมเกี่ยวและสหสัมพันธ์



การแจกแจงแบบปกติตัวแปรเดียว

ก่อนอื่นมาทวนเรื่องการแจกแจงแบบปกติตัวแปรเดียวก่อน

การแจกแจงแบบปกติ (正态分布, normal distribution) หรืออาจเรียกว่า การแจกแจงแบบเกาส์ (高斯分布, Gaussian distribution) เป็นการแจกแจงในรูปแบบที่ใช้บ่อยที่สุด

การแจกแจงแบบปกติเป็นการแจกแจงที่ขึ้นอยู่กับพารามิเตอร์ ๒​ ตัวที่ต้องกำหนดคือ
- μ คือ จุดกึ่งกลาง
- σ คือ ส่วนเบี่ยงเบนมาตรฐาน

โดยความหนาแน่นของการแจกแจงจะเป็นไปตามสมการนี้
..(1)

หากเขียนเป็นฟังก์ชันในไพธอนแล้ววาดกราฟจะเป็นแบบนี้
import numpy as np
import matplotlib.pyplot as plt
def pakati(x,mu,sigma):
    return np.exp(-(x-mu)**2/(2*sigma**2))/np.sqrt(2*np.pi)/sigma
x = np.linspace(0,2,101)
y = pakati(x,1,0.3)
plt.plot(x,y,'r')
plt.show()

กราฟเป็นรูประฆังคว่ำ ค่าจะมากสุดที่จุด μ และค่อยๆลดลงเมื่อยิ่งห่างออกไป อัตราการลดลงจะเร็วแค่ไหนขึ้นอยู่กับ σ



หากสร้างค่าสุ่มให้แจกแจงแบบปกติโดยกำหนดค่า μ และ σ ที่ค่าหนึ่ง เมื่อนำมาหาค่าเฉลี่ยก็จะได้ประมาณ μ ส่วนเบี่ยงเบนมาตรฐานก็จะได้ประมาณ σ

ในไพธอนสามารถสร้างการแจกแจงแบบปกติได้ด้วยฟังก์ชัน np.random.normal พารามิเตอร์คือ (μ, σ, จำนวนที่ต้องการสุ่ม)
s = np.random.normal(1,0.3,10000)
plt.hist(s,100,color='r')
plt.show()

print(s.mean()) # ได้ 0.9987304195735669
print(s.std()) # ได้ 0.2962901261818392


หากลองสร้างสองตัวแปรที่แจกแจงแบบปกติพร้อมกันแล้วหาความแปรปรวนร่วมเกี่ยวดู
x = np.random.normal(12,2,100000)
y = np.random.normal(27,1,100000)
plt.gca(aspect=1,xlim=[x.min(),x.max()],ylim=[y.min(),y.max()])
plt.scatter(x,y,c='r',edgecolors='k',alpha=0.01)
plt.show()
print(np.cov(x,y))

จะได้ว่าค่าในแนวทแยงก็คือค่า σ กำลังสอง ส่วนขวาบนซ้ายล่างเข้าใกล้ 0 นั่นเพราะทั้ง ๒ ตัวแปรถูกสุ่มอย่างอิสระไม่ได้เกี่ยวข้องต่อกันเลย

แต่ว่าในธรรมชาติโดยทั่วไปแล้ว ตัวแปรสองตัวที่มีการแจกแจงแบบปกติอาจมีความเกี่ยวพันกัน เช่นพอตัวนึงมากอีกตัวนึงก็มีแนวโน้มมากหรือน้อยตาม

กรณีแบบนี้จะต้องสร้างตัวแปรทั้งสองพร้อมกันด้วยการแจกแจงแบบปกติหลายตัวแปร



การแจกแจงแบบปกติหลายตัวแปร

หากมีตัวแปร x และ y ที่ต้องการให้แจกแจงแบบปกติโดยมีความเกี่ยวเนื่องกัน อาจเขียนสมการการแจกแจงได้ดังนี้
..(2)

โดย ρxy คือค่าสหสัมพันธ์ของ x และ y
..(3)

จะเห็นว่าเดิมทีที่เป็นแค่ (x-μx)2 ได้ถูกเปลี่ยนมาเป็นก้อนที่ซับซ้อนมากขึ้น มีทั้งส่วนร่วมของ x และ y อยู่ด้วยกัน

สมการอาจเขียนเป็นรูปทั่วไปด้วยเมทริกซ์ได้แบบนี้
..(4)

โดย x ในที่นี้ใช้ตัวหนา หมายถึงเป็นปริมาณเวกเตอร์

Σ คือเมทริกซ์ความแปรปรวนร่วมเกี่ยว โดย |Σ| คือดีเทอร์มิแนนต์ (determinant) ของเมทริกซ์ Σ

กรณีที่มีตัวแปรแค่ x และ y นั้น Σ จะเป็นแบบนี้
..(5)

กรณีที่มีหลายตัวแปรจะได้ว่า
..(6)

ส่วน x-μ ในที่นี้คือเวกเตอร์ที่แสดงถึงระยะห่างจากจุดกึ่งกลางของตัวแปรต่างๆ เมื่อนำมาคูณประกบกันกับเมทริกซ์ความแปรปรวนร่วมเกี่ยวในลักษณะแบบนี้สิ่งที่ได้ก็คือสิ่งที่เรียกว่า ระยะทางมหาลโนพิส (Mahalanobis distance)
..(7)

ชื่อที่ฟังดูยืดยาวจำยากนี้ถูกตั้งตามนักคณิตศาสตร์ชาวอินเดียชื่อ ประศานตะ จันทระ มหาลโนพิส (प्रशान्त चन्द्र महालनोबिस)

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

ในไพธอนเมทริกซ์จะเขียนอยู่ในรูปตัวแปรชนิดอาเรย์ของ numpy ส่วนการคูณกันของเว็กเตอร์ จะใช้คำสั่ง .dot()

ส่วนการหาดีเทอร์มิแนนท์ก็มีคำสั่ง np.linalg.det() การหาเมทริกซ์ผกผันก็ใช้ np.linalg.inv()

เราอาจเขียนโปรแกรมให้คำนวณตามสมการ (4) ได้ดังนี้
def pakati(X,mu,cov):
    d = X-mu
    k = len(d)
    det = np.linalg.det(cov)
    inv = np.linalg.inv(cov)
    m = d.T.dot(inv).dot(d) # กำลังสองของระยะห่าง mahalanobis
    return np.exp(-m/2)/np.sqrt((2*np.pi)**k*det)

นอกจากนี้ใน scipy ได้เตรียมฟังก์ชันแบบนี้ไว้ คือ scipy.stats.multivariate_normal.pdf ใช้อันนี้ได้เลย

ลองนำมาใช้วาดแผนภาพแสดงการแจกแจงในสองมิติดู โดยแสดงผลดูทั้งแบบสามมิติและแบบคอนทัวร์
from mpl_toolkits.mplot3d import Axes3D
from scipy.stats import multivariate_normal

mu = np.array([0,0])
cov = np.array([[1.2,0.8],[0.8,0.9]])
mx,my = np.meshgrid(np.linspace(-4,4,51),np.linspace(-4,4,51))
mX = np.stack([mx,my],2)
mz = multivariate_normal.pdf(mX,mu,cov)
ax = plt.figure(figsize=[6,6]).add_axes([0,0,1,1],projection='3d')
ax.plot_surface(mx,my,mz,rstride=1,cstride=1,alpha=0.2,edgecolor='k',cmap='rainbow')
plt.figure().gca(aspect=1)
plt.contourf(mx,my,mz,40,cmap='rainbow')
plt.colorbar()
plt.show()




ส่วนการสร้างค่าสุ่มตามการแจกแจงแบบปกติหลายตัวแปรในไพธอนมีฟังก์ชัน np.random.multivariate_normal

วิธีใช้คือ np.random.multivariate_normal(จุดกึ่งกลาง, เมทริกซ์ความแปรปรวนร่วมเกี่ยว, จำนวนที่ต้องการสุ่ม)

ตัวอย่าง
mu = np.array([10,5])
cov = np.array([[1.8,1.1],[1.1,1.5]])
X = np.random.multivariate_normal(mu,cov,20000)
plt.figure().gca(aspect=1)
plt.scatter(X[:,0],X[:,1],edgecolors='k',color='r',alpha=0.01)
plt.show()




การสร้างเมทริกซ์ความแปรปรวนร่วมเกี่ยวด้วยวิธีเคอร์เนล

เพื่อที่จะสร้างข้อมูลสุ่มที่มีการแจกแจงแบบปกติหลายตัวแปร เราจำเป็นต้องมีเมทริกซ์ความแปรปรวนร่วมเกี่ยว

ปกติถ้าเราสร้างเมทริกซ์ความแปรปรวนร่วมเกี่ยวขึ้นมาจากการคำนวณ เช่น มีชุดข้อมูลอยู่แล้วนำมาคำนวณด้วยฟังก์ชัน np.cov แบบนั้นแน่นอนว่าเมทริกซ์ที่ได้สามารถนำมาใช้ได้

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

วิธีหนึ่งที่มักใช้ในการสร้างเมทริกซ์ความแปรปรวนร่วมเกี่ยวให้ได้ตามที่ต้องการคือ สร้างจากเคอร์เนล (kernel) คือการนำสมาชิกของตัวแปรชุดนึงทั้งหมดมาจับคู่กระทำอะไรบางอย่างกันแล้วเขียนให้อยู่ในรูปเมทริกซ์
..(8)

โดยในที่นี้ k คือฟังก์ชันที่นำตัวแปร ๒ ตัวมาใช้คำนวณ มักเรียกกันว่าเคอร์เนล

เมทริกซ์ที่ได้จากการจับคู่ทำอะไรกันในลักษณะนี้มีชื่อเรียกว่า เมทริกซ์กรัม (格拉姆矩阵, Gramian matrix) เมทริกซ์ความแปรปรวนร่วมเกี่ยวนั้นโดยทั่วไปถือว่าเป็นเมทริกซ์กรัมชนิดหนึ่ง

ฟังก์ชันที่มักใช้เป็นเคอร์เนลนั้นมีอยู่หลากหลายรูปแบบ

เช่นแบบที่ง่ายที่สุดคือ มีค่าค่าหนึ่งเฉพาะเมื่อ x และ y เท่ากัน คือ
..(9)

โดยในที่นี้ a คือค่าส่วนเบี่ยงเบนมาตรฐาน

เขียนในไพธอนได้แบบนี้
def kernel(x,a):
    x1,x2 = np.meshgrid(x,x)
    return a**2*(x1==x2)

เมื่อใช้คำนวณก็จะได้ผลที่เป็นแบบที่เรียบง่ายที่สุดก็คือเมทริกซ์แนวทแยง คือมีค่าแต่ในแนวทแยงนอกนั้น 0 ซึ่งจะหมายถึงว่าตัวแปรทั้งหมดมีการแจกแจงในตัวเองที่เป็นอิสระจากตัวแปรอื่นใด

..(10)

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

เช่น
def kernel(x,a):
    x1,x2 = np.meshgrid(x,x)
    return a**2*(x1==x2)

x = np.linspace(0,100,101)
mu = np.zeros(101)
cov = kernel(x,1)
y = np.random.multivariate_normal(mu,cov)
plt.plot(x,y,'r-o')
plt.show()




เคอร์เนลกำลังสองเอ็กซ์โพเนนเชียล

ในการสร้างข้อมูลที่เป็นอนุกรมเวลาที่ค่าในจุดนึงสัมพันธ์กับจุดก่อนและหลัง เคอร์เนลแบบหนึ่งที่นิยมใช้กันมากก็คือกำลังสองเอ็กซ์โพเนนเชียล (squared-exponential)
..(11)
..(12)

โดย a เป็นตัวกำหนดขนาดการแจกแจง (ส่วนเบี่ยงเบนมาตรฐาน) ส่วน s กำหนดขนาดระยะที่ความเกี่ยวพันจะส่งผลไปถึง

สามารถสร้างฟังก์ชันเคอร์เนลตามสมการได้ดังนี้
def kernel(x,a,s):
    x1,x2 = np.meshgrid(x,x)
    return a**2*np.exp(-0.5*((x1-x2)/s)**2)

แล้วลองนำเคอร์เนลมาใช้สร้างเมทริกซ์ความแปรปรวนร่วมเกี่ยวสำหรับข้อมูลอนุกรมเวลา
x = np.linspace(0,100,101)
cov = kernel(x,2,5)
plt.imshow(cov,cmap='afmhot')
plt.colorbar()
plt.show()

จะได้ค่าสูงสุดที่แนวทแยง โดยมีค่าเป็น a ยกกำลังสอง และค่อยๆลดลงเมื่อห่างออกไป



จากนั้นก็นำมาใช้ สุ่มสร้างอนุกรมเวลา โดยลองเทียบค่า s ต่างกัน สร้างดูสักแบบละ ๕ เส้น
plt.figure(figsize=[5,6])
x = np.linspace(0,100,201)
mu = np.zeros(201)
for i in range(1,4):
    s = 3**i
    cov = kernel(x,10,s)
    y = np.random.multivariate_normal(mu,cov,5)
    plt.subplot(3,1,i)
    plt.plot(x,y.T)
    plt.ylabel('s=%d'%s)
plt.tight_layout()
plt.show()

จะเห็นว่าแต่ละเส้นที่ได้มีลักษณะที่มีการเปลี่ยนแปลงไปทีละเล็กน้อย ค่าในจุดนึงจะขึ้นอยู่กับจุดก่อนหลัง ไม่ได้โดดไปไกลนัก และยิ่งค่า s มากขึ้น ความเปลี่ยนแปลงก็จะยิ่งน้อย



นอกจากนี้เรายังอาจเอาเคอร์เนลมารวมกันได้ เช่น
..(13)
def kernel(x,a1,s,a2):
    x1,x2 = np.meshgrid(x,x)
    return a1**2*np.exp(-0.5*((x1-x2)/s)**2) + a2**2*(x1==x2)

แบบนี้จะได้การเปลี่ยนแปลงที่โดยรวมแล้วมีลักษณะค่อยเป็นค่อยไปตามเวลา แต่ก็มีการเปลี่ยนแปลงขึ้นลงจากความไม่แน่นอนภายในตัว
x = np.linspace(0,100,501)
mu = np.zeros(501)
cov = kernel(x,20,15,2)
y = np.random.multivariate_normal(mu,cov,8)
plt.plot(x,y.T)
plt.show()




การเปลี่ยนแปลงแบบเป็นคาบ

บางครั้งเราอาจต้องการสุ่มสร้างค่าบางอย่างที่มีการเปลี่ยนแปลงแบบเป็นคาบ (periodic) คือเมื่อถึงจุดนึงจะกลับมามีค่าเท่าเดิม กรณีแบบนี้นิยมใช้เคอร์เนลแบบไซน์กำลังสองเอ็กซ์โพเนนเชียล (exp-sine-squared kernel)
..(14)
def kernel(x,a,tau,l):
    x1,x2 = np.meshgrid(x,x)
    return a**2*np.exp(-2*(np.sin((x1-x2)*np.pi/tau)/l)**2)

ในที่นี้ τ คือคาบ คือเป็นตัวกำหนดว่าไปไกลแค่ไหนจะกลับมาวนซ้ำ

ส่วน s คือตัวกำหนดว่าจะมีความเปลี่ยนแปลงขึ้นลงบ่อยแค่ไหนในคาบ ยิ่งน้อยยิ่งขึ้นลงมาก

ลองดูค่าในเมทริกซ์ความแปรปรวนร่วมเกี่ยวที่ได้จากการคำนวณ
x = np.linspace(0,100,101)
cov = kernel(x,1,20,1)
plt.imshow(cov,cmap='afmhot',vmin=0,vmax=1)
plt.colorbar()
plt.show()


ลองนำมาใช้แล้ววาดกราฟดู
x = np.linspace(0,100,501)
mu = np.zeros(501)
plt.figure(figsize=[6,6])
for i in range(1,4):
    l = 4**(i-2)
    cov = kernel(x,10,40,l)
    y = np.random.multivariate_normal(mu,cov,5)
    plt.subplot(3,1,i)
    plt.plot(x,y.T)
    plt.ylabel(r'$\lambda$=%.2f'%l)
plt.tight_layout()
plt.show()




การเปลี่ยนแปลงแบบเกือบเป็นคาบ

บางครั้งเราอาจต้องการได้ค่าที่มีการเปลี่ยนแปลงเป็นคาบ แต่ไม่ได้เหมือนเดิมตลอดทุกคาบ แต่ในแต่ละคาบมีการเปลี่ยนแปลงไปทีละน้อย ลักษณะแบบนั้นเรียกว่า quasi-periodic ซึ่งอาจแปลเป็นไทยว่า "เกือบเป็นคาบ"

เราสามารถทำเคอร์เนลในลักษณะนี้ได้โดยเอาเคอร์เนลแบบไซน์กำลังสองเอ็กซ์โพเนนเชียลมาคูณกับกำลังสองเอ็กซ์โพเนนเชียล กลายเป็นแบบนี้
..(15)

ตัวอย่าง
def kernel(x,a,tau,l,s):
    x1,x2 = np.meshgrid(x,x)
    dx = x1-x2
    return a**2*np.exp(-2*(np.sin(dx*np.pi/tau)/l)**2-0.5*(dx/s)**2)

x = np.linspace(0,100,501)
mu = np.zeros(501)
plt.figure(figsize=[6,6])
for i in range(1,4):
    s = 3**i*10
    cov = kernel(x,10,20,2,s)
    y = np.random.multivariate_normal(mu,cov,5)
    plt.subplot(3,1,i)
    plt.plot(x,y.T)
    plt.ylabel(r's=%.2f'%s)
plt.tight_layout()
plt.show()


ค่าจะกลับมาเหมือนเดิมทุก 20 แต่ไม่ได้เหมือนกันเป๊ะ แต่ยิ่ง s มาก ความแตกต่างของแต่ละรอบก็จะยิ่งน้อยลง



นี่เป็นแค่คัวอย่างส่วนหนึ่งเท่านั้น นอกจากนี้ยังมีรูปแบบการใช้อยู่อีกมากมายสามารถไปลองศึกษาเพิ่มเติมกันได้

สรุปเนื้อหาทั้งหมดนี้น่าจะทำให้เข้าใจวิธีการสร้างชุดข้อมูลจากการสุ่มค่าด้วยฟังก์ชัน np.random.multivariate_normal ตามที่ต้องการกันได้ไม่มากก็น้อยแล้ว


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

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

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

หมวดหมู่

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

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

สารบัญ

รวมคำแปลวลีเด็ดจากญี่ปุ่น
python
-- numpy
-- matplotlib

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

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



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

  ค้นหาบทความ

  บทความแนะนำ

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

บทความแต่ละเดือน

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月

2016年

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

2015年

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

ค้นบทความเก่ากว่านั้น

ไทย

日本語

中文