# banyaknya amatan
n <- 1000
# membangkitan U(0,1)
u <- runif(n)
# fungsi inverse CDF
set.seed(10) # membuat seed
x <- -(log(1-u)/3) #log merupakan ln di R
head(x)
## [1] 0.70798113 0.08811037 0.28793867 0.01425366 0.44501461 0.12618293
# memeriksa banyaknya amatan
length(x)
## [1] 1000
set.seed(10) # membuat seed
y <- rexp(n,rate = 3)
head(y)
## [1] 0.004985469 0.306740402 0.250719646 0.525013950 0.077219539 0.362224334
# memeriksa banyaknya amatan
length(y)
## [1] 1000
par(mfrow=c(1,2))
hist(x,main="Exp dari Inverse Transform")
hist(y,main="Exp dari fungsi rexp")
# banyaknya amatan
n=1000
set.seed(10)
y = runif(n,0,1)
set.seed(30)
u = runif(n,0,1)
h <- function(y) 6*y*(1-y)
derivH <- optimize(f=h,interval=c(0,1),maximum=T)
derivH$objective
## [1] 1.5
f <- function(y) 6*y*(1-y)
g <- function(y) 1
nilaic <- derivH$objective
#kriteria penerimaan
criteria <- u < f(y)/(nilaic*g(y))
head(criteria)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE
#jadikan Y sebagai X
x <- y[criteria]
head(x)
## [1] 0.50747820 0.30676851 0.42690767 0.69310208 0.08513597 0.22543662
# banyaknya amatan
length(x)
## [1] 668
z <- rbeta(n,shape1 = 2,shape2 = 2)
par(mfrow=c(1,2))
hist(x,main="beta(2,2) dari AR")
hist(z,main="beta(2,2) dari rbeta")
# banyaknya amatan
n=1500
set.seed(10)
y = rexp(n,rate = 2/3)
set.seed(30)
u = runif(n,0,1)
h <- function(y) (-1) * (3/(2*gamma(3/2))) * (y)^(1/2) * exp(-y/3)
derivH <- optim(par = 0,fn=h,lower = 0,method = "L-BFGS-B")
derivH$value
## [1] -1.257317
f <- function(y) (3/(2*gamma(3/2))) * (y)^(1/2) * exp(-y)
g <- function(y) (2/3) * exp(-(2/3) * y)
nilaic <- -derivH$value #jangan lupa tambahkan negatif
#kriteria penerimaan
criteria <- u < f(y)/(nilaic*g(y))
head(criteria)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE
#jadikan Y sebagai X
x <- y[criteria]
head(x)
## [1] 0.02243461 1.38033181 1.12823841 2.36256278 0.34748792 1.63000950
# banyaknya amatan
length(x)
## [1] 1393
#karena >1000 ambil 1000 saja
x <- x[1:1000]
z <- rgamma(n,shape = 3/2,scale = 1)
par(mfrow=c(1,2))
hist(x,main="gamma(3/2,1) dari AR")
hist(z,main="beta(3/2,1) dari rgamma")
#banyaknya amatan
n <- 1000
set.seed(15)
z <- rnorm(n,mean = 0,sd = 1)
y <- rchisq(n,df = 2)
x <- z/sqrt(y/2)
head(x)
## [1] 0.4249926 2.2981978 -0.3956900 0.8938340 0.9930945 -1.4325087
set.seed(5)
v <- rt(n,df=2)
par(mfrow=c(1,2))
hist(x,main="t(2) dari Direct")
hist(v,main="t(2) dari rt")
#banyaknya amatan
n <- 1000
set.seed(12)
epsilon <- rnorm(n,mean = 0,sd = 1)
set.seed(10)
X <- replicate(3,rnorm(n,mean=0,sd=1.5))
head(X)
## [,1] [,2] [,3]
## [1,] 0.02811926 1.5750205 -0.46179747
## [2,] -0.27637881 0.4291389 1.13712837
## [3,] -2.05699582 0.3608472 -0.86079512
## [4,] -0.89875157 1.2490578 -1.40811669
## [5,] 0.44181769 -0.3344748 -0.04154899
## [6,] 0.58469145 0.4325163 -1.59937298
betas <- matrix(c(10,3,5,7),nrow = 4,ncol = 1)
Xgab <- cbind(1,X)
Y <- Xgab%*%betas+epsilon
dataRegresi <- data.frame(Y,X)
colnames(dataRegresi) <- c("Y","X1","X2","X3")
head(dataRegresi)
## Y X1 X2 X3
## 1 13.246310 0.02811926 1.5750205 -0.46179747
## 2 20.853626 -0.27637881 0.4291389 1.13712837
## 3 -1.349062 -2.05699582 0.3608472 -0.86079512
## 4 2.772212 -0.89875157 1.2490578 -1.40811669
## 5 7.364594 0.44181769 -0.3344748 -0.04154899
## 6 2.448749 0.58469145 0.4325163 -1.59937298
modelRegresi <- lm(Y~.,data = dataRegresi)
summary(modelRegresi)
##
## Call:
## lm(formula = Y ~ ., data = dataRegresi)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0364 -0.6052 -0.0207 0.6087 3.1774
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.97484 0.03033 328.8 <2e-16 ***
## X1 2.97225 0.02042 145.6 <2e-16 ***
## X2 4.97282 0.01939 256.5 <2e-16 ***
## X3 7.00493 0.01989 352.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9589 on 996 degrees of freedom
## Multiple R-squared: 0.9955, Adjusted R-squared: 0.9954
## F-statistic: 7.275e+04 on 3 and 996 DF, p-value: < 2.2e-16
#banyaknya amatan
n <- 1000
ulangan <- 100
# membuat tempat penyimpanan
# koefisien tiap ulangan
koefisien_ulangan <- NULL
set.seed(13)
for(i in seq(ulangan)){
epsilon <- rnorm(n,mean = 0,sd = 1)
X <- replicate(3,rnorm(n,mean=0,sd=1.5))
Xgab <- cbind(1,X)
Y <- Xgab%*%betas+epsilon
dataRegresi <- data.frame(Y,X)
colnames(dataRegresi) <- c("Y","X1","X2","X3")
modelRegresi <- lm(Y~.,data = dataRegresi)
# rbind digunakan untuk menggabungkan baris
koefisien_ulangan <- rbind(koefisien_ulangan,coef(modelRegresi))
}
head(koefisien_ulangan)
## (Intercept) X1 X2 X3
## [1,] 9.998451 2.997131 4.961536 7.019022
## [2,] 9.999899 2.981410 5.012151 6.977783
## [3,] 9.971732 2.986735 4.992031 7.019945
## [4,] 9.940730 3.013790 5.006416 7.003125
## [5,] 10.029303 3.033154 5.009336 6.972909
## [6,] 10.008416 3.019092 5.010800 6.999371
#menghitung mean dari setiap koefisien
colMeans(koefisien_ulangan)
## (Intercept) X1 X2 X3
## 9.997236 3.002950 5.000115 7.001103
#menghitung sd dari setiap koefisien
apply(koefisien_ulangan,2,sd)
## (Intercept) X1 X2 X3
## 0.02972359 0.01960207 0.01923794 0.02379381