Grafik dan Bilangan Acak
Grafik/Plot
Bentuk Plot
#dasar plot
x <- 1:40
y <- rnorm(40,5,1)
plot(x,y, type="p") # "p" for points
plot(x,y, type="l") # “l" for lines
plot(x,y, type="b") # “b" for both
plot(x,y, type="c") # "c" for for the lines part alone of “b"
plot(x,y, type="o") # “o" for overplotted
plot(x,y, type="n") # “n" for no plotting
plot(x,y, type="h") # “h" for histogram like
plot(x,y, type="s") # “s" for stair steps
plot(x,y, type="S") # “S" for other steps
plot(x,y, type="n") # “n" for no plotting
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 =2, pch=16)
#menambahkan amatan
x1 <- 41:50
y1 <- rnorm(10,5,1)
points(x1,y1,cex=2)
#menambahkan amatan
x1 <- 41:50
y1 <- rnorm(10,5,1)
points(x1,y1,cex=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)
Latihan 1
Could you explain what are these programs do for ?
• plot(sin, -pi, 2*pi)
• plot(table(rpois(100,5)),type=“h”,col =“red”,lwd=1,main=“rpois(100,lambda=5)”)
plot(sin,-pi, 2*pi)
> Script plot(sin, -pi, 2*pi) membuat garifk antara nilai x yang merupakan nilai dari -phi dengan sin(-phi) dari rentang axis X 0 sampai dengan 2 kali phi
plot(table(rpois(100,5)),type="h",col ="red",lwd=1,main="rpois(100,lambda=5)")
Script
plot(table(rpois(100,5)),type="h",col ="red",lwd=1,main="rpois(100,lambda=5)")membuat bar plot dari sebaran poisson dengan lambda = 5
Latihan 2
• Create some programs to make the graph below
a1 <- 1:25
a2 <- rnorm(25,4,2)
plot(a1, a2,pch="w",main="W")
Latihan 3
plot(a1,a2,type="n",main="W")
text(a1,a2,labels=paste("w",1:25,sep=""),col=rainbow(25), cex=0.8)
Latihan 4
• Create some programs to make the graph below, using 100 observation of Chi-square (4)
x <- rchisq(100,df=4)
hist(x,freq=FALSE,ylim=c(0,0.2))
curve(dchisq(x,df=4),col=2,lty=2,lwd=2,add=TRUE)
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)
Latihan 5
• Create 4 graph in one window from 100 random numbers which follow N(5,1), with format described bellow
windows()
yb <- rnorm(100,5,1)
split.screen(c(2,2))## [1] 1 2 3 4screen(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)
Pembangkitan Bilangan Acak
_
x <- rnorm(10) #x~N(0,1)
x## [1] -0.52296887 1.38055352 -1.10752438 0.69853952 -1.28688901 0.06337517
## [7] 1.93947651 0.06196071 -0.11692580 1.01794751x1 <- rnorm(10,3,2) #x1~N(3,sd=2)
x1## [1] 4.7211727 2.1141527 1.2948607 7.8784763 -0.1513073 0.9038301
## [7] 1.8865003 2.3865610 6.2497321 4.1430213x2 <- rbinom(10,1,0.4) #x2~bernoulli(0.4)
x2## [1] 0 1 1 1 0 0 0 1 0 1mencari nilai peluang kumulatif peubah acak
p1 <- pnorm(1.645) #P(Z<1.645)=0.95
p1## [1] 0.9500151p2 <- pnorm(1.96) #P(Z<1.975)=0.975
p2## [1] 0.9750021p3 <- pnorm(-1.96)
p3## [1] 0.0249979p4 <- pf(15,df1=10,df2=15)
p4## [1] 0.9999955mencari nilai kuantil dari peluang peubah acak
q1 <- qnorm(0.975)
q1## [1] 1.959964q2 <- qnorm(0.95,2,1) #X~N(2,1), P(X<x)=0.95
q2## [1] 3.644854mencari nilai density peubah acak plot density normal
a <- seq(-4,4,length=1000)
da <- dnorm(a)
plot(a,da)
Invers Transform Method
_
Latihan 1
_
• Membangkitkan bilangan acak Eksponensial
#Eksponensial(lambda=3)
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")
Acceptance-Rejection Method
_
_
Latihan 2
_
ar <- function(n,x0,x1,f) {
xx <- seq(x0,x1,length=10000)
c <- max(f(xx))
terima <- 0; iterasi <- 0
hasil <- numeric(n)
while(terima<n) {
x <- runif(1,x0,x1)
y1<- runif(1,0, c)
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,1,f)
yyy## $hasil
## [1] 0.8647212 0.7751099 0.8382877 0.7707715 0.5355970 0.8613824 0.2036477
## [8] 0.7979930 0.7438394 0.3443435 0.9837322 0.6935082 0.6331153 0.8880315
## [15] 0.7690405 0.6483695 0.8795432 0.9360689 0.7233519 0.7620444 0.9868082
## [22] 0.8760261 0.7240640 0.8140516 0.5588949 0.8900940 0.7456896 0.8480646
## [29] 0.8703302 0.8223331 0.8508123 0.7709219 0.8953595 0.5803863 0.5982260
## [36] 0.9235285 0.7367755 0.6898170 0.8301572 0.9293209 0.9095163 0.5347576
## [43] 0.3478601 0.8759762 0.7286815 0.8749293 0.6988356 0.8312562 0.5572238
## [50] 0.6647687 0.7400502 0.9806898 0.3800746 0.7553169 0.5184889 0.8879149
## [57] 0.9177773 0.8084086 0.8537441 0.4232184 0.7604306 0.3405763 0.3886568
## [64] 0.4774175 0.5387605 0.9485434 0.7124685 0.9081691 0.9457656 0.7716899
## [71] 0.6946655 0.5368832 0.8481593 0.8242752 0.5123742 0.3152032 0.9924487
## [78] 0.9327120 0.9892809 0.6283590 0.5254605 0.8810815 0.5291748 0.5765517
## [85] 0.7231807 0.8761180 0.3995670 0.8986123 0.9335217 0.7859216 0.7784128
## [92] 0.6955333 0.9060413 0.9916424 0.4729846 0.9770567 0.9386110 0.9959093
## [99] 0.8543663 0.8309547
##
## $iterasi
## [1] 322Direct Transformation
_
_
Analisis Regresi
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)
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.3124model2 <- 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-06model3 <- 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-10model4 <- 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-07model5 <- 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-09R2 <- 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.3298839coef(model3)## (Intercept) X1 X2
## -0.6288095 2.7295126 0.7245911confint(model3)## 2.5 % 97.5 %
## (Intercept) -28.4234482 27.1658292
## X1 1.8854322 3.5735929
## X2 0.5979779 0.8512044cbind(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.8512044anova(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 ' ' 1par(mfrow=c(2,2)) ; plot(model3)