Pengolahan Objek

Operasi aljabar sederhana vektor numerik: Operatur yang biasa digunakan adalah +, -,, /, ^, %%, %o%

Pertama kita ingin menampilkan suatu vektor dengan elemen 2,6,9,5 seperti berikut dan diberi nama objek x1

x1 <- c(2,6,9,5)
x1
## [1] 2 6 9 5

objek lain yaitu vektor x2, dengan memanggil elemen 1:4 sehingga muncul hasil terdapat 4 elemen 1,2,3,4

x2 <- 1:4
x2
## [1] 1 2 3 4

Kemudian kita ingin menjumlahkan vektor x1 dengan elemen 1:2 menggunakan operator +. 1:2 merupakan vektor yang dengan elemennya yaitu 1 dan 2 sehingga apabila ingin dijumlahkan dengan vektor x1, secara otomatis mengalami recycle disebabkan panjang vektor tidak sama, sehingga diperoleh vektor x3 berikut

x3 <-x1 + 1:2
x3
## [1]  3  8 10  7

vektor x4 di bawah ini merupakan hasil jumlahan vektor x1 dan vektor 1:3 yang memiliki 3 elemen. Panjang vektor dari keduanya berbeda sehingga pada saat melakukan operasi bilangan secara otomatis akan di-recycle sehingga diperoelh hasil sebagai berikut

x4 <-x1 + 1:3
## Warning in x1 + 1:3: longer object length is not a multiple of shorter object
## length
x4
## [1]  3  8 12  6

kemudian kita ingin mengalikan setiap elemen yang ada pada vektor x1 dengan vektor x2 dengan menggunakan operator *. Sehingga diperoleh hasil perkalian seperti berikut

x5 <- x1*x2
x5
## [1]  2 12 27 20

sedangkan untuk operator %*% digunakan untuk perkalian vektor yang setara dengan x1'x1 dan menghasilkan skalar

x6 <- x1 %*% x1 #setara x'x
x6
##      [,1]
## [1,]  146

dan operatur %o$ digunakan untuk perkalian vektor yang setara dengan perkalian x1x1' sehingga hasil yang diperoleh adalah matrik berdimensi nxp

x7 <- x1 %o% x1 #setara xx'
x7
##      [,1] [,2] [,3] [,4]
## [1,]    4   12   18   10
## [2,]   12   36   54   30
## [3,]   18   54   81   45
## [4,]   10   30   45   25

Operasi dasar vektor karakter: nchar(???), paste(???), substr(???), substring(???)

Berikut ini merupakan vektor yang berisi karakter c(Institut Pertanian Bogor)

y1 <- c("Institut Pertanian Bogor")

y1
## [1] "Institut Pertanian Bogor"

kemudian kita ingin menghitung panjang karakter yang ada pada vektor y1 menggunakan fungsi nchar seperti di bawah ini

n1<-nchar(y1)
n1
## [1] 24

Contoh lain adalah kita memiliki vektor karakter y2 sebanyak 4 elemen

y2<-c("Adam","Pramesti","Fatih","Ririn")
y2
## [1] "Adam"     "Pramesti" "Fatih"    "Ririn"

kemudian dihitung karakter masing-masing elemen sehingga diperoleh

nchar(y2)
## [1] 4 8 5 5

Selanjutnya kita menjalankan objek y3 menggunakan fungsi substr dengan argumen(y1,20,24) artinya dengan menggunakan argumen y1 ("Institut Pertanian Bogor") kita ingin mengambil karakter dimulai dari karakter ke k=20 hingga karakter k=24. diperoleh hasil sebagai berikut

y3<-substr(y1,20,24)
y3
## [1] "Bogor"

kemudian kita ingin menjalankan fungsi substring dengan argumen (y1,15) artinya dengan menggunakan argumen y1 kita ingin mengambil karakter mulai dari karakter ke-k=15 hingga akhir, dan substring (y1,4,8) artinya dengan menggunakan argumen y1 kita ingin mengambil karakter mulai dari k=4 hingga batas k=8

substring(y1,15) #???nian Bogor???
## [1] "nian Bogor"
substring(y1,4,8) #???titut???
## [1] "titut"

atau dapat juga kita mengambil karakternya menggunakan [] dan diikuti urutan karakternya. Contohnya seperti berikut

y4<-substring(y2[2],2)
y4
## [1] "ramesti"

Operasi dasar matriks: t(???), %*%, %o%, diag(???), solve(???), eigen(???)

Misalkan kita memiliki matriks Z1 berukuran 2x3 dengan elemen 1:6

Z1 <- matrix(1:6,2,3)
Z1
##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6

kemudian kita memiliki matriks Z2berukuran 3x2 dengan elemen 1:6

Z2 <- matrix(1:6,3,2,byrow=T)
Z2
##      [,1] [,2]
## [1,]    1    2
## [2,]    3    4
## [3,]    5    6

kita juga mempunyai matriks Z3 dengan ukuran 2x2 dengan elemen 6:9

Z3 <- matrix(6:9,2,2)
Z3
##      [,1] [,2]
## [1,]    6    8
## [2,]    7    9

kemudian kita mengalikan matriks Z1 dan Z2 menggunakan operasi %*% sehingga diperoleh matriks hasil kali berukuran 2x2

Z4 <- Z1 %*% Z2
Z4
##      [,1] [,2]
## [1,]   35   44
## [2,]   44   56

kemudian, apabila kita mengalikan dua matrik menggunakan operasi *, maka kita mengalikan bilangan pada urutan baris dan kolom pada masing-masing matriks tersebut.

Z5 <- Z3 * Z4
Z5
##      [,1] [,2]
## [1,]  210  352
## [2,]  308  504

untuk menghitung invers menggunakan operasi solve()

INVZ <- solve(Z4) #invers
INVZ
##           [,1]      [,2]
## [1,]  2.333333 -1.833333
## [2,] -1.833333  1.458333

mengalikan suatu matriks dengan inversnya menghasilkan matriks identitas

INVZ %*% Z4 #identitas
##      [,1]         [,2]
## [1,]    1 2.842171e-14
## [2,]    0 1.000000e+00
h <- c(5,11)
h
## [1]  5 11
p <- solve(Z4,h) #solusi SPL Zp=h
p
## [1] -8.500  6.875

Menghitung Vektor eigen

e <- eigen(Z4) #eigen value & eigen vector dr Z4
e
## eigen() decomposition
## $values
## [1] 90.7354949  0.2645051
## 
## $vectors
##           [,1]       [,2]
## [1,] 0.6196295 -0.7848945
## [2,] 0.7848945  0.6196295

menghitung nilai eigen

e$values #akses eigen values
## [1] 90.7354949  0.2645051
e[[2]] #akses eigen vectors
##           [,1]       [,2]
## [1,] 0.6196295 -0.7848945
## [2,] 0.7848945  0.6196295

Struktur Kendali

Misalkan kita ingin melakukan perulangan dengan i ada pada sekuens 1:5 kemudian kita ingin menampilkan hasil pangkat dari i dimana expresinya yaitu i^2, maka diperoleh hasil sebagai berikut

for (i in 1:5) print(i^2)
## [1] 1
## [1] 4
## [1] 9
## [1] 16
## [1] 25

Misalkan kita memiliki objek x dengan sekuens 1:100

x<-1:100
x
##   [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
##  [19]  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36
##  [37]  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54
##  [55]  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72
##  [73]  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90
##  [91]  91  92  93  94  95  96  97  98  99 100

kemudian diambil sampel acak dari x tersebut sebanyak 5, diperoleh

acak1<-sample(x,5)
acak1
## [1] 29 41 57 53 79

kemudian kita melakukan perulangan objek acak1 menggunakan fungsi perulangan for sebanyak 10 pengulangan

for (i in 1:10)print(acak1)
## [1] 29 41 57 53 79
## [1] 29 41 57 53 79
## [1] 29 41 57 53 79
## [1] 29 41 57 53 79
## [1] 29 41 57 53 79
## [1] 29 41 57 53 79
## [1] 29 41 57 53 79
## [1] 29 41 57 53 79
## [1] 29 41 57 53 79
## [1] 29 41 57 53 79

fungsi pengulangan lainnya adalah while, dimana kita ingin melakukan pengulangan dengan kondisi i<=5 dan ekspresi i^2,

i<-1
while (i<=5) {
print(i^2)
i=i+1
}
## [1] 1
## [1] 4
## [1] 9
## [1] 16
## [1] 25

misalkan kita mempunyai y dengan sebaran seragam, kita ingin melakukan perulangan sebanyak 20 kali dengan for (kondisi i< 0.5) maka 100*i else i/100

y=runif(20)
for (i in y) {
if(i<0.5){
print(100*i)
} else print(i/100)
}
## [1] 0.009068921
## [1] 48.67142
## [1] 0.006709216
## [1] 41.56238
## [1] 21.9584
## [1] 22.26169
## [1] 33.16541
## [1] 34.6405
## [1] 20.72828
## [1] 0.005058022
## [1] 34.97792
## [1] 0.005328688
## [1] 0.00539941
## [1] 0.009482588
## [1] 19.76537
## [1] 10.24741
## [1] 0.007208348
## [1] 0.009049243
## [1] 9.278486
## [1] 0.007275262
z=0
while(z<=10) {
y=runif(20)
z=sum(y)
print(z)
}
## [1] 9.979289
## [1] 9.014451
## [1] 11.10092

Latihan 1

berikut adalah sintaks yang digunakan untuk melakukan perulangan menggunakan fungsi for

a<-0
for(i in 1:5)
{ b<-a+i
print(b)
a<-b
}
## [1] 1
## [1] 3
## [1] 6
## [1] 10
## [1] 15

Berikut adalah sintaks yang digunakan untuk melakukan perulangan menggunakan fungsi while

i<-1
z<-1
while(z<15)
{ y<-z+i
z<-y
print(z)
i<-i+1
}
## [1] 2
## [1] 4
## [1] 7
## [1] 11
## [1] 16
i<-1
m<-2
repeat
{ m<-m+i
print(m)
i<-i+1
if(m>15)
break
}
## [1] 3
## [1] 5
## [1] 8
## [1] 12
## [1] 17

Grafik

misalkan kita mempunyai x dengan sekuens 1:40 dan y mengikuti sebaran normal dengan n=40, rata-rata=5 dan simpangan baku sigma=1. kemudian kita membuat plot dari x dan y dengan menggunakan type="p" untuk memunculkan titik-titik datanya seperti di bawah ini

#dasar plot
x <- 1:40
y <- rnorm(40,5,1)
plot(x,y,type="p")

tipe lain yang dapat digunakan adalah "o" yaitu titik dan garis overlaid. seperti tampilan di bawah ini.

plot(x,y,type="o")

kemudian ada juga tipe "n" yang merupakan nothing, dapat dilihat seperti pada tampilan di bawah ini

plot(x,y,type="n")

kemudian kita juga dapat menampilkan plot dengan berbagai aksesoris misalnya ingin menambahkan label pada sumbu x dan y menggunakan fungsi xlab dan ylab, dan menambahkan judul plot dengan menggunakan fungsi main. Adapun untuk mengubah warna titik-titiknya dapat menggunakan fungsi col="". Kemudian menggunakan pch=16 untuk menambah plotting karakter dengan size dari pch itu sendiri menggunakan cex.Jika kita ingin membatasi sumbu dengan suatu batasan nilai, maka gunakan xlim untuk membatasi sumbu x dan ylim untuk membatasi sumbu y.

plot(x,y,type="p",xlab="Sumbu x",ylab="Sumbu y",main="Bilangan Acak Normal",col=2,pch=16)

plot(x,y,type="p",xlab="Sumbu x",ylab="Sumbu y",main="Bilangan Acak Normal",col=rainbow(40),pch=16,cex=2,xlim=c(0,50),ylim=c(2.5,7.5))

Kita juga dapat menambah amatan dan garis, panah, dan juga tulisan pada plot. Dapat dilihat pada tampilan di bawah ini

plot(x,y,type="p",xlab="Sumbu x",ylab="Sumbu y",main="Bilangan Acak Normal",col=rainbow(40),pch=16,cex=2,xlim=c(0,50),ylim=c(2.5,7.5))
#menambahkan amatan
x1 <- 41:50
y1 <- rnorm(10,5,1)
points(x1,y1,cex=2)
#menambahkan garis
x2 <- rep(40.5,20)
y2 <- seq(min(c(y,y1)),max(c(y,y1)),length=20)
lines(x2,y2,col=4,lwd=2,lty=2)
abline(h=mean(y),col="red",lwd=2.5)
abline(a=2,b=1/10,col="maroon3",lwd=2,lty=2)
#menambahkan tanda panah
arrows(x0=30,y0=3.5,x1=40,y1=mean(y)-.1,lwd=2)
#menambahkan tulisan
text(x=29,y=3.3,labels="Titik potong",cex=0.7)
text(x=3,y=7.3,labels="Data awal",cex=0.7)
text(x=46,y=7.3,labels="Data baru",cex=0.7)

plot(sin,-pi, 2*pi) #membuat grafik sin

plot(table(rpois(100,5)),type="h",col="red",
lwd=1,main="rpois(100,lambda=5)")     # membuat plot dari sebaran poisson dengan lambda=5

a1 <- 1:25
a2 <- rnorm(25,4,2)
plot(a1,a2,pch="w",xlab="x", ylab="y",main="W")   #membuat plot dengan data bangkitan acak sebaran normal, dimana titik-titik plot merupakan suatu katakter

plot(a1,a2,type="n",xlab="x", ylab="y", main="W")
text(a1,a2,labels=paste("w",1:25,sep=""),col=rainbow(25),cex=0.8)

x <- rchisq(100,df=4)    #membangkitkan data berdasarkan sebaran chi square
hist(x,freq=FALSE,ylim=c(0,0.2)) #membuat histogram dengan memberikan batas 0.2 pada sumbu y 
curve(dchisq(x,df=4),col=2,lty=2,lwd=2,add=TRUE)

Menggabungkan grafik dengan menggunakan par(mfrow=c(2,2))

par(mfrow=c(2,2))
plot(1:40,y,type="p",xlab="Sumbu x",ylab="Sumbu y",
main="Bilangan Acak Normal",col=2,pch=16)
plot(sin,-pi, 2*pi)
plot(table(rpois(100,5)),type="h",col="red",
lwd=1,main="rpois(100,lambda=5)")
plot(a1,a2,type="n",main="W")
text(a1,a2,labels=paste("w",1:25,sep=""),
col=rainbow(25),cex=0.8)

par(mfcol=c(2,2))
plot(1:40,y,type="p",xlab="Sumbu x",ylab="Sumbu y",
main="Bilangan Acak Normal",col=2,pch=16)
plot(sin,-pi, 2*pi)
plot(table(rpois(100,5)),type="h",col="red",lwd=1,
main="rpois(100,lambda=5)")
plot(a1,a2,type="n",main="W")
text(a1,a2,labels=paste("w",1:25,sep=""),
col=rainbow(25),cex=0.8)

windows()
yb <- rnorm(100,5,1)

menggabungkan grafik menggunakan split.screen

split.screen(c(2,2))
## [1] 1 2 3 4
screen(3)
boxplot(yb)
title("Boxplot Bilangan Acak Normal",cex.main=0.7)
screen(4)
xb <- 1:100
plot(xb,yb,type="l",lwd=2,col="blue")
title("Line Plot Bilangan Acak Normal",cex.main=0.7)
screen(2)
hist(yb,freq=FALSE,main=NULL,ylim=c(0,0.5))
x <- yb
curve(dnorm(x,5,1),col="red",lty=2,lwd=2,add=TRUE)
title("Histogram Bilangan Acak Normal",cex.main=0.7)
screen(1)
plot(xb,yb,pch=16,col=rainbow(100))
title("Scatter Plot Bilangan Acak Normal",cex.main=0.7)

close.screen()
## [1] 1 2 3 4

Pembangkitan Bilangan Acak

Berikut pembangkitan bilangan acak sebaran normal baku sebanyak 10 bilangan dengan nilai tengah miyu=0 dan simpangan baku sigma=1

#pembangkitan bilangan acak
x <- rnorm(10) #x~N(0,1)
x
##  [1]  0.29442951  1.24910197 -1.39668097 -2.43912602 -0.39595887 -0.63946352
##  [7]  0.13939974 -0.63994744  0.06098726  1.60755220

kemudian kita ingin membangkitkan bilangan acak sebaran normal sebanyak 10 bilangan dengan nilai tengah miyu=3 dan simpangan baku sigma=1 seperti di bawah ini

x1<- rnorm(10,3,2) #x1~N(3,sd=2)
x1
##  [1] 4.0199526 0.2922498 1.8736544 0.5766411 3.7038196 3.5613237 1.5681530
##  [8] 6.7531127 1.6728830 2.3761784

Selanjutnya, untuk membangkitkan bilangan acak yang mengikuti sebaran binomial, kita menggunakan fungsi rbinom. Berikut kita ingin membangkitkan bilangan acak mengikuti sebaran binomial dengan banyaknya amatan 10 bilangan dengan pecobaan sebanyak 1 kali dan memiliki nilai peluang=0.4

x2<- rbinom(10,1,0.4) #x2~bernoulli(0.4)
x2
##  [1] 1 1 0 1 0 1 0 0 1 1

untuk menghitung nilai peluang peubah acak yang mengikuti sebaran normal menggunakan fungsi pnorm. Di bawah ini kita ingin menghitung nilai peluang pada tingkat kepercayaan 90%, nilai taraf nyata alfa=10%. Kita peroleh Z(alfa/2)=z(0.95)=1.645.

#mencari nilai peluang peubah acak
p1 <- pnorm(1.645) #P(Z<1.645)=0.95
p1
## [1] 0.9500151

berikut nilai peluang peubah acak mengikuti sebaran normal dengan tingkat kepercayaan 95%, nilai taraf nyata alfa=5%

p2 <- pnorm(1.96) #P(Z<1.975)=0.975
p2
## [1] 0.9750021
p3 <- pnorm(-1.96)
p3
## [1] 0.0249979

Untuk menentukan peluang dari sebaran F, kita menggunakan fungsi pf, misalkan n=15, df1=10, df2=15, maka diperoleh nilai peluang berikut

p4 <- pf(15,df1=10,df2=15)
p4
## [1] 0.9999955

untuk mencari nilai kuantil dari peluang peubah acak menggunakan fungsi qnorm

#mencari nilai kuantil dari peluang peubah acak
q1 <- qnorm(0.975)
q1
## [1] 1.959964
q2 <- qnorm(0.95,2,1) #X~N(2,1), P(X<x)=0.95
q2
## [1] 3.644854
#mencari nilai density (fungsi kepekatan peluang) peubah acak
##plot density normal
a <- seq(-4,4,length=1000)
da <- dnorm(a)
plot(a,da)

??? Invers Transform Method ??? Acceptance-Rejection Method ??? Direct Transformation

latihan 1

Inverse Transform Method

#Eksponensial(lambda=3)
set.seed(5)
eks <- function(n,lambda){
U <- runif(n)
X <- -log(1-U)/lambda
return(X)
}
yy1 <- eks(1000,3) #inverse transform method
yy2 <- rexp(1000,rate=3) #fungsi bawaan R
par(mfrow=c(1,2))
hist(yy1,main="Exp dari Inverse Transform")
hist(yy2,main="Exp dari fungsi rexp")

Latihan 2

Acceptance-Rejection Method

ar <- function(n,x0,x1,f) {
xx <- seq(x0,x1,length=10000)
F <- max(f(xx))
terima <-0; 
iterasi<-0
hasil <- numeric(n)
while(terima<n) {
x <- runif(1,x0,x1)
y1<- runif(1,0,F)
y2<- f(x)
if(y1<=y2) {
terima <- terima+1
hasil[terima]<-x}
iterasi <- iterasi+1}
list(hasil=hasil,iterasi=iterasi)
}
set.seed(10)
f <- function(x) {3*x^2}
yyy <- ar(100,0,100,f
)
yyy; par(mfrow=c(1,1))
## $hasil
##   [1] 86.47212 77.51099 83.82877 77.07715 53.55970 86.13824 20.36477 79.79930
##   [9] 74.38394 34.43435 98.37322 69.35082 63.31153 88.80315 76.90405 64.83695
##  [17] 87.95432 93.60689 72.33519 76.20444 98.68082 87.60261 72.40640 81.40516
##  [25] 55.88949 89.00940 74.56896 84.80646 87.03302 82.23331 85.08123 77.09219
##  [33] 89.53595 58.03863 59.82260 92.35285 73.67755 68.98170 83.01572 92.93209
##  [41] 90.95163 53.47576 34.78601 87.59762 72.86815 87.49293 69.88356 83.12562
##  [49] 55.72238 66.47687 74.00502 98.06898 38.00746 75.53169 51.84889 88.79149
##  [57] 91.77773 80.84086 85.37441 42.32184 76.04306 34.05763 38.86568 47.74175
##  [65] 53.87605 94.85434 71.24685 90.81691 94.57656 77.16899 69.46655 53.68832
##  [73] 84.81593 82.42752 51.23742 31.52032 99.24487 93.27120 98.92809 62.83590
##  [81] 52.54605 88.10815 52.91748 57.65517 72.31807 87.61180 39.95670 89.86123
##  [89] 93.35217 78.59216 77.84128 69.55333 90.60413 99.16424 47.29846 97.70567
##  [97] 93.86110 99.59093 85.43663 83.09547
## 
## $iterasi
## [1] 322
hist(yyy$hasil, main="f(x)=3*x^2")

Direct Transformation

x <- runif(50)
y <- 2*x+3
y
##  [1] 4.595651 4.168212 3.865580 3.368448 3.685894 3.168929 4.232466 3.471525
##  [9] 3.816635 4.432023 3.654097 3.025882 4.017800 3.645034 4.035468 4.220792
## [17] 4.369285 4.777494 4.474976 4.379829 3.045420 3.458973 3.489782 4.082726
## [25] 3.896056 4.563540 3.730657 3.726081 3.958200 3.236057 4.076578 3.582958
## [33] 4.210074 3.726026 4.420964 3.973593 4.783335 3.698059 4.224458 3.630964
## [41] 3.646383 3.750261 4.467030 3.371470 3.469069 4.693028 4.894643 4.280374
## [49] 3.537450 4.790113
x=rbinom(10,2,0.3)
y=3*x-2
y
##  [1]  1  1  4 -2 -2 -2 -2  1  1 -2

Analisis Regresi

Misalnya kita mempunyai data bangkitan untuk variavbel X1 dari sebaran seragam, variabel X2 dari sebaran seragam, serta variabel respon Y seperti berikut

set.seed(123)
X1 <- runif(25,10,25)
X2 <- runif(25,90,200)
Y <- 10 + 2.3*X1 + 0.7*X2 + rnorm(25,0,9)
Y
##  [1] 145.2980 172.6297 157.2390 138.4862 151.0582 175.5627 181.0410 188.0289
##  [9] 184.1787 121.0427 171.9988 175.0267 135.4824 137.5021 113.9618 131.7872
## [17] 134.5388 117.9193 155.2350 151.5403 127.2692 134.2194 149.7767 157.8019
## [25] 183.9243
model1 <- lm(Y~X1)
summary(model1)
## 
## Call:
## lm(formula = Y ~ X1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.991 -17.738  -2.632  17.896  33.171 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  131.846     19.741   6.679 8.18e-07 ***
## X1             1.049      1.015   1.033    0.312    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.42 on 23 degrees of freedom
## Multiple R-squared:  0.04433,    Adjusted R-squared:  0.002775 
## F-statistic: 1.067 on 1 and 23 DF,  p-value: 0.3124

Hasil output di atas merupakan model regresi yang dapat dibentuk dari variabel respon Y dan variabel X1. Diperoleh hasil bahwa nilai p sebesar 0.312 lebih besar dari taraf nyata alfa 0.05. Sehingga dapat dikatakan bahwa variabel X1 tidak berpengaruh nyata terhadap variabel respon Y.

model2 <- lm(Y~X2)
summary(model2)
## 
## Call:
## lm(formula = Y ~ X2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.804  -4.850  -1.606  10.011  20.569 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 70.83913   13.86879   5.108 3.57e-05 ***
## X2           0.58213    0.09767   5.960 4.46e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.38 on 23 degrees of freedom
## Multiple R-squared:  0.607,  Adjusted R-squared:  0.5899 
## F-statistic: 35.52 on 1 and 23 DF,  p-value: 4.464e-06

Hasi di atas merupakan suatu model regresi yang dibentuk dari variabel X2 dan variabel respon Y. Diperoleh hasil bahwa nilai p sebesar 4.464e-06 lebih kecil dari taraf nyata alfa 0.05. Sehingga dapat dikatakan bahwa variabel X2 berpengaruh nyata terhadap variabel respon Y.

model3 <- lm(Y~X1+X2)
summary(model3)
## 
## Call:
## lm(formula = Y ~ X1 + X2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.8291  -4.5994  -0.4576   5.0602  20.5307 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.62881   13.40229  -0.047    0.963    
## X1           2.72951    0.40701   6.706 9.67e-07 ***
## X2           0.72459    0.06105  11.869 4.91e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.427 on 22 degrees of freedom
## Multiple R-squared:  0.8709, Adjusted R-squared:  0.8592 
## F-statistic: 74.21 on 2 and 22 DF,  p-value: 1.66e-10
model4 <- lm(Y~X1:X2)
summary(model4)
## 
## Call:
## lm(formula = Y ~ X1:X2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.929  -7.009  -1.430   2.856  37.374 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 80.645221  10.481428   7.694 8.32e-08 ***
## X1:X2        0.027491   0.003929   6.997 3.94e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.97 on 23 degrees of freedom
## Multiple R-squared:  0.6804, Adjusted R-squared:  0.6665 
## F-statistic: 48.96 on 1 and 23 DF,  p-value: 3.942e-07
model5 <- lm(Y~X1*X2)
summary(model5)
## 
## Call:
## lm(formula = Y ~ X1 * X2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.4425  -4.5808  -0.6702   4.6431  20.2409 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 11.625896  37.941991   0.306   0.7623  
## X1           2.063573   1.967526   1.049   0.3062  
## X2           0.635870   0.263687   2.411   0.0251 *
## X1:X2        0.004905   0.014165   0.346   0.7326  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.6 on 21 degrees of freedom
## Multiple R-squared:  0.8716, Adjusted R-squared:  0.8533 
## F-statistic: 47.53 on 3 and 21 DF,  p-value: 1.548e-09
R2 <- matrix(c(summary(model1)$r.squared,
summary(model1)$adj.r.squared,
summary(model2)$r.squared,
summary(model2)$adj.r.squared,
summary(model3)$r.squared,
summary(model3)$adj.r.squared,
summary(model4)$r.squared,
summary(model4)$adj.r.squared,
summary(model5)$r.squared,
summary(model5)$adj.r.squared), 5, byrow=T)
colnames(R2)<-c("R2","R2.adj"); R2*100
##             R2     R2.adj
## [1,]  4.432592  0.2774876
## [2,] 60.699195 58.9904640
## [3,] 87.090357 85.9167529
## [4,] 68.036265 66.6465370
## [5,] 87.163648 85.3298839
coef(model3)
## (Intercept)          X1          X2 
##  -0.6288095   2.7295126   0.7245911
confint(model3)
##                   2.5 %     97.5 %
## (Intercept) -28.4234482 27.1658292
## X1            1.8854322  3.5735929
## X2            0.5979779  0.8512044
cbind(coef(model3), confint(model3))
##                              2.5 %     97.5 %
## (Intercept) -0.6288095 -28.4234482 27.1658292
## X1           2.7295126   1.8854322  3.5735929
## X2           0.7245911   0.5979779  0.8512044
anova(model3)
## Analysis of Variance Table
## 
## Response: Y
##           Df  Sum Sq Mean Sq  F value    Pr(>F)    
## X1         1   536.4   536.4   7.5538   0.01173 *  
## X2         1 10002.4 10002.4 140.8614 4.911e-11 ***
## Residuals 22  1562.2    71.0                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow=c(2,2)); plot(model3)

Latihan Mandiri (Bahan Ajar Perkuliahan)

help(runif)
## starting httpd help server ... done
set.seed(1)
rnorm(10)
##  [1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078 -0.8204684
##  [7]  0.4874291  0.7383247  0.5757814 -0.3053884
rnorm(10,mean=5,sd=1)
##  [1] 6.511781 5.389843 4.378759 2.785300 6.124931 4.955066 4.983810 5.943836
##  [9] 5.821221 5.593901
a<-runif(10)
a
##  [1] 0.8209463 0.6470602 0.7829328 0.5530363 0.5297196 0.7893562 0.0233312
##  [8] 0.4772301 0.7323137 0.6927316

Bangkitkan X~N(10,1)

langkah pertama: bangkitkan vektor r langkah kedua :Gunakan fungsi cdf(fungsi sebaran kumulatif) dari sebaran yang diinginkan untuk mendapatkan X~F^(-1) (r) inverse dari r

r<-runif(100)
x<-qnorm(r,mean=10,sd=1) #x merupakan inverse (qnorm) atau fungsi kuantil dari r dengan mean=10 dan sd 1
hist(x,prob=TRUE)
sbx<-seq(6,14,.01) #sbx merupakan sebaran dari x (ini yang dimaksud dengan suatu sebaran yang diinginkan) F(x)
f.massa.p<-dnorm(sbx,mean=10,sd=1) #dnorm itu adalah fungsi massa peluang atau fungsi kepekatan peluang dilambangkan f(x)
lines(sbx,f.massa.p)

inverse transform method

peubah acak kontinyu bangkitkan X~f(x)=3x^2, 0

peubah acak diskrit bangkitkan X~Binom(5,0.5) #n=5 p=0.5

set.seed(1)
u=runif(100)  #sebaran seragam berukuran n
x= qbinom(u,size=5, prob=.5)
barplot(table(x))

u=runif(100)
x=ifelse(u<=0.5,0,
         ifelse(u<=0.8,1,2))
barplot(table(x))

Acceptance-rejection method

n=1000
j=k=0
y=numeric(n)
while(k<n) {
  u=runif(1)
  j=j+1
  x= runif(1)
  if (3*x*(1-x)>u){
    k=k+1 #terima x
    y[k]=x
  }
}
hist(y,prob=TRUE)
sbx=seq(0,1,.01)
lines(sbx,dbeta(sbx,2,2))

j
## [1] 2031

Direct Transformation

misal: bangkitkan bilangan acak ~ Chisq(df=1)

y= rnorm(1000)
x=y^2
hist(x, prob=T)
sbx<- seq(0,13,0.01)
lines(sbx,dchisq(sbx,df=1), col="red")

Multivariate normal distribution

rmvn.eigen <- function (n , mu , Sigma ) {
 d <- length ( mu )
 ev <- eigen ( Sigma , symmetric = TRUE )
 lambda <- ev$ values
 V <- ev$ vectors
 R <- V %*% diag ( sqrt ( lambda ) ) %*% t( V )
 Z <- matrix ( rnorm ( n*d ) , nrow = n , ncol = d )
 X <- Z %*% R + matrix ( mu , n , d , byrow = T )
 X
 }
 mu <- c(0 , 0)
 Sigma <- matrix (c(1 , .9 , .9 , 1) , nrow =2 , ncol =2)
 X <- rmvn.eigen (1000 , mu , Sigma )
 cor(X)
##           [,1]      [,2]
## [1,] 1.0000000 0.8969975
## [2,] 0.8969975 1.0000000
cov(X)
##           [,1]      [,2]
## [1,] 1.0139683 0.8926498
## [2,] 0.8926498 0.9766870
var(X)
##           [,1]      [,2]
## [1,] 1.0139683 0.8926498
## [2,] 0.8926498 0.9766870
library(MASS)
X1=mvrnorm(1000,mu,Sigma)
var(X1)
##           [,1]      [,2]
## [1,] 1.0181576 0.8980061
## [2,] 0.8980061 0.9888291
X2=mvrnorm(1000,mu,Sigma, empirical=TRUE)
var(X2)
##      [,1] [,2]
## [1,]  1.0  0.9
## [2,]  0.9  1.0

membangkitkan bilangan acak untuk model

b0 <- 1; b1 <- 1
b0hat <- NULL ; b1hat <- NULL
for ( i in 1:100) {
eps <- rnorm (10)
X <- runif (10 ,5 ,10)
Y <- b0 + b1*X + eps
obj <- lm( Y~X )
b0hat <- c( b0hat , obj$ coefficients [1])
b1hat <- c( b1hat , obj$ coefficients [2])
 }
hasil <- matrix (c( mean ( b0hat ) , sd( b0hat ) , mean ( b1hat ) ,sd(
b1hat ) ) ,nrow =2 , ncol =2)
rownames ( hasil ) <- c(" mean ","sd")
colnames ( hasil ) <- c("b0", "b1")
hasil
##               b0        b1
##  mean  0.9931473 1.0047069
## sd     1.7392996 0.2181047
library ( MASS )
b0 <- 1; b1 <- 1; b2 <- 1
b0hat <- NULL ; b1hat <- NULL ; b2hat <- NULL
Sigma <- matrix (c (1 ,0.9 ,0.9 ,1) ,nrow =2 , ncol =2)
mu <- c(1 ,1)
for ( i in 1:100) {
eps <- rnorm (10)
X <- mvrnorm (10 , mu , Sigma )
Y <- b0 + b1*X [ ,1] * b2*X [ ,2] + eps
obj <- lm( Y~X )
b0hat <- c( b0hat , obj$ coefficients [1])
b1hat <- c( b1hat , obj$ coefficients [2])
b2hat <- c( b2hat , obj$ coefficients [3])
 }
hasil <- matrix (c( mean ( b0hat ) , sd( b0hat ) , mean ( b1hat ) , sd(b1hat ) , mean ( b2hat ) , sd( b2hat ) ) , nrow =2 , ncol =3)
rownames ( hasil ) <- c(" mean ","sd")
colnames ( hasil ) <- c("b0", "b1", "b2")
hasil
##               b0        b1       b2
##  mean  0.6394752 0.9253045 1.135690
## sd     0.9045239 1.2860522 1.233582