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



[python] วิเคราะห์การถดถอยพหุนาม
เขียนเมื่อ 2016/12/19 10:43
ก่อนหน้านี้ได้เขียนถึงการวิเคราะห์การถดถอยเชิงเส้นในหนึ่งมิติไป https://phyblas.hinaboshi.com/20161210

แต่ว่าการวิเคราะห์การถดถอยเชิงเส้นนั้นสามารถแก้ได้แต่ปัญหาที่เป็นเชิงเส้น นั่นคืออยู่ในรูป z = w*x + b

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


แบบนี้ย่อมไม่ควรจะใช้การถดถอยเชิงเส้นได้

หากต้องการปัญหาที่ผลเฉลยไม่น่าจะเป็นเชิงเส้น ก็ต้องปรับวิธีใหม่ซึ่งซับซ้อนขึ้น

รูปแบบคำตอบที่ง่ายรองลงมาจากเชิงเส้นธรรมดาก็คืออยู่ในรูปพหุนาม นั่นคือเป็นลักษณะแบบนี้
z = w0 + w1*x + w2*x**2 + w3*x**3 + ...

รูปแบบสมการลักษณะนี้ดูแล้วก็จะคล้ายกับสมการเชิงเส้นหลายมิติ ซึ่งเขียนถึงไปใน https://phyblas.hinaboshi.com/20161212

เพียงแต่ว่าเปลี่ยนจาก x1 เป็น x**1 ส่วน x2 เป็น x**2 และ xn เป็น x**n

นั่นหมายความว่าถ้าใช้แนวคิดเดียวกันเราก็สามารถแก้ปัญหาพหุนามได้ โดยให้ถือว่า x, x**2 และ x**3 ล้วนเป็นตัวแปรต้นคนละตัว โดย x1 เป็น x และ x2 เป็น x**2 แบบนี้

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

ก่อนอื่นเริ่มจากเฉลยว่าข้อมูลที่วาดในรูปข้างต้นมาจากโค้ดนี้
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0,1,101)
z = 1+3*x-7*x**2+5*x**3+0.1*np.random.randn(101)
plt.scatter(x,z)
plt.show()

ซึ่งดูจากตรงนี้จะเห็นว่าคำตอบที่ควรจะเป็นก็คือ w0=5, w1=21, w2=31 และ w3=1

เช่นเคย คราวนี้เราก็จะมาดูว่าวิธีการที่เราใช้จะสามารถหาคำตอบได้ตามนี้หรือเปล่า

ก่อนอื่นที่น่าพูดถึงคือวิธีการสร้างข้อมูลตัวแปรต้น โดยเดิมทีเรามีแค่ x ตัวเดียว แต่ว่าต้องสร้างค่า x ที่ยกกำลังต่างๆตั้งแต่ 0 ไปจนถึงดีกรีที่ต้องการ

อาจใช้วิธีการดังนี้
deg = 3 # ดีกรีที่ต้องการ
X = np.stack([x**i for i in range(deg+1)],1)
print(X[::20])

ได้
[[ 1.     0.     0.     0.   ]
 [ 1.     0.2    0.04   0.008]
 [ 1.     0.4    0.16   0.064]
 [ 1.     0.6    0.36   0.216]
 [ 1.     0.8    0.64   0.512]
 [ 1.     1.     1.     1.   ]]

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

ใช้วิธีการนี้เราสามารถเขียนคลาสสำหรับการวิเคราะห์การถดถอยพหุนามแบบตัวแปรเดียวได้ดังนี้
class ThotthoiPhahunam1D:
    def __init__(self,eta,deg=2):
        self.eta = eta
        self.deg = deg

    def rianru(self,x,z,d_yut=1e-7,n_thamsam=100000):
        X = np.stack([x**i for i in range(self.deg+1)],1) # ทำให้เป็นพหุนาม
        self.w = np.zeros(X.shape[1])
        dw = np.zeros(X.shape[1])
        h = np.dot(X,self.w)
        self.sse = [self.ha_sse(h,z)]
        for i in range(n_thamsam):
            eee = 2*(z-h)*self.eta
            dw = np.dot(eee,X)
            self.w += dw
            h = np.dot(X,self.w)
            self.sse += [self.ha_sse(h,z)]
            if(np.all(np.abs(dw)<d_yut)):
                break

    def thamnai(self,x):
        X = np.stack([x**i for i in range(self.deg+1)],1)
        return np.dot(X,self.w)

    def ha_sse(self,h,z):
        return ((h-z)**2).sum()

โดยในนี้เราให้กำหนดค่าดีกรีในพารามิเตอร์ deg ตอนที่สร้างออบเจ็กต์พร้อมกับค่าอัตราการเรียนรู้

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

ลองนำมาใช้กับค่า x และ z ของโจทย์ข้อนี้ที่เราเตรียมเอาไว้ โดยสร้างออบเจ็กต์จากนั้นใช้เมธอด rianru
eta = 0.005
deg = 3
tp = ThotthoiPhahunam1D(eta,deg)
tp.rianru(x,z)

จากนั้นลองดูค่าน้ำหนักที่ได้
print(tp.w)

ก็จะได้ค่าน้ำหนักใกล้เคียงกับที่เราใส่ตอนสร้างข้อมูลขึ้น
[ 1.06250237  2.51213577 -5.8278461   4.21149162]

ทีนี้ลองวาดกราฟดูเทียบกับจุดข้อมูลที่ใช้เรียนรู้ โดยใช้เมธอด thamnai เพื่อหาค่าที่ทำนายโดยใช้ค่าน้ำหนักที่ได้
plt.scatter(x,z)
x = np.linspace(0,1,101)
h = tp.thamnai(x)
plt.plot(x,h,'b')
plt.show()



เส้นที่ทำนายขึ้นลากผ่านกลางจุดข้อมูลอย่างสวยงาม แสดงให้เห็นว่าการคำนวณถูกต้อง



ต่อมาลองคิดกรณีที่ปัญหาประกอบด้วย ๒ ตัวแปรต้นขึ้นไป ซึ่งจะซับซ้อนไปอีก

เช่นลองกรณีของสองมิติดังตัวอย่างนี้
x = np.random.uniform(0,4,500)
y = np.random.uniform(0,4,500)
z = 22 + 4*x + 7*y + x*x - 6*x*y + y*y + np.random.randn(500)

ซึ่งประกอบด้วย ๖ พจน์ คือ 1,x,y,xx,xy,yy ตามลำดับ

กรณีแบบนี้เรียกว่าเป็นพหุนามดีกรี ๒ คือผลรวมเลขชี้กำลังสูงสุดไม่เกิน ๒ ถ้ามี x*x*x หรือ x*y*y อยู่ด้วยก็จะเป็นดีกรี ๓ เป็นต้น

ลองนำมาวาดภาพการกระจายในสามมิติดู
from mpl_toolkits.mplot3d import Axes3D
plt.figure(figsize=[8,8])
ax = plt.axes([0,0,1,1],projection='3d',xlim=[0,4],ylim=[0,4],zlim=[z.min(),z.max()])
for i in range(len(z)):
    ax.plot([x[i],x[i]],[y[i],y[i]],[z[i],0],c=plt.get_cmap('jet')((z[i]-z.min())/(z.max()-z.min())),alpha=0.2)
ax.scatter(x,y,z,c=z,edgecolor='k')
plt.show()



วิธีการสร้างตัวแปรต้นสำหรับมาใช้วิเคราะห์ถดถอยพหุนามนั้น สำหรับกรณีสองตัวแปรเราอาจเขียนแจกแจงพหุนามได้ดังนี้
deg = 2
xy_ = np.stack([x**(k-i)*y**i for k in range(deg+1) for i in range(k+1)],1)
print(xy_)

***ใครงงกับรูปแบบการเขียนโค้ดตรงนี้อาจลองอ่านทวนเนื้อหา ไพธอนเบื้องต้น บทที่ ๙: การทำซ้ำด้วย for ตรงส่วนของการใช้ for ซ้อน for

ได้
[[  1.00000000e+00   2.99032562e+00   1.22752793e+00   8.94204731e+00
    3.67070821e+00   1.50682481e+00]
 [  1.00000000e+00   1.51777983e+00   2.29233098e+00   2.30365561e+00
    3.47925373e+00   5.25478133e+00]
 [  1.00000000e+00   9.15271919e-01   1.12181126e-01   8.37722686e-01
    1.02676235e-01   1.25846051e-02]
 ...,
 [  1.00000000e+00   1.39509132e+00   2.08564151e+00   1.94627979e+00
    2.90966036e+00   4.34990050e+00]
 [  1.00000000e+00   3.75774494e+00   2.03009611e+00   1.41206470e+01
    7.62858336e+00   4.12129020e+00]
 [  1.00000000e+00   1.40287884e-02   7.70100843e-01   1.96806904e-04
    1.08035818e-02   5.93055309e-01]]

โดยกรณีนี้เป็นดีกรี ๒ ค่าที่ได้แต่ละหลักจะเป็นพจน์ของ 1,x,y,xx,xy,yy ตามลำดับ

และถ้าหากเป็นดีกรี ๓ ก็จะมีพจน์ xxx,xxy,xyy,yyy เพิ่มเข้ามาอีก รวมเป็น ๑๐ หลัก เป็นต้น

การเขียนแบบนี้ดูแล้วยุ่งยากพอสมควร อีกทั้งถ้ามี ๓ ตัวแปรการเขียนก็จะยิ่งยุ่งยากขึ้นไปอีก

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

ในมอดูลย่อย preprocessing ของ sklearn มีคลาสที่ชื่อ PolynomialFeatures อยู่ มีไว้สำหรับทำในสิ่งเดียวกับที่เราเพิ่งทำไป

ขออธิบายด้วยการยกตัวอย่างการใช้งานจริงก่อน

การใช้ลองเขียนตามนี้
from sklearn.preprocessing import PolynomialFeatures as Pf
xy = np.stack([x,y],1)
xy_ = Pf(deg).fit_transform(xy)

ก็จะได้ xy_ แบบเดียวกับที่ใช้วิธีอย่างในตอนแรก

การใช้ PolynomialFeatures จะคล้ายกับ StandardScaler ที่กล่าวถึงไปในเรื่องของการทำข้อมูลให้เป็นมาตรฐาน https://phyblas.hinaboshi.com/20161124

คือจะใช้เมธอด fit_transform ในลักษณะเดียวกัน นั่นคือเป็นการย่อขั้นตอนการใช้ เพื่อให้ไม่ต้องใช้ fit แล้วต่อด้วย transform อีกที หากเขียนแบบเต็มๆน่าจะเป็นแบบนี้
pf = Pf(deg)
pf.fit(xy)
xy_ = pf.transform(xy)

เมื่อใช้ PolynomialFeatures แบบนี้แล้วจะใช้กับปัญหาที่มีตัวแปรต้นกี่ตัวก็ได้ เช่นถ้าลองใช้กับ ๓​ ตัวแปรดู แจกแจงถึงดีกรี ๒
x1 = np.linspace(0,2,3)
x2 = np.linspace(0,3,3)
x3 = np.linspace(0,5,3)
print(Pf(2).fit_transform(np.stack([x1,x2,x3],1)))

ได้
[[  1.     0.     0.     0.     0.     0.     0.     0.     0.     0.  ]
 [  1.     1.     1.5    2.5    1.     1.5    2.5    2.25   3.75   6.25]
 [  1.     2.     3.     5.     4.     6.    10.     9.    15.    25.  ]]

ผลที่ได้มี ๑๐ พจน์ ตามนี้ 1, x1, x2, x3, x1x1, x1x2, x1x3, x2x2, x2x3, x3x3

เมื่อเข้าใจวิธีการแปลงตัวแปรต้นให้เป็นพหุนามแล้ว ก็ลองสร้างแบบจำลองขึ้นดู เขียนได้ดังนี้
import numpy as np
from sklearn.preprocessing import PolynomialFeatures as Pf

class ThotthoiPhahunam:
    def __init__(self,eta,deg=2):
        self.eta = eta
        self.deg = deg

    def rianru(self,X,z,d_yut=1e-7,n_thamsam=100000):
        X = Pf(self.deg).fit_transform(X)
        self.w = np.zeros(X.shape[1])
        dw = np.zeros(X.shape[1])
        h = np.dot(X,self.w)
        self.sse = [self.ha_sse(h,z)]
        for i in range(n_thamsam):
            eee = 2*(z-h)*self.eta
            dw = np.dot(eee,X)
            self.w += dw
            h = np.dot(X,self.w)
            self.sse += [self.ha_sse(h,z)]
            if(np.all(np.abs(dw)<d_yut)):
                break

    def thamnai(self,X):
        X = Pf().fit_transform(X)
        return np.dot(X,self.w)

    def ha_sse(self,h,z):
        return ((h-z)**2).sum()

ลักษณะคล้ายกับกรณีตัวแปรเดียวแต่ว่าค่า X ที่ต้องใส่ตอนที่เรียนรู้คราวนี้เป็นอาเรย์สองมิติ และใช้ PolynomialFeatures ในการทำให้เป็นพหุนาม โดยค่า deg จะกำหนดดีกรีที่ต้องการแปลง

นำมาใช้กับ x,y,z ที่เพิ่งเตรียมไว้จากตัวอย่างที่แล้ว แล้วลองวาดผลที่ได้ดู
X = np.stack([x,y],1)
eta = 0.00001
deg = 2
tp = ThotthoiPhahunam(eta,deg)
tp.rianru(X,z)

plt.figure(figsize=[8,8])
ax = plt.axes([0,0,1,1],projection='3d',xlim=[0,4],ylim=[0,4],zlim=[z.min(),z.max()])
mx,my = np.meshgrid(np.linspace(0,4,11),np.linspace(0,4,11))
mX = np.stack([mx.ravel(),my.ravel()],1)
mz = tp.thamnai(mX).reshape(11,11)
ax.plot_surface(mx,my,mz,rstride=1,cstride=1,alpha=0.2,color='b',edgecolor='k')
h = tp.thamnai(X)
for i in range(len(z)):
    ax.plot([x[i],x[i]],[y[i],y[i]],[h[i],z[i]],'k')
ax.scatter(x,y,z,c=np.abs(z-h),edgecolor='k',cmap='jet')
plt.show()



พื้นผิวค่าที่ได้ลากผ่านกลางระหว่างจุดแสดงว่าแทนคำตอบได้ดี

ลองดูค่าน้ำหนักที่ได้
print(tp.w)

ได้ค่าออกมาใกล้เคียงกับ [22,4,7,1,-6,1] ซึ่งเป็นคำตอบที่ควรเป็นจริงๆ
[ 21.89569089   4.39975668   6.65906795   0.91058478  -6.01294355
   1.08767486]



อ้างอิง


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

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

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

หมวดหมู่

-- คอมพิวเตอร์ >> ปัญญาประดิษฐ์
-- คอมพิวเตอร์ >> เขียนโปรแกรม >> python >> matplotlib
-- คอมพิวเตอร์ >> เขียนโปรแกรม >> python >> numpy

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

สารบัญ

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

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

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



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

  ค้นหาบทความ

  บทความแนะนำ

หลักการเขียนทับศัพท์ภาษาจีนกลาง
g ในภาษาญี่ปุ่นออกเสียง "ก" หรือ "ง" กันแน่
ค้นพบระบบดาวเคราะห์ ๘ ดวง เบื้องหลังความสำเร็จคือปัญญาประดิษฐ์ (AI)
หอดูดาวโบราณปักกิ่ง ตอนที่ ๑: แท่นสังเกตการณ์และสวนดอกไม้
พิพิธภัณฑ์สถาปัตยกรรมโบราณปักกิ่ง
บ้านเก่าของจางเสวียเหลียงในเทียนจิน
เที่ยวจิ่นโจว ๓ วัน ๒ คืน 23 - 25 พ.ค. 2015
เที่ยวเมืองตานตง ล่องเรือในน่านน้ำเกาหลีเหนือ
บันทึกการเที่ยวสวีเดน 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月

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

ไทย

日本語

中文