Regresi dan Interpolasi
Modul #3 Praktikum
AS2205 Astronomi Komputasi
Oleh Dr. Muhamad Irfan Hakim
Program Studi Astronomi Fakultas Matematika dan Ilmu Pengetahuan Alam Institut Teknologi Bandung 2014
Daftar Isi Daftar Isi
i
1 Pendahuluan
1
2 Praktikum
4
2.1
Interpolasi Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2.2
Regresi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
i
Bab 1
Pendahuluan Untuk pekerjaan komputasi regresi dan interpolasi kali ini, ada beberapa hal yang perlu kita perhatikan: 1. Pekerjaan pertama praktikum ini adalah interpolasi, pekerjaan kedua adalah regresi linear dan kuadratik. 2. Pada interpolasi Newton, proses perhitungan interpolasi terdiri atas dua algoritma (a) algoritma perhitungan koefisien interpolan Newton m ← len(xData) # jumlah titik data a ← yData.copy() # salin koordinat-y dari data for k in range(1, m): for i in range(k, m): a[i] ← (a[i] - a[k-1])/(xData[i] - xData[k-1]) atau m ← len(xData) # jumlah titik data a ← yData.copy() # salin koordinat-y dari data for k in range(1, m): a[k:m] ← (a[k:m]-a[k-1])/(xData[k:m]-xData[k-1]) (b) algoritma perhitungan interpolan Newton n ← len(xData) - 1 # derajat maksimum polinom p ← a[n] # nilai awal polinom for k in range(1, n + 1): p ← a[n - k] + (x - xData[n - k])*p 3. Untuk dapat memanfaatkan modul python secara maksimal, modul array dan arange dapat dimanfaatkan. Contoh: from numpy import array,arange cx = array([2., 3., 4.])
# deklarasi cx sebagai array
1
Bab 1. Pendahuluan
# arange = array range, mirip dengan pemanfaatan range # namun dapat menggunakan bilangan riil for x in arange(0.0, 1.0+0.1, 0.1): print(x) # cek!! Ini dapat dimanfaatkan misalnya untuk tabulasi fungsi, atau plot data. 4. Pada interpolasi, variabel xData dan yData perlu dipasok sebagai input. Akan dipilih input dengan cara baca file. Misalkan file “data.csv” (6 baris): 0.15 2.30 3.15 4.85 6.25 7.95
4.79867 4.49013 4.22430 3.47313 2.66674 1.51909
File ini dapat dibaca sebagai berikut (silakan coba): import numpy as np fdata = open("data.csv", "r") # buka data.csv u/ dibaca m = input(’Number of data, m = ’) # jumlah baris di data.csv xData = np.zeros(m) yData = np.zeros(m)
# nilai awal nol u/ xData # nilai awal nol u/ yData
for j in range(0, m): # # # # #
data kolom-1 dimulai dari byte 0, dibaca 4 karakter data kolom-2 dimulai dari byte 5, dibaca 7 karakter tiap baris ada 12 karakter (termasuk karakter ganti baris (invisible character)). Maka, pembacaan berikutnya dimulai dari byte 13 (dan seterusnya).
fdata.seek(0+13*j); x = fdata.read(4); xData[j] = x fdata.seek(5+13*j); y = fdata.read(7); yData[j] = y fdata.close() # print(xData) # print(yData)
# jangan lupa tutup file # u/ debug # u/ debug
5. Untuk regresi, diperlukan solusi SPL. Dengan melihat rujukan modul praktikum sebelumnya, solusi via algoritma milik numpy dapat dilihat dalam contoh di bawah ini: # Solusi numerik SPL: # BA = C Institut Teknologi Bandung
2
Bab 1. Pendahuluan
# menggunakan algoritma SPL dari numpy import numpy as np from numpy.linalg import * # matriks B untuk contoh B = np.array([ [ 2., 3., -1.], [ 4., 4., -3.], [-2., 3., -1.] ]) # C A #
matriks C untuk contoh = np.array([5., 3., 1.]) = np.linalg.solve(B, C) print(A)
Institut Teknologi Bandung
# solusi # u/ debug
3
Bab 2
Praktikum 2.1
Interpolasi Newton
Lakukan langkah-langkah berikut: 1. Siapkan file “viscos0.csv” (7 baris) di bawah ini (dengan editor sederhana atau misalnya dengan spreadsheet Excel). Ini adalah relasi temperatur terhadap viskositas kinematik pada suatu fluida. 0.0 21.1 37.8 54.4 71.1 87.8 100.0
1.790 1.130 0.696 0.519 0.338 0.321 0.296
2. Untuk dapat melihat hasil yang lebih baik dalam grafik, plot data dan interpolasi akan dilakukan. Keseluruhan modul-modul yang diperlukan dipanggil sebagai berikut: import matplotlib.pyplot as plt import numpy as np from numpy import array,arange 3. Secara ringkas, ada beberapa urutan blok pada program yang harus Anda kerjakan (a) blok pemanggilan modul (lihat sebelumnya) (b) blok buka file, baca data, simpan data, dan tutup file (pelajari contoh sebelumnya) (c) blok perhitungan koefisien interpolan (d) blok penyiapan plot, contoh: 4
2.2. Regresi
Bab 2. Praktikum
plt.grid(True) plt.xlabel(’label sumbu x’) plt.ylabel(’label sumbu y’) plt.xlim(min(minx, 0.),1.1*maxx) plt.ylim(min(miny, 0.),1.1*maxy)
# # # # #
ya, pakai grid temperatur viskositas min & max u/ plot
(e) blok perhitungan interpolan Newton yang dibuat bersarang dalam loop untuk tabulasi atau plot menggunakan arange. Contoh: for x in arange(minx, maxx + h, h): # algoritma ASLINYA mulai dari baris di bawah ini p = a[n] # nilai awal polinom for k in range(1, n + 1): p = a[n - k] + (x - xData[n - k])*p # SAMPAI DI SINI y = p plt.plot(x,y,’b.’)
# black ’b’ plot w/ dot ’.’
(f) Jangan lupa, plot titik data dari file (kali ini pakai range), dan tunjukkan grafik for j in range(0,m): # red ’r’ plot w/ circle ’o’ plt.plot(xData[j], yData[j],’ro’,linewidth=5) plt.title(’Tulis judul grafik bila perlu’) plt.show()
2.2
Regresi
Untuk regresi (linear, polinomial), lakukan debugging blok demi blok pada program berikut, lalu berikan ulasan, manakah hasil regresi yang lebih baik? import matplotlib.pyplot as plt import numpy as np from numpy import array,arange from math import * from numpy.linalg import * ####################################################### # kinematic viscosity of water varies with temperatur T # T <--> xData # viscosity <--> yData ####################################################### fdata = open("viscos0.csv", "r") m = input(’Number of data, m = ’) xData = np.zeros(m) Institut Teknologi Bandung
5
Bab 2. Praktikum
2.2. Regresi
yData = np.zeros(m) sumx = 0.; sumx2 = 0.; sumx3 = 0.; sumx4 = 0. sumy = 0.; sumxy = 0.; sumx2y = 0. st = 0.; sr1 = 0.; sr2 = 0. # 2.2 11.098 # sample data for j in range(0,m): fdata.seek(0+11*j); xdt = fdata.read(3); xData[j] = xdt fdata.seek(4+11*j); ydt = fdata.read(6); yData[j] = ydt sumx += float(xdt) sumx2 += float(xdt)*float(xdt) sumx3 += float(xdt)*float(xdt)*float(xdt) sumx4 += float(xdt)*float(xdt)*float(xdt)*float(xdt) sumy += float(ydt) sumxy += float(xdt)*float(ydt) sumx2y += float(xdt)*float(xdt)*float(ydt) fdata.close() xmean = sumx/m ymean = sumy/m # # # # # # # #
print(m) print(sumx) print(sumx2) print(sumx3) print(sumx4) print(sumy) print(sumxy) print(sumx2y)
# linear regression a1 = (m*sumxy - sumx*sumy)/(m*sumx2 - sumx*sumx) a0 = ymean - a1*xmean # quadratic regression B = np.array([ [m, sumx, sumx2], [sumx, sumx2, sumx3], [sumx, sumx3, sumx4] ]) C A
= np.array([sumy, sumxy, sumx2y]) = np.linalg.solve(B, C)
print(’[A] = ’) print(A)
Institut Teknologi Bandung
6
2.2. Regresi
Bab 2. Praktikum
for j in range(0,m): st += (yData[j] sr1 += (yData[j] sr2 += (yData[j] (yData[j] syx1 syx2 r12 r22
= = = =
-
ymean)*(yData[j] - ymean) a0 - a1*xData[j])*(yData[j] - a0 - a1*xData[j]) A[0] - A[1]*xData[j] - A[2]*xData[j])*\ A[0] - A[1]*xData[j] - A[2]*xData[j])
sqrt(sr1/(m - 2)) sqrt(sr2/(m - 3)) (st - sr1)/st (st - sr2)/st
n = m - 1 minx maxx miny maxy
= = = =
min(xData[0],xData[n]) max(xData[0],xData[n]) min(yData[0],yData[n]) max(yData[0],yData[n])
print(’------------’) print(’ xi yi’) print(’------------’) for j in range(0,m): print(’%.1f %.3f’ % (xData[j], yData[j])) print(’------------’) print(’a[0] print(’a[1] print(’syx1 print(’syx2 print(’r1^2 print(’r2^2
= = = = = =
%+.5f’ %+.5f’ %+.5f’ %+.5f’ %+.5f’ %+.5f’
% % % % % %
a0) a1) syx1) syx2) r12) r22)
plt.grid(True) plt.xlabel(’label sumbu x’) plt.ylabel(’label sumbu x’) plt.xlim(min(minx, 0.),1.1*maxx) plt.ylim(min(miny, 0.),1.1*maxy) h = input(’plot interval, h = ’) der = input(’linear (1) / quadratic (2) ? der = ’) ################################# # ################################# if der==1: for x in arange(minx,maxx+h,h): y1 = a0 + a1*x plt.plot(x,y1,’b.’)
Institut Teknologi Bandung
7
Bab 2. Praktikum
2.2. Regresi
plt.title(’judul’) elif der==2: for x in arange(minx,maxx+h,h): y2 = A[0] + x*(A[1] + A[2]*x) plt.plot(x,y2,’b+’) plt.title(’judul’) for j in range(0,m): plt.plot(xData[j], yData[j],’ro’,linewidth=5) plt.show()
Institut Teknologi Bandung
8
