Bab 6 Fitting functions to data

Jika anda memiliki ide tentang bentuk fungsi pada sebua model dan anda perlu memilih parameter yang akan membuat fungsi model cocok untuk observasi. Proses pemilihan parameter untuk mencocokkan pengamatan disebut model fitting.

Sebagai contoh, data dalam file “utilities.csv” mencatat suhu rata-rata setiap bulan(dalam parameter fahrenheit) serta penggunaan gas alam bulanan(dalam kaki rubik,ccf)

library(mosaicCalc)
## Warning: package 'mosaicCalc' was built under R version 4.2.2
## Loading required package: mosaic
## Warning: package 'mosaic' was built under R version 4.2.2
## Registered S3 method overwritten by 'mosaic':
##   method                           from   
##   fortify.SpatialPolygonsDataFrame ggplot2
## 
## The 'mosaic' package masks several functions from core packages in order to add 
## additional features.  The original behavior of these functions should not be affected by this.
## 
## Attaching package: 'mosaic'
## The following objects are masked from 'package:dplyr':
## 
##     count, do, tally
## The following object is masked from 'package:Matrix':
## 
##     mean
## The following object is masked from 'package:ggplot2':
## 
##     stat
## The following objects are masked from 'package:stats':
## 
##     binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
##     quantile, sd, t.test, var
## The following objects are masked from 'package:base':
## 
##     max, mean, min, prod, range, sample, sum
## Loading required package: mosaicCore
## Warning: package 'mosaicCore' was built under R version 4.2.2
## 
## Attaching package: 'mosaicCore'
## The following objects are masked from 'package:dplyr':
## 
##     count, tally
## 
## Attaching package: 'mosaicCalc'
## The following object is masked from 'package:stats':
## 
##     D
Utils <- read.csv("http://www.mosaic-web.org/go/datasets/utilities.csv")
gf_point(ccf ~ temp, data = Utils) %>%
  gf_labs(y = "Natural gas usage (ccf/month)", 
          x = "Average outdoor temperature (F)")

Banyak jenis fungsi yang berbeda dapat digunakan untuk mewakili data ini. Salah satu yang paling sederhana dan paling umum digunakan dalam pemodelan adalah fungsi garis lurus f(x)=Ax+B . Dalam fungsi f(x) , variabel x singkatan dari input, sedangkan A dan B adalah parameter. Penting untuk diingat apa nama input dan output saat menyesuaikan model dengan data – Anda perlu mengatur agar namanya cocok dengan data yang sesuai.

Dengan data utilitas, masukannya adalah suhu, suhu. Keluaran yang akan dimodelkan adalah ccf. Untuk menyesuaikan fungsi model dengan data, tuliskan rumus dengan nama input, parameter, dan output yang sesuai di tempat yang tepat:

f <- fitModel(ccf ~ A * temp + B, data = Utils)

Keluaran dari fitModel() adalah fungsi dengan bentuk matematika yang sama seperti yang Anda tentukan di argumen pertama (di sini, ccf ~ A * temp + B) dengan nilai numerik spesifik yang diberikan ke parameter untuk membuat fungsi paling cocok dengan data. Bagaimana cara fitModel()mengetahui besaran mana dalam bentuk matematika yang merupakan variabel dan mana yang merupakan parameter? Apa pun yang terkandung dalam data yang digunakan untuk pemasangan adalah variabel (di sini temp); hal-hal lain (di sini, A dan B) adalah parameter.

gf_point(ccf ~ temp, data = Utils) %>%
  slice_plot(f(temp) ~ temp)

Anda dapat menambahkan fungsi lain ke dalam campuran dengan mudah. Misalnya, Anda mungkin berpikir itu sqrt(temp)berhasil di sana.

f2 <- fitModel(
  ccf ~ A * temp + B + C *sqrt(temp),
  data = Utils)
gf_point(
  ccf ~ temp, data = Utils) %>%
  slice_plot(f2(temp) ~ temp)

Contoh ini hanya melibatkan satu variabel masukan. Sepanjang ilmu alam dan sosial, teknik yang sangat penting dan banyak digunakan adalah menggunakan banyak variabel dalam sebuah proyeksi. Sebagai ilustrasi, lihat data “used-hondas.csv”harga mobil Honda bekas.

Hondas <- read.csv("http://www.mosaic-web.org/go/datasets/used-hondas.csv")
head(Hondas)
##   Price Year Mileage Location Color Age
## 1 20746 2006   18394  St.Paul  Grey   1
## 2 19787 2007       8  St.Paul Black   0
## 3 17987 2005   39998  St.Paul  Grey   2
## 4 17588 2004   35882  St.Paul Black   3
## 5 16987 2004   25306  St.Paul  Grey   3
## 6 16987 2005   33399  St.Paul Black   2

Seperti yang Anda lihat, kumpulan data menyertakan variabel Price, Age, dan Mileage. Tampaknya masuk akal untuk berpikir bahwa harga akan bergantung pada jarak tempuh dan usia mobil. Inilah model yang sangat sederhana yang menggunakan kedua variabel:

carPrice1 <- fitModel(
  Price ~ A + B * Age + C * Mileage, data = Hondas
)

Anda dapat memplot fungsi yang sesuai:

contour_plot(
  carPrice1(Age = age, Mileage = miles) ~ age + miles,
  domain(age=2:8, miles=range(0, 60000)))

6.0.1 Latihan

6.0.1.1 Latihan 1

Teks bagian ini menjelaskan model carPrice1()dengan Umur dan Jarak Tempuh sebagai jumlah masukan dan harga (dalam USD) sebagai keluaran. Klaim dibuat bahwa harga dapat dilihat sebagai fungsi dari Agedan Mileage. Mari kita buat grafik itu lagi.

contour_plot(
  carPrice1(Age = age, Mileage = miles) ~ age + miles,
  domain(age = range(0, 8), miles = range(0, 60000)))

  1. Apa yang dikatakan kontur vertikal tentang harga sebagai fungsi dari Agedan Mileage? a.Harga sangat tergantung pada kedua variabel. b.Harga tergantung Agetapi tidak Mileage. c.Harga tergantung Mileagetapi tidak Age. d.Harga tidak banyak bergantung pada salah satu variabel.

JAWABAN: Setiap kontur sesuai dengan harga yang berbeda. Saat Anda melacak secara horizontal dengan Age, Anda berpindah dari satu kontur ke kontur lainnya. Namun saat Anda melacak secara vertikal dengan Mileage, Anda tidak melewati kontur. Artinya harga tidak bergantung pada Mileage, karena perubahan Mileagetidak menyebabkan perubahan harga. Tapi harga tidak berubah dengan Age.

  1. Grafik fungsi yang sama yang ditampilkan di badan teks memiliki kontur yang miring ke bawah dari kiri ke kanan. Apa yang dikatakan tentang harga sebagai fungsi dari Age dan Mileage? a.Harga sangat tergantung pada kedua variabel. b.Harga tergantung Agetapi tidak Mileage. c.Harga tergantung Mileagetapi tidak Age. d.Harga tidak banyak bergantung pada salah satu variabel.

JAWABAN: Saat Anda menjiplak secara horizontal, dengan Age, Anda berpindah dari kontur ke kontur: harga berubah. Jadi harga tergantung Age. Hal yang sama berlaku ketika Anda menjiplak secara vertikal, dengan Mileage. Jadi harga juga tergantung Mileage.

  1. Fungsi yang sama digambarkan baik di badan teks maupun di latihan ini. Tapi grafiknya sangat berbeda! Jelaskan mengapa ada perbedaan dan sebutkan grafik mana yang benar.

JAWABAN: Perhatikan tanda centang pada sumbu. Dalam grafik di badan teks, Ageberjalan dari dua hingga delapan tahun. Namun dalam grafik latihan, Agehanya berjalan dari nol hingga satu tahun. Demikian pula, grafik di badan teks Mileageberjalan dari 0 hingga 60.000 mil, tetapi dalam grafik latihan, Mileageberjalan dari 0 hingga 1. Kedua grafik menunjukkan fungsi yang sama, jadi keduanya “benar”. Tapi grafik latihan itu menyesatkan secara visual. Tidak mengherankan jika harga tidak banyak berubah dari 0 mil menjadi 1 mil, tetapi harganya (agak) berubah dari 0 tahun menjadi 1 tahun

6.1 Kurva dan model linier

Kata linier dalam “model linier” mengacu pada “kombinasi linier”, bukan “garis lurus”. Seperti yang akan Anda lihat, Anda dapat membuat kurva rumit dengan mengambil kombinasi fungsi linier, dan menggunakan operasi proyeksi aljabar linier untuk mencocokkan kurva ini sedekat mungkin dengan data. Proses pencocokan itu disebut “pas”.

Sebagai contoh, data dalam file “utilities.csv”merekam suhu rata-rata setiap bulan (dalam derajat Fahrenheit) serta penggunaan gas alam bulanan (dalam kaki kubik, ccf). Ada, seperti yang Anda duga, hubungan yang kuat antara keduanya.

Utilities = read.csv("http://www.mosaic-web.org/go/datasets/utilities.csv")
gf_point(ccf ~ temp, data = Utilities)

Banyak jenis fungsi yang berbeda dapat digunakan untuk mewakili data ini. Salah satu yang paling sederhana dan paling umum digunakan dalam pemodelan adalah fungsi garis lurus. Dalam hal aljabar linier, ini adalah kombinasi linier dari fungsi f1(T)=1 dan f2(T)=T . Secara konvensional, tentu saja, fungsi garis lurus ditulis f(T)=b+mT . (Mungkin Anda lebih suka menulisnya seperti ini: f(x)=mx+b . Hal yang sama.) Notasi konvensional ini hanyalah penamaan skalar sebagai m dan b yang akan berpartisipasi dalam kombinasi linier. Untuk menemukan skalar numerik yang paling cocok dengan data — untuk “menyesuaikan fungsi” dengan data — dapat dilakukan dengan project( )operator aljabar linier.

project(ccf ~ temp + 1, data = Utilities)
## (Intercept)        temp 
##  253.098208   -3.464251

Operator project( )memberikan nilai skalar. Fungsi pemasangan terbaik itu sendiri dibangun dengan menggunakan nilai skalar ini untuk menggabungkan fungsi yang terlibat.

model_fun = makeFun( 253.098 - 3.464*temp ~ temp)
gf_point(ccf ~ temp, data=Utils) %>%
  slice_plot(model_fun(temp) ~ temp)

Anda dapat menambahkan fungsi lain ke dalam campuran dengan mudah. Misalnya, Anda mungkin berpikir itu sqrt(T)berhasil di sana. Cobalah!

project(ccf ~ temp  + sqrt(temp) + 1, data = Utils)
## (Intercept)        temp  sqrt(temp) 
##  447.029273    1.377666  -63.208025
mod2 <- makeFun(447.03 + 1.378*temp - 63.21*sqrt(temp) ~ temp)
gf_point(ccf ~ temp, data=Utils) %>% # the data
  slice_plot(mod2(temp) ~ temp) %>%
  gf_labs(x = "Temperature (F)", 
          y = "Natural gas used (ccf)")

Operator project( )mengambil serangkaian vektor. Saat menyesuaikan fungsi ke data, vektor ini berasal dari kumpulan data sehingga perintah harus mengacu pada nama besaran seperti yang muncul di kumpulan data, misalnya, ccf atau temp. Anda diperbolehkan melakukan operasi pada besaran tersebut, misalnya pada contoh sqrtdi atas, untuk membuat vektor baru. The ~digunakan untuk memisahkan vektor “target” dari kumpulan satu atau lebih vektor tempat proyeksi dibuat. Dalam notasi matematika tradisional, operasi ini akan ditulis sebagai persamaan yang melibatkan matriks A terdiri dari sekumpulan vektor (→v1,→v2,…,→vp)=A , vektor sasaran →b , dan himpunan koefisien yang tidak diketahui →x . Persamaan yang menghubungkan jumlah ini ditulis A⋅→x≈→b . Dalam notasi ini, proses “solving” untuk →x implisit. Notasi komputer mengatur ulang ini menjadi

x=→b∼→v1+→v2+…+→vp.

Setelah Anda menyelesaikan proyeksi dan menemukan koefisien, Anda dapat membuat fungsi matematika yang sesuai dengan menggunakan koefisien dalam ekspresi matematika untuk membuat fungsi. Seperti semua fungsi, nama yang Anda gunakan untuk argumen adalah masalah pilihan pribadi, meskipun masuk akal untuk menggunakan nama yang mengingatkan Anda tentang apa yang diwakili oleh fungsi tersebut.

Sepanjang ilmu alam dan sosial, teknik yang sangat penting dan banyak digunakan adalah menggunakan banyak variabel dalam sebuah proyeksi. Sebagai ilustrasi, lihat data “used-hondas.csv”harga mobil Honda bekas.

Hondas = read.csv("http://www.mosaic-web.org/go/datasets/used-hondas.csv")
head(Hondas)
##   Price Year Mileage Location Color Age
## 1 20746 2006   18394  St.Paul  Grey   1
## 2 19787 2007       8  St.Paul Black   0
## 3 17987 2005   39998  St.Paul  Grey   2
## 4 17588 2004   35882  St.Paul Black   3
## 5 16987 2004   25306  St.Paul  Grey   3
## 6 16987 2005   33399  St.Paul Black   2

kumpulan data menyertakan variabel Price, Age, dan Mileage. Tampaknya masuk akal untuk berpikir bahwa harga akan bergantung pada jarak tempuh dan usia mobil. Inilah model yang sangat sederhana yang menggunakan kedua variabel:

project(Price ~ Age + Mileage + 1, data = Hondas)
##   (Intercept)           Age       Mileage 
##  2.133049e+04 -5.382931e+02 -7.668922e-02

Anda dapat memplotnya sebagai fungsi matematika

car_price <- makeFun(21330-5.383e2*age-7.669e-2*miles ~ age & miles)
contour_plot(car_price(age, miles) ~ age + miles,
  domain(age=range(2, 8), miles=range(0, 60000))) %>%
  gf_labs(title = "Miles per gallon")

Model yang agak lebih canggih mungkin menyertakan apa yang disebut “interaksi” antara usia dan jarak tempuh, menyadari bahwa pengaruh usia mungkin berbeda tergantung pada jarak tempuh

project(Price ~ Age + Mileage + Age*Mileage + 1, data = Hondas)
##   (Intercept)           Age       Mileage   Age:Mileage 
##  2.213744e+04 -7.494928e+02 -9.413962e-02  3.450033e-03
car_price2 <- makeFun(22137 - 7.495e2*age - 9.414e-2*miles +
                         3.450e-3*age*miles ~ age & miles)
contour_plot(
  car_price2(Age, Mileage) ~ Age + Mileage,  
  domain(Age = range(0, 10), Mileage = range(0, 100000))) %>%
  gf_labs(title = "Price of car (USD)")

6.1.1 Latihan

6.1.1.1 Latihan 1: Pemasangan polinimial

Sebagian besar mahasiswa mengambil kursus aljabar yang mencakup banyak tentang polinomial, dan polinomial sangat sering digunakan dalam pemodelan. (Mungkin, mereka digunakan lebih sering daripada yang seharusnya. Dan guru aljabar mungkin kecewa mendengar bahwa model polinomial yang paling penting adalah model orde rendah, misalnya, f(x,Y)=A+bx+cY+dxY daripada menjadi kubik atau kuartik, dll.) Menyesuaikan polinomial dengan data adalah masalah aljabar linier: menyusun vektor yang sesuai untuk mewakili berbagai kekuatan. Misalnya, inilah cara menyesuaikan model kuadrat dengan variabel ccfversus dalam file data:temp “utilities.csv”

Utilities = read.csv("http://www.mosaic-web.org/go/datasets/utilities.csv")
project(ccf ~ 1 + temp + I(temp^2), data = Utilities)
##  (Intercept)         temp    I(temp^2) 
## 317.58743630  -6.85301947   0.03609138

Anda mungkin bertanya-tanya, I( )untuk apa? Ternyata ada notasi yang berbeda untuk statistik dan matematika, dan ^ memiliki arti yang agak berbeda dalam rumus R daripada eksponensial sederhana. Memberitahu I( )software untuk mengambil eksponensial secara harfiah dalam arti matematis. Koefisien memberi tahu kita bahwa model kuadrat yang paling pas dari ccf versus temp adalah:

ccfQuad <- makeFun(317.587 - 6.853*T + 0.0361*T^2 ~ T)
gf_point(ccf ~ temp, data = Utilities) %>%
  slice_plot(ccfQuad(temp) ~ temp)

Untuk mencari nilai model ini pada temperatur tertentu, evaluasi saja fungsinya. (Dan perhatikan bahwa ccfQuad( )didefinisikan dengan variabel input T.)

ccfQuad(T=72)
## [1] 11.3134
  1. Sesuaikan polinomial orde-3 versus dengan data utilitas. Berapa nilai model ini untuk suhu 32 derajat? {87.103.128, 142 .143.168.184} MENJAWAB:
project(ccf ~ 1 + temp + I(temp^2) + I(temp^3), data = Utils)
##   (Intercept)          temp     I(temp^2)     I(temp^3) 
##  2.550709e+02 -1.427408e+00 -9.643482e-02  9.609511e-04
ccfCubic <- 
  makeFun(2.551e2 - 1.427*T -
          9.643e-2*T^2 + 9.6095e-4*T^3 ~ T)
gf_point(ccf ~ temp, data = Utils) %>%
  slice_plot(ccfCubic(temp) ~ temp) 

ccfCubic(32)
## [1] 142.1801

6.2 Fungsi dengan parameter nonlinier

Teknik aljabar linier dapat digunakan untuk mencari kombinasi linier terbaik dari suatu himpunan fungsi. Namun, seringkali, ada parameter dalam fungsi yang muncul secara nonlinier. Menemukan parameter nonlinier tidak dapat dilakukan secara langsung menggunakan aljabar linier, meskipun metode aljabar linier memang membantu menyederhanakan situasi. Untungnya, gagasan bahwa jarak antar fungsi dapat diukur berfungsi dengan baik ketika ada parameter nonlinear yang terlibat. Jadi kita akan terus menggunakan “jumlah residu kuadrat” saat mengevaluasi seberapa dekat perkiraan fungsi dengan sekumpulan data.

6.3 Fungsi eksponensial

Families <- read.csv("http://www.mosaic-web.org/go/datasets/Income-Housing.csv")
gf_point(TwoVehicles ~ Income, data = Families)

Pola data menunjukkan “pembusukan” eksponensial terhadap hampir 100% keluarga yang memiliki dua kendaraan. Bentuk matematis dari fungsi eksponensial ini adalah A dan xp(kY)+C . A dan C adalah parameter linier yang tidak diketahui. k adalah parameter nonlinear yang tidak diketahui – ini akan menjadi negatif untuk peluruhan eksponensial. Aljabar linier memungkinkan kita menemukan parameter linier terbaik A dan C untuk mencocokkan data. Tapi bagaimana menemukan k ?

Anda dapat melihat dari datanya sendiri bahwa “waktu paruh” adalah sekitar $25.000. Parameter k adalah sesuai dengan paruh, itu di(0.5)/setengah hidup , jadi inilah tebakan yang bagus untuk k adalah di(0.5)/25000 , itu adalah

kguess <- log(0.5) / 25000
kguess
## [1] -2.772589e-05

Dimulai dengan tebakan tersebut, Anda dapat menemukan nilai terbaik dari parameter linier A dan C melalui teknik aljabar linier:

project( TwoVehicles ~ 1 + exp(Income*kguess), data = Families)
##          (Intercept) exp(Income * kguess) 
##             110.4263            -101.5666

Pastikan bahwa Anda benar-benar memahami arti dari pernyataan di atas. Itu TIDAK berarti itu TwoVehicles jumlahnya 1+exp−Penghasilan×tebak . Sebaliknya, itu berarti Anda sedang mencari kombinasi linier dari kedua fungsi tersebut 1 dan exp−Penghasilan×tebak yang cocok TwoVehiclessedekat mungkin. Nilai yang dikembalikan oleh memberi tahu Anda akan seperti apa kombinasi ini: berapa banyak 1 dan berapa banyak exp−Penghasilan×tebak untuk menambahkan bersama-sama untuk mendekati TwoVehicles.

Anda dapat membuat fungsi yang merupakan kombinasi linier terbaik dengan menambahkan kedua fungsi secara eksplisit:

f <- makeFun( 110.43 - 101.57*exp(Income * k) ~ Income, k = kguess)
gf_point(TwoVehicles ~ Income, data = Families) %>%
  slice_plot(f(Income) ~ Income)

Grafik berjalan sangat dekat dengan titik data. Tapi Anda juga bisa melihat nilai numerik dari fungsi untuk setiap pendapatan:

f(Income = 10000)
## [1] 33.45433
f(Income = 50000)
## [1] 85.0375

Sangat informatif untuk melihat nilai fungsi untuk Incomelevel tertentu dalam data yang digunakan untuk pemasangan, yaitu data frame Families:

Results <- Families %>% 
  dplyr::select(Income, TwoVehicles) %>%
  mutate(model_val = f(Income = Income),
         resids = TwoVehicles - model_val)
Results
##   Income TwoVehicles model_val     resids
## 1   3914        17.3  19.30528 -2.0052822
## 2  10817        34.3  35.17839 -0.8783904
## 3  21097        56.4  53.84097  2.5590313
## 4  34548        75.3  71.45680  3.8432013
## 5  51941        86.6  86.36790  0.2320981
## 6  72079        92.9  96.66273 -3.7627306

Residual adalah perbedaan antara nilai model ini dan nilai sebenarnya dari TwoVehicleskumpulan data.

Kolom residsmemberikan sisa untuk setiap baris. Tapi Anda juga bisa memikirkan residskolom sebagai vektor . Ingatlah bahwa panjang kuadrat vektor adalah jumlah residu kuadrat

sum(Results$resids^2)
## [1] 40.32358

Panjang persegi residsvektor ini merupakan cara penting untuk mengukur seberapa cocok model dengan data.

6.5 Mengoptimalkan tebakan

Perlu diingat bahwa jumlah residu kuadrat adalah fungsi dari k. Nilai di atas hanya untuk tebakan khusus kita $k = $ kguess. Daripada hanya menggunakan satu tebakan untuk k, Anda dapat melihat berbagai kemungkinan. Untuk melihat semuanya pada saat yang sama, mari kita gambarkan jumlah residu kuadrat sebagai fungsi dari k. Kita akan melakukan ini dengan membuat fungsi yang menghitung jumlah sisa kuadrat untuk nilai tertentu dari k.

sum_square_resids <- Vectorize(function(k) {
  sum((Families$TwoVehicles - f(Income=Families$Income, k)) ^ 2)
})
slice_plot(
   sum_square_resids(k) ~ k, 
   domain(k = range(log(0.5)/40000,log(0.5)/20000)))

Ini adalah perintah komputer yang agak rumit, tetapi grafiknya langsung. Anda dapat melihat bahwa nilai “terbaik” dari k, yaitu nilai dari kyang membuat jumlah residu kuadrat sekecil mungkin, sudah dekat k=−2.8×10-5 tidak terlalu jauh dari tebakan awal, seperti yang terjadi. (Itu karena waktu paruh sangat mudah diperkirakan.)

Untuk melanjutkan penjelajahan Anda dalam penyesuaian kurva nonlinier, Anda akan menggunakan fungsi tujuan khusus yang melakukan sebagian besar pekerjaan untuk Anda sambil memungkinkan Anda mencoba berbagai nilaikkdengan menggerakkan slider.