########### Computational Statistics Final Exam ###########
############### 1. Newton Raphson Algorithm ###############
setwd('D:\\R Code')
# Input data
x <- c(3.17693741992116,
10,
5.18462506122887,
2.56333450908598,
3.40210871770978,
10,
2.18819495498594,
0.317795732254817,
10,
9.15132800889364)
summary(x)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.3178 2.7170 4.2930 5.5980 9.7880 10.0000
# Functions
# Original loglike function
loglike <- function(theta, dat){
y <- dat
m = 7
l <- rep(0, length(theta) )
for (i in 1:length(theta)){
l[i] = m*log(theta[i]) - theta[i] * sum(y)
}
return(l)
}
# First derivative
loglike1 <- function(theta, dat){
y <- dat
m = 7
l <- rep(0, length(theta) )
for (i in 1:length(theta)){
l[i] = m/theta[i] - sum(y)
}
return(l)
}
# Second derivative
loglike2 <- function(theta, dat){
y <- dat
m = 7
l <- rep(0, length(theta) )
for (i in 1:length(theta)){
l[i] = -m/((theta[i])^2)
}
return(l)
}
Newton = function(theta, dat, iter=50){
for(i in 1:iter){
theta = theta - loglike1(theta, dat)/loglike2(theta,dat)
print(c(i,theta))
}
return(theta)
}
# Plot loglike, loglike1, loglike2
theta <- seq(0.005,1,by=0.01) # theta > 0
plot(theta, loglike(theta, dat=x), type='l', xlab="theta", main="Loglike Function")
plot(theta, loglike1(theta, dat=x), type='l', xlab="theta", main="First Derivative of Loglike Function" )
plot(theta, loglike2(theta, dat=x), type='l', xlab="theta", main="Second Derivative of Loglike Function" )
min(loglike1(theta,dat=x))
## [1] -48.94915
theta <- seq(0.005,50,by=0.01) # theta > 0
plot(theta, loglike(theta, dat=x), type='l', xlab="theta", main="Loglike Function")
plot(theta, loglike1(theta, dat=x), type='l', xlab="theta", main="First Derivative of Loglike Function" )
plot(theta, loglike2(theta, dat=x), type='l', xlab="theta", main="Second Derivative of Loglike Function" )
min(loglike1(theta,dat=x))
## [1] -55.84431
# Find MLE using the Newton-Raphson Method with different starting points
Newton(mean(x), dat=x)
## [1] 1.0000 -239.4725
## [1] 2.0 -459127.2
## [1] 3.000000e+00 -1.685911e+12
## [1] 4.000000e+00 -2.273201e+25
## [1] 5.000000e+00 -4.132796e+51
## [1] 6.000000e+00 -1.366018e+104
## [1] 7.000000e+00 -1.492386e+209
## [1] 8 -Inf
## [1] 9 -Inf
## [1] 10 -Inf
## [1] 11 -Inf
## [1] 12 -Inf
## [1] 13 -Inf
## [1] 14 -Inf
## [1] 15 -Inf
## [1] 16 -Inf
## [1] 17 -Inf
## [1] 18 -Inf
## [1] 19 -Inf
## [1] 20 -Inf
## [1] 21 -Inf
## [1] 22 -Inf
## [1] 23 -Inf
## [1] 24 -Inf
## [1] 25 -Inf
## [1] 26 -Inf
## [1] 27 -Inf
## [1] 28 -Inf
## [1] 29 -Inf
## [1] 30 -Inf
## [1] 31 -Inf
## [1] 32 -Inf
## [1] 33 -Inf
## [1] 34 -Inf
## [1] 35 -Inf
## [1] 36 -Inf
## [1] 37 -Inf
## [1] 38 -Inf
## [1] 39 -Inf
## [1] 40 -Inf
## [1] 41 -Inf
## [1] 42 -Inf
## [1] 43 -Inf
## [1] 44 -Inf
## [1] 45 -Inf
## [1] 46 -Inf
## [1] 47 -Inf
## [1] 48 -Inf
## [1] 49 -Inf
## [1] 50 -Inf
## [1] -Inf
Newton(0.01, dat=x)
## [1] 1.00000000 0.01920022
## [1] 2.00000000 0.03545208
## [1] 3.00000000 0.06085218
## [1] 4.00000000 0.09208875
## [1] 5.0000000 0.1163538
## [1] 6.0000000 0.1244323
## [1] 7.0000000 0.1250321
## [1] 8.000000 0.125035
## [1] 9.000000 0.125035
## [1] 10.000000 0.125035
## [1] 11.000000 0.125035
## [1] 12.000000 0.125035
## [1] 13.000000 0.125035
## [1] 14.000000 0.125035
## [1] 15.000000 0.125035
## [1] 16.000000 0.125035
## [1] 17.000000 0.125035
## [1] 18.000000 0.125035
## [1] 19.000000 0.125035
## [1] 20.000000 0.125035
## [1] 21.000000 0.125035
## [1] 22.000000 0.125035
## [1] 23.000000 0.125035
## [1] 24.000000 0.125035
## [1] 25.000000 0.125035
## [1] 26.000000 0.125035
## [1] 27.000000 0.125035
## [1] 28.000000 0.125035
## [1] 29.000000 0.125035
## [1] 30.000000 0.125035
## [1] 31.000000 0.125035
## [1] 32.000000 0.125035
## [1] 33.000000 0.125035
## [1] 34.000000 0.125035
## [1] 35.000000 0.125035
## [1] 36.000000 0.125035
## [1] 37.000000 0.125035
## [1] 38.000000 0.125035
## [1] 39.000000 0.125035
## [1] 40.000000 0.125035
## [1] 41.000000 0.125035
## [1] 42.000000 0.125035
## [1] 43.000000 0.125035
## [1] 44.000000 0.125035
## [1] 45.000000 0.125035
## [1] 46.000000 0.125035
## [1] 47.000000 0.125035
## [1] 48.000000 0.125035
## [1] 49.000000 0.125035
## [1] 50.000000 0.125035
## [1] 0.125035
Newton(0.1, dat=x)
## [1] 1.0000000 0.1200224
## [1] 2.000000 0.124834
## [1] 3.0000000 0.1250347
## [1] 4.000000 0.125035
## [1] 5.000000 0.125035
## [1] 6.000000 0.125035
## [1] 7.000000 0.125035
## [1] 8.000000 0.125035
## [1] 9.000000 0.125035
## [1] 10.000000 0.125035
## [1] 11.000000 0.125035
## [1] 12.000000 0.125035
## [1] 13.000000 0.125035
## [1] 14.000000 0.125035
## [1] 15.000000 0.125035
## [1] 16.000000 0.125035
## [1] 17.000000 0.125035
## [1] 18.000000 0.125035
## [1] 19.000000 0.125035
## [1] 20.000000 0.125035
## [1] 21.000000 0.125035
## [1] 22.000000 0.125035
## [1] 23.000000 0.125035
## [1] 24.000000 0.125035
## [1] 25.000000 0.125035
## [1] 26.000000 0.125035
## [1] 27.000000 0.125035
## [1] 28.000000 0.125035
## [1] 29.000000 0.125035
## [1] 30.000000 0.125035
## [1] 31.000000 0.125035
## [1] 32.000000 0.125035
## [1] 33.000000 0.125035
## [1] 34.000000 0.125035
## [1] 35.000000 0.125035
## [1] 36.000000 0.125035
## [1] 37.000000 0.125035
## [1] 38.000000 0.125035
## [1] 39.000000 0.125035
## [1] 40.000000 0.125035
## [1] 41.000000 0.125035
## [1] 42.000000 0.125035
## [1] 43.000000 0.125035
## [1] 44.000000 0.125035
## [1] 45.000000 0.125035
## [1] 46.000000 0.125035
## [1] 47.000000 0.125035
## [1] 48.000000 0.125035
## [1] 49.000000 0.125035
## [1] 50.000000 0.125035
## [1] 0.125035
Newton(0.2, dat=x)
## [1] 1.00000000 0.08008957
## [1] 2.0000000 0.1088788
## [1] 3.0000000 0.1229474
## [1] 4.0000000 0.1250001
## [1] 5.000000 0.125035
## [1] 6.000000 0.125035
## [1] 7.000000 0.125035
## [1] 8.000000 0.125035
## [1] 9.000000 0.125035
## [1] 10.000000 0.125035
## [1] 11.000000 0.125035
## [1] 12.000000 0.125035
## [1] 13.000000 0.125035
## [1] 14.000000 0.125035
## [1] 15.000000 0.125035
## [1] 16.000000 0.125035
## [1] 17.000000 0.125035
## [1] 18.000000 0.125035
## [1] 19.000000 0.125035
## [1] 20.000000 0.125035
## [1] 21.000000 0.125035
## [1] 22.000000 0.125035
## [1] 23.000000 0.125035
## [1] 24.000000 0.125035
## [1] 25.000000 0.125035
## [1] 26.000000 0.125035
## [1] 27.000000 0.125035
## [1] 28.000000 0.125035
## [1] 29.000000 0.125035
## [1] 30.000000 0.125035
## [1] 31.000000 0.125035
## [1] 32.000000 0.125035
## [1] 33.000000 0.125035
## [1] 34.000000 0.125035
## [1] 35.000000 0.125035
## [1] 36.000000 0.125035
## [1] 37.000000 0.125035
## [1] 38.000000 0.125035
## [1] 39.000000 0.125035
## [1] 40.000000 0.125035
## [1] 41.000000 0.125035
## [1] 42.000000 0.125035
## [1] 43.000000 0.125035
## [1] 44.000000 0.125035
## [1] 45.000000 0.125035
## [1] 46.000000 0.125035
## [1] 47.000000 0.125035
## [1] 48.000000 0.125035
## [1] 49.000000 0.125035
## [1] 50.000000 0.125035
## [1] 0.125035
Newton(0.25, dat=x)
## [1] 1.0000000000 0.0001399607
## [1] 2.0000000000 0.0002797647
## [1] 3.0000000000 0.0005589034
## [1] 4.000000000 0.001115309
## [1] 5.000000000 0.002220669
## [1] 6.000000000 0.004401897
## [1] 7.000000000 0.008648824
## [1] 8.0000000 0.0166994
## [1] 9.00000000 0.03116846
## [1] 10.00000000 0.05456732
## [1] 11.00000000 0.08532056
## [1] 12.0000000 0.1124206
## [1] 13.0000000 0.1237624
## [1] 14.000000 0.125022
## [1] 15.000000 0.125035
## [1] 16.000000 0.125035
## [1] 17.000000 0.125035
## [1] 18.000000 0.125035
## [1] 19.000000 0.125035
## [1] 20.000000 0.125035
## [1] 21.000000 0.125035
## [1] 22.000000 0.125035
## [1] 23.000000 0.125035
## [1] 24.000000 0.125035
## [1] 25.000000 0.125035
## [1] 26.000000 0.125035
## [1] 27.000000 0.125035
## [1] 28.000000 0.125035
## [1] 29.000000 0.125035
## [1] 30.000000 0.125035
## [1] 31.000000 0.125035
## [1] 32.000000 0.125035
## [1] 33.000000 0.125035
## [1] 34.000000 0.125035
## [1] 35.000000 0.125035
## [1] 36.000000 0.125035
## [1] 37.000000 0.125035
## [1] 38.000000 0.125035
## [1] 39.000000 0.125035
## [1] 40.000000 0.125035
## [1] 41.000000 0.125035
## [1] 42.000000 0.125035
## [1] 43.000000 0.125035
## [1] 44.000000 0.125035
## [1] 45.000000 0.125035
## [1] 46.000000 0.125035
## [1] 47.000000 0.125035
## [1] 48.000000 0.125035
## [1] 49.000000 0.125035
## [1] 50.000000 0.125035
## [1] 0.125035
Newton(0.26, dat=x)
## [1] 1.00000000 -0.02064862
## [1] 2.00000000 -0.04470721
## [1] 3.0000000 -0.1053998
## [1] 4.0000000 -0.2996477
## [1] 5.000000 -1.317404
## [1] 6.00000 -16.51535
## [1] 7.000 -2214.475
## [1] 8 -39224658
## [1] 9.000000e+00 -1.230514e+16
## [1] 1.000000e+01 -1.210994e+33
## [1] 1.100000e+01 -1.172876e+67
## [1] 1.200000e+01 -1.100203e+135
## [1] 1.300000e+01 -9.680856e+270
## [1] 14 -Inf
## [1] 15 -Inf
## [1] 16 -Inf
## [1] 17 -Inf
## [1] 18 -Inf
## [1] 19 -Inf
## [1] 20 -Inf
## [1] 21 -Inf
## [1] 22 -Inf
## [1] 23 -Inf
## [1] 24 -Inf
## [1] 25 -Inf
## [1] 26 -Inf
## [1] 27 -Inf
## [1] 28 -Inf
## [1] 29 -Inf
## [1] 30 -Inf
## [1] 31 -Inf
## [1] 32 -Inf
## [1] 33 -Inf
## [1] 34 -Inf
## [1] 35 -Inf
## [1] 36 -Inf
## [1] 37 -Inf
## [1] 38 -Inf
## [1] 39 -Inf
## [1] 40 -Inf
## [1] 41 -Inf
## [1] 42 -Inf
## [1] 43 -Inf
## [1] 44 -Inf
## [1] 45 -Inf
## [1] 46 -Inf
## [1] 47 -Inf
## [1] 48 -Inf
## [1] 49 -Inf
## [1] 50 -Inf
## [1] -Inf
Newton(0.5, dat=x)
## [1] 1.0000000 -0.9994402
## [1] 2.000000 -9.987688
## [1] 3.0000 -817.7834
## [1] 4 -5350295
## [1] 5.000000e+00 -2.289412e+14
## [1] 6.00000e+00 -4.19195e+29
## [1] 7.000000e+00 -1.405402e+60
## [1] 8.000000e+00 -1.579682e+121
## [1] 9.000000e+00 -1.995758e+243
## [1] 10 -Inf
## [1] 11 -Inf
## [1] 12 -Inf
## [1] 13 -Inf
## [1] 14 -Inf
## [1] 15 -Inf
## [1] 16 -Inf
## [1] 17 -Inf
## [1] 18 -Inf
## [1] 19 -Inf
## [1] 20 -Inf
## [1] 21 -Inf
## [1] 22 -Inf
## [1] 23 -Inf
## [1] 24 -Inf
## [1] 25 -Inf
## [1] 26 -Inf
## [1] 27 -Inf
## [1] 28 -Inf
## [1] 29 -Inf
## [1] 30 -Inf
## [1] 31 -Inf
## [1] 32 -Inf
## [1] 33 -Inf
## [1] 34 -Inf
## [1] 35 -Inf
## [1] 36 -Inf
## [1] 37 -Inf
## [1] 38 -Inf
## [1] 39 -Inf
## [1] 40 -Inf
## [1] 41 -Inf
## [1] 42 -Inf
## [1] 43 -Inf
## [1] 44 -Inf
## [1] 45 -Inf
## [1] 46 -Inf
## [1] 47 -Inf
## [1] 48 -Inf
## [1] 49 -Inf
## [1] 50 -Inf
## [1] -Inf
Newton(2, dat=x)
## [1] 1.00000 -27.99104
## [1] 2.000 -6322.215
## [1] 3 -319686379
## [1] 4.000000e+00 -8.173662e+17
## [1] 5.000000e+00 -5.343204e+36
## [1] 6.000000e+00 -2.283347e+74
## [1] 7.000000e+00 -4.169771e+149
## [1] 8.00000e+00 -1.39057e+300
## [1] 9 -Inf
## [1] 10 -Inf
## [1] 11 -Inf
## [1] 12 -Inf
## [1] 13 -Inf
## [1] 14 -Inf
## [1] 15 -Inf
## [1] 16 -Inf
## [1] 17 -Inf
## [1] 18 -Inf
## [1] 19 -Inf
## [1] 20 -Inf
## [1] 21 -Inf
## [1] 22 -Inf
## [1] 23 -Inf
## [1] 24 -Inf
## [1] 25 -Inf
## [1] 26 -Inf
## [1] 27 -Inf
## [1] 28 -Inf
## [1] 29 -Inf
## [1] 30 -Inf
## [1] 31 -Inf
## [1] 32 -Inf
## [1] 33 -Inf
## [1] 34 -Inf
## [1] 35 -Inf
## [1] 36 -Inf
## [1] 37 -Inf
## [1] 38 -Inf
## [1] 39 -Inf
## [1] 40 -Inf
## [1] 41 -Inf
## [1] 42 -Inf
## [1] 43 -Inf
## [1] 44 -Inf
## [1] 45 -Inf
## [1] 46 -Inf
## [1] 47 -Inf
## [1] 48 -Inf
## [1] 49 -Inf
## [1] 50 -Inf
## [1] -Inf
Newton(10, dat=x)
## [1] 1.0000 -779.7761
## [1] 2 -4864604
## [1] 3.00000e+00 -1.89262e+14
## [1] 4.000000e+00 -2.864805e+29
## [1] 5.000000e+00 -6.563848e+59
## [1] 6.000000e+00 -3.445764e+120
## [1] 7.00000e+00 -9.49597e+241
## [1] 8 -Inf
## [1] 9 -Inf
## [1] 10 -Inf
## [1] 11 -Inf
## [1] 12 -Inf
## [1] 13 -Inf
## [1] 14 -Inf
## [1] 15 -Inf
## [1] 16 -Inf
## [1] 17 -Inf
## [1] 18 -Inf
## [1] 19 -Inf
## [1] 20 -Inf
## [1] 21 -Inf
## [1] 22 -Inf
## [1] 23 -Inf
## [1] 24 -Inf
## [1] 25 -Inf
## [1] 26 -Inf
## [1] 27 -Inf
## [1] 28 -Inf
## [1] 29 -Inf
## [1] 30 -Inf
## [1] 31 -Inf
## [1] 32 -Inf
## [1] 33 -Inf
## [1] 34 -Inf
## [1] 35 -Inf
## [1] 36 -Inf
## [1] 37 -Inf
## [1] 38 -Inf
## [1] 39 -Inf
## [1] 40 -Inf
## [1] 41 -Inf
## [1] 42 -Inf
## [1] 43 -Inf
## [1] 44 -Inf
## [1] 45 -Inf
## [1] 46 -Inf
## [1] 47 -Inf
## [1] 48 -Inf
## [1] 49 -Inf
## [1] 50 -Inf
## [1] -Inf
Newton(50, dat=x)
## [1] 1.0 -19894.4
## [1] 2 -3165451186
## [1] 3.000000e+00 -8.013821e+19
## [1] 4.000000e+00 -5.136268e+40
## [1] 5.000000e+00 -2.109909e+82
## [1] 6.000000e+00 -3.560377e+165
## [1] 7 -Inf
## [1] 8 -Inf
## [1] 9 -Inf
## [1] 10 -Inf
## [1] 11 -Inf
## [1] 12 -Inf
## [1] 13 -Inf
## [1] 14 -Inf
## [1] 15 -Inf
## [1] 16 -Inf
## [1] 17 -Inf
## [1] 18 -Inf
## [1] 19 -Inf
## [1] 20 -Inf
## [1] 21 -Inf
## [1] 22 -Inf
## [1] 23 -Inf
## [1] 24 -Inf
## [1] 25 -Inf
## [1] 26 -Inf
## [1] 27 -Inf
## [1] 28 -Inf
## [1] 29 -Inf
## [1] 30 -Inf
## [1] 31 -Inf
## [1] 32 -Inf
## [1] 33 -Inf
## [1] 34 -Inf
## [1] 35 -Inf
## [1] 36 -Inf
## [1] 37 -Inf
## [1] 38 -Inf
## [1] 39 -Inf
## [1] 40 -Inf
## [1] 41 -Inf
## [1] 42 -Inf
## [1] 43 -Inf
## [1] 44 -Inf
## [1] 45 -Inf
## [1] 46 -Inf
## [1] 47 -Inf
## [1] 48 -Inf
## [1] 49 -Inf
## [1] 50 -Inf
## [1] -Inf
############### 2. Logistic Regression with IRLS ###############
# Iteratively Reweighted Least Squares
IRLS = function(beta.init, iter, X, y, X.nrow){
beta.values = matrix(0, iter+1, 4)
cov.values = matrix(0, iter, 16)
beta.values[1,] = beta.init
W = matrix(0, X.nrow, X.nrow)
for(i in 1:iter){
pi = exp(X %*% beta.init)/(1 + exp(X %*% beta.init))
for(j in 1:X.nrow) W[j,j] = pi[j]
cov.matrix = solve(t(X) %*% W %*% X)
beta.init = beta.init + cov.matrix %*% t(X) %*% (y - pi)
beta.values[i+1,] = beta.init
cov.values[i,] = c(cov.matrix[1,],cov.matrix[2,],cov.matrix[3,],cov.matrix[4,])
}
list(beta=beta.values[2:(iter+1),], cov=cov.values)
}
# Read Data
data = read.table("logistic.txt");
X = matrix(1,nrow(data),4); y = rep(NA,nrow(data));
for(i in 1:nrow(data)){
X[i,2] = data[i,2]
X[i,3] = data[i,3]
X[i,4] = data[i,4]
y[i] = data[i,1]
}
# Running Function
beta.init = c(-4.0,0.0,0.5,-1.0); iter = 100; X.nrow = nrow(data);
result = IRLS(beta.init, iter, X, y, X.nrow);
result$beta
## [,1] [,2] [,3] [,4]
## [1,] -35.111223 0.036000268 0.7628247 17.94253279
## [2,] -36.236029 0.036393376 0.9111047 17.82662722
## [3,] -37.344215 0.036774537 1.0583220 17.70873760
## [4,] -38.432693 0.037142161 1.2042287 17.58843247
## [5,] -39.498939 0.037494897 1.3486865 17.46529783
## [6,] -40.541567 0.037831803 1.4917666 17.33898393
## [7,] -41.560286 0.038152351 1.6337399 17.20919824
## [8,] -42.555009 0.038456185 1.7749317 17.07562372
## [9,] -43.524971 0.038742801 1.9155879 16.93783744
## [10,] -44.468703 0.039011468 2.0558718 16.79532157
## [11,] -45.384815 0.039261478 2.1959528 16.64757318
## [12,] -46.272844 0.039492532 2.3360754 16.49423034
## [13,] -47.133485 0.039704980 2.4765398 16.33512460
## [14,] -47.968170 0.039899808 2.6176150 16.17024164
## [15,] -48.778414 0.040078428 2.7594557 15.99963913
## [16,] -49.565341 0.040242390 2.9020753 15.82337888
## [17,] -50.329555 0.040393070 3.0453851 15.64149967
## [18,] -51.071274 0.040531369 3.1892786 15.45402695
## [19,] -51.790587 0.040657549 3.3337116 15.26099944
## [20,] -52.487708 0.040771355 3.4787329 15.06249439
## [21,] -53.163172 0.040872444 3.6244541 14.85864069
## [22,] -53.817746 0.040960709 3.7709790 14.64960638
## [23,] -54.451937 0.041036182 3.9183145 14.43555549
## [24,] -55.065392 0.041098740 4.0662955 14.21659841
## [25,] -55.656468 0.041147852 4.2145367 13.99275333
## [26,] -56.221813 0.041182252 4.3624070 13.76390750
## [27,] -56.756007 0.041199526 4.5090306 13.52978084
## [28,] -57.251629 0.041195992 4.6533325 13.28992342
## [29,] -57.699954 0.041167107 4.7941203 13.04376754
## [30,] -58.092035 0.041108247 4.9301823 12.79071718
## [31,] -58.419673 0.041015386 5.0603716 12.53023178
## [32,] -58.675818 0.040885341 5.1836400 12.26186886
## [33,] -58.854410 0.040715605 5.2990228 11.98528424
## [34,] -58.950173 0.040504017 5.4056479 11.70021720
## [35,] -58.958944 0.040248686 5.5028433 11.40649152
## [36,] -58.878510 0.039948345 5.5902937 11.10403733
## [37,] -58.708930 0.039602591 5.6680857 10.79289711
## [38,] -58.451520 0.039211403 5.7365415 10.47318155
## [39,] -58.107098 0.038774318 5.7959281 10.14499518
## [40,] -57.674593 0.038289885 5.8462104 9.80837778
## [41,] -57.150703 0.037755666 5.8869647 9.46328939
## [42,] -56.530567 0.037168727 5.9174597 9.10964001
## [43,] -55.808809 0.036526236 5.9368024 8.74733887
## [44,] -54.980083 0.035825702 5.9440245 8.37632938
## [45,] -54.039120 0.035064823 5.9381099 7.99660219
## [46,] -52.980842 0.034241396 5.9180336 7.60820644
## [47,] -51.800796 0.033353509 5.8828315 7.21127399
## [48,] -50.495478 0.032399781 5.8316252 6.80604657
## [49,] -49.061938 0.031379244 5.7635248 6.39288721
## [50,] -47.496706 0.030290863 5.6774183 5.97227714
## [51,] -45.794680 0.029133054 5.5717448 5.54482064
## [52,] -43.948762 0.027903580 5.4443852 5.11127777
## [53,] -41.950994 0.026600128 5.2928260 4.67262563
## [54,] -39.795478 0.025221614 5.1146831 4.23013404
## [55,] -37.481848 0.023769565 4.9083670 3.78543511
## [56,] -35.017299 0.022248660 4.6735338 3.34056179
## [57,] -32.416568 0.020666228 4.4111876 2.89794856
## [58,] -29.700060 0.019031046 4.1234216 2.46042412
## [59,] -26.890944 0.017351892 3.8128981 2.03124242
## [60,] -24.013487 0.015636980 3.4824775 1.61418184
## [61,] -21.095395 0.013895636 3.1355151 1.21370411
## [62,] -18.174571 0.012142040 2.7770128 0.83513900
## [63,] -15.306624 0.010398382 2.4151377 0.48488662
## [64,] -12.567875 0.008695828 2.0619905 0.17064216
## [65,] -10.053601 0.007077938 1.7327899 -0.09849666
## [66,] -7.874799 0.005611536 1.4435776 -0.31238050
## [67,] -6.144123 0.004389173 1.2080704 -0.46185422
## [68,] -4.929645 0.003490496 1.0342509 -0.54542654
## [69,] -4.194261 0.002920693 0.9200908 -0.57707337
## [70,] -3.804949 0.002604657 0.8529566 -0.58083162
## [71,] -3.616331 0.002444572 0.8164416 -0.57560382
## [72,] -3.528387 0.002366960 0.7973648 -0.56985890
## [73,] -3.487508 0.002329729 0.7875368 -0.56575388
## [74,] -3.468208 0.002311738 0.7824773 -0.56322543
## [75,] -3.458893 0.002302917 0.7798618 -0.56177241
## [76,] -3.454297 0.002298525 0.7785030 -0.56096736
## [77,] -3.451988 0.002296307 0.7777942 -0.56053040
## [78,] -3.450811 0.002295175 0.7774233 -0.56029608
## [79,] -3.450206 0.002294591 0.7772288 -0.56017135
## [80,] -3.449892 0.002294289 0.7771267 -0.56010524
## [81,] -3.449728 0.002294132 0.7770730 -0.56007031
## [82,] -3.449643 0.002294050 0.7770448 -0.56005188
## [83,] -3.449598 0.002294007 0.7770300 -0.56004217
## [84,] -3.449574 0.002293984 0.7770222 -0.56003706
## [85,] -3.449562 0.002293973 0.7770181 -0.56003437
## [86,] -3.449556 0.002293966 0.7770160 -0.56003296
## [87,] -3.449552 0.002293963 0.7770148 -0.56003221
## [88,] -3.449550 0.002293961 0.7770142 -0.56003182
## [89,] -3.449549 0.002293961 0.7770139 -0.56003162
## [90,] -3.449549 0.002293960 0.7770138 -0.56003151
## [91,] -3.449549 0.002293960 0.7770137 -0.56003145
## [92,] -3.449549 0.002293960 0.7770136 -0.56003142
## [93,] -3.449548 0.002293960 0.7770136 -0.56003140
## [94,] -3.449548 0.002293960 0.7770136 -0.56003140
## [95,] -3.449548 0.002293960 0.7770136 -0.56003139
## [96,] -3.449548 0.002293960 0.7770136 -0.56003139
## [97,] -3.449548 0.002293960 0.7770136 -0.56003139
## [98,] -3.449548 0.002293960 0.7770136 -0.56003139
## [99,] -3.449548 0.002293960 0.7770136 -0.56003139
## [100,] -3.449548 0.002293960 0.7770136 -0.56003139
result$cov
## [,1] [,2] [,3] [,4] [,5]
## [1,] 20.1183171 -4.072544e-03 -4.61036173 -0.81335512 -4.072544e-03
## [2,] 0.2460934 -5.762827e-05 -0.05453021 -0.00983559 -5.762827e-05
## [3,] 0.2482481 -5.888648e-05 -0.05475095 -0.01005653 -5.888648e-05
## [4,] 0.2507770 -6.020783e-05 -0.05503809 -0.01031425 -6.020783e-05
## [5,] 0.2535742 -6.154285e-05 -0.05537671 -0.01060006 -6.154285e-05
## [6,] 0.2564216 -6.282698e-05 -0.05572879 -0.01089698 -6.282698e-05
## [7,] 0.2591002 -6.400280e-05 -0.05605436 -0.01118849 -6.400280e-05
## [8,] 0.2615337 -6.505011e-05 -0.05633835 -0.01147014 -6.505011e-05
## [9,] 0.2637951 -6.598910e-05 -0.05659171 -0.01174986 -6.598910e-05
## [10,] 0.2659928 -6.685133e-05 -0.05683248 -0.01203710 -6.685133e-05
## [11,] 0.2681700 -6.764948e-05 -0.05707032 -0.01233289 -6.764948e-05
## [12,] 0.2702907 -6.837016e-05 -0.05730464 -0.01262873 -6.837016e-05
## [13,] 0.2722953 -6.898744e-05 -0.05753257 -0.01291329 -6.898744e-05
## [14,] 0.2741603 -6.948031e-05 -0.05775699 -0.01317962 -6.948031e-05
## [15,] 0.2759180 -6.984342e-05 -0.05798834 -0.01342771 -6.984342e-05
## [16,] 0.2776392 -7.009385e-05 -0.05823994 -0.01366251 -7.009385e-05
## [17,] 0.2794006 -7.027926e-05 -0.05852022 -0.01389041 -7.027926e-05
## [18,] 0.2812570 -7.048266e-05 -0.05882626 -0.01411624 -7.048266e-05
## [19,] 0.2832262 -7.080562e-05 -0.05914373 -0.01434197 -7.080562e-05
## [20,] 0.2852809 -7.131822e-05 -0.05945341 -0.01456628 -7.131822e-05
## [21,] 0.2873490 -7.201103e-05 -0.05973964 -0.01478495 -7.201103e-05
## [22,] 0.2893607 -7.281792e-05 -0.05999708 -0.01499411 -7.281792e-05
## [23,] 0.2913253 -7.369357e-05 -0.06023490 -0.01519522 -7.369357e-05
## [24,] 0.2933513 -7.464630e-05 -0.06047591 -0.01539626 -7.464630e-05
## [25,] 0.2956163 -7.572018e-05 -0.06075259 -0.01560967 -7.572018e-05
## [26,] 0.2983501 -7.699125e-05 -0.06110356 -0.01585156 -7.699125e-05
## [27,] 0.3018166 -7.856894e-05 -0.06156863 -0.01614093 -7.856894e-05
## [28,] 0.3062528 -8.056004e-05 -0.06218139 -0.01649552 -8.056004e-05
## [29,] 0.3117949 -8.301339e-05 -0.06296273 -0.01692584 -8.301339e-05
## [30,] 0.3184406 -8.589470e-05 -0.06391761 -0.01743188 -8.589470e-05
## [31,] 0.3260740 -8.911353e-05 -0.06503617 -0.01800452 -8.911353e-05
## [32,] 0.3345497 -9.258316e-05 -0.06630248 -0.01863116 -9.258316e-05
## [33,] 0.3437711 -9.627024e-05 -0.06770357 -0.01930110 -9.627024e-05
## [34,] 0.3536728 -1.001925e-04 -0.06922437 -0.02000618 -1.001925e-04
## [35,] 0.3640902 -1.043390e-04 -0.07082861 -0.02073590 -1.043390e-04
## [36,] 0.3746576 -1.085808e-04 -0.07244950 -0.02147282 -1.085808e-04
## [37,] 0.3849375 -1.127330e-04 -0.07401208 -0.02219670 -1.127330e-04
## [38,] 0.3947010 -1.167298e-04 -0.07547263 -0.02289547 -1.167298e-04
## [39,] 0.4040562 -1.206876e-04 -0.07683912 -0.02357128 -1.206876e-04
## [40,] 0.4133473 -1.248237e-04 -0.07815926 -0.02423718 -1.248237e-04
## [41,] 0.4229433 -1.293353e-04 -0.07948902 -0.02490911 -1.293353e-04
## [42,] 0.4330542 -1.343053e-04 -0.08086364 -0.02559794 -1.343053e-04
## [43,] 0.4437043 -1.396955e-04 -0.08229388 -0.02630673 -1.396955e-04
## [44,] 0.4548597 -1.454288e-04 -0.08378474 -0.02703377 -1.454288e-04
## [45,] 0.4665259 -1.514596e-04 -0.08534773 -0.02777641 -1.514596e-04
## [46,] 0.4787182 -1.577612e-04 -0.08699472 -0.02853092 -1.577612e-04
## [47,] 0.4913971 -1.642794e-04 -0.08872961 -0.02929039 -1.642794e-04
## [48,] 0.5044951 -1.709336e-04 -0.09055583 -0.03004461 -1.709336e-04
## [49,] 0.5180468 -1.776809e-04 -0.09249866 -0.03078299 -1.776809e-04
## [50,] 0.5323135 -1.845829e-04 -0.09462602 -0.03149723 -1.845829e-04
## [51,] 0.5478012 -1.918176e-04 -0.09705170 -0.03218073 -1.918176e-04
## [52,] 0.5651172 -1.996195e-04 -0.09990860 -0.03282554 -1.996195e-04
## [53,] 0.5846322 -2.081456e-04 -0.10328101 -0.03341891 -2.081456e-04
## [54,] 0.6060941 -2.173193e-04 -0.10712792 -0.03393973 -2.173193e-04
## [55,] 0.6286032 -2.268089e-04 -0.11128295 -0.03435524 -2.268089e-04
## [56,] 0.6509906 -2.361793e-04 -0.11553249 -0.03462146 -2.361793e-04
## [57,] 0.6721873 -2.450426e-04 -0.11968808 -0.03468619 -2.450426e-04
## [58,] 0.6915003 -2.531413e-04 -0.12364640 -0.03448975 -2.531413e-04
## [59,] 0.7087901 -2.603755e-04 -0.12743379 -0.03396329 -2.603755e-04
## [60,] 0.7243050 -2.666693e-04 -0.13118469 -0.03303003 -2.666693e-04
## [61,] 0.7381337 -2.716274e-04 -0.13504765 -0.03161654 -2.716274e-04
## [62,] 0.7496914 -2.742269e-04 -0.13908569 -0.02968022 -2.742269e-04
## [63,] 0.7578838 -2.729691e-04 -0.14327045 -0.02725453 -2.729691e-04
## [64,] 0.7622446 -2.666521e-04 -0.14762132 -0.02451291 -2.666521e-04
## [65,] 0.7642818 -2.553766e-04 -0.15235193 -0.02185054 -2.553766e-04
## [66,] 0.7676326 -2.411124e-04 -0.15779853 -0.01990419 -2.411124e-04
## [67,] 0.7762047 -2.272355e-04 -0.16409808 -0.01930341 -2.272355e-04
## [68,] 0.7903814 -2.166379e-04 -0.17070155 -0.02012179 -2.166379e-04
## [69,] 0.8057094 -2.099789e-04 -0.17641568 -0.02164941 -2.099789e-04
## [70,] 0.8175817 -2.063126e-04 -0.18039092 -0.02304308 -2.063126e-04
## [71,] 0.8251571 -2.045245e-04 -0.18275980 -0.02397204 -2.045245e-04
## [72,] 0.8297302 -2.037897e-04 -0.18409794 -0.02450769 -2.037897e-04
## [73,] 0.8324690 -2.035545e-04 -0.18485074 -0.02479995 -2.035545e-04
## [74,] 0.8340790 -2.035114e-04 -0.18527220 -0.02495630 -2.035114e-04
## [75,] 0.8349999 -2.035247e-04 -0.18550518 -0.02503928 -2.035247e-04
## [76,] 0.8355135 -2.035462e-04 -0.18563218 -0.02508313 -2.035462e-04
## [77,] 0.8357944 -2.035632e-04 -0.18570061 -0.02510625 -2.035632e-04
## [78,] 0.8359460 -2.035741e-04 -0.18573716 -0.02511843 -2.035741e-04
## [79,] 0.8360271 -2.035806e-04 -0.18575658 -0.02512483 -2.035806e-04
## [80,] 0.8360702 -2.035843e-04 -0.18576686 -0.02512820 -2.035843e-04
## [81,] 0.8360930 -2.035864e-04 -0.18577228 -0.02512998 -2.035864e-04
## [82,] 0.8361050 -2.035875e-04 -0.18577514 -0.02513091 -2.035875e-04
## [83,] 0.8361114 -2.035881e-04 -0.18577664 -0.02513140 -2.035881e-04
## [84,] 0.8361147 -2.035884e-04 -0.18577743 -0.02513166 -2.035884e-04
## [85,] 0.8361165 -2.035886e-04 -0.18577785 -0.02513179 -2.035886e-04
## [86,] 0.8361174 -2.035887e-04 -0.18577807 -0.02513186 -2.035887e-04
## [87,] 0.8361179 -2.035887e-04 -0.18577819 -0.02513190 -2.035887e-04
## [88,] 0.8361181 -2.035887e-04 -0.18577825 -0.02513192 -2.035887e-04
## [89,] 0.8361183 -2.035887e-04 -0.18577828 -0.02513193 -2.035887e-04
## [90,] 0.8361184 -2.035888e-04 -0.18577830 -0.02513194 -2.035888e-04
## [91,] 0.8361184 -2.035888e-04 -0.18577830 -0.02513194 -2.035888e-04
## [92,] 0.8361184 -2.035888e-04 -0.18577831 -0.02513194 -2.035888e-04
## [93,] 0.8361184 -2.035888e-04 -0.18577831 -0.02513194 -2.035888e-04
## [94,] 0.8361184 -2.035888e-04 -0.18577831 -0.02513194 -2.035888e-04
## [95,] 0.8361184 -2.035888e-04 -0.18577831 -0.02513194 -2.035888e-04
## [96,] 0.8361184 -2.035888e-04 -0.18577831 -0.02513194 -2.035888e-04
## [97,] 0.8361184 -2.035888e-04 -0.18577831 -0.02513194 -2.035888e-04
## [98,] 0.8361184 -2.035888e-04 -0.18577831 -0.02513194 -2.035888e-04
## [99,] 0.8361184 -2.035888e-04 -0.18577831 -0.02513194 -2.035888e-04
## [100,] 0.8361184 -2.035888e-04 -0.18577831 -0.02513194 -2.035888e-04
## [,6] [,7] [,8] [,9] [,10]
## [1,] 1.850942e-05 -2.182722e-03 2.403031e-04 -4.61036173 -2.182722e-03
## [2,] 2.290188e-07 -2.536181e-05 3.478401e-06 -0.05453021 -2.536181e-05
## [3,] 2.300583e-07 -2.530629e-05 3.630541e-06 -0.05475095 -2.530629e-05
## [4,] 2.311326e-07 -2.525155e-05 3.797484e-06 -0.05503809 -2.525155e-05
## [5,] 2.322237e-07 -2.520570e-05 3.975452e-06 -0.05537671 -2.520570e-05
## [6,] 2.333054e-07 -2.517818e-05 4.158741e-06 -0.05572879 -2.517818e-05
## [7,] 2.343517e-07 -2.517649e-05 4.341821e-06 -0.05605436 -2.517649e-05
## [8,] 2.353481e-07 -2.520206e-05 4.522444e-06 -0.05633835 -2.520206e-05
## [9,] 2.362983e-07 -2.525124e-05 4.702449e-06 -0.05659171 -2.525124e-05
## [10,] 2.372180e-07 -2.532016e-05 4.884889e-06 -0.05683248 -2.532016e-05
## [11,] 2.381205e-07 -2.540779e-05 5.070245e-06 -0.05707032 -2.540779e-05
## [12,] 2.390084e-07 -2.551532e-05 5.254993e-06 -0.05730464 -2.551532e-05
## [13,] 2.398756e-07 -2.564398e-05 5.433156e-06 -0.05753257 -2.564398e-05
## [14,] 2.407177e-07 -2.579403e-05 5.598984e-06 -0.05775699 -2.579403e-05
## [15,] 2.415410e-07 -2.596491e-05 5.748876e-06 -0.05798834 -2.596491e-05
## [16,] 2.423670e-07 -2.615486e-05 5.882208e-06 -0.05823994 -2.615486e-05
## [17,] 2.432320e-07 -2.635871e-05 6.001623e-06 -0.05852022 -2.635871e-05
## [18,] 2.441827e-07 -2.656556e-05 6.112928e-06 -0.05882626 -2.656556e-05
## [19,] 2.452631e-07 -2.676018e-05 6.223652e-06 -0.05914373 -2.676018e-05
## [20,] 2.464876e-07 -2.692975e-05 6.339475e-06 -0.05945341 -2.692975e-05
## [21,] 2.478188e-07 -2.707075e-05 6.460784e-06 -0.05973964 -2.707075e-05
## [22,] 2.491941e-07 -2.718826e-05 6.583951e-06 -0.05999708 -2.718826e-05
## [23,] 2.505877e-07 -2.728885e-05 6.706471e-06 -0.06023490 -2.728885e-05
## [24,] 2.520256e-07 -2.737565e-05 6.829284e-06 -0.06047591 -2.737565e-05
## [25,] 2.535598e-07 -2.744834e-05 6.955889e-06 -0.06075259 -2.744834e-05
## [26,] 2.552665e-07 -2.750427e-05 7.092695e-06 -0.06110356 -2.750427e-05
## [27,] 2.572548e-07 -2.754015e-05 7.249697e-06 -0.06156863 -2.754015e-05
## [28,] 2.596384e-07 -2.755512e-05 7.437781e-06 -0.06218139 -2.755512e-05
## [29,] 2.624801e-07 -2.755245e-05 7.663723e-06 -0.06296273 -2.755245e-05
## [30,] 2.657632e-07 -2.753889e-05 7.927060e-06 -0.06391761 -2.753889e-05
## [31,] 2.694170e-07 -2.752243e-05 8.221240e-06 -0.06503617 -2.752243e-05
## [32,] 2.733778e-07 -2.750950e-05 8.537735e-06 -0.06630248 -2.750950e-05
## [33,] 2.776351e-07 -2.750211e-05 8.869874e-06 -0.06770357 -2.750211e-05
## [34,] 2.822156e-07 -2.749703e-05 9.214176e-06 -0.06922437 -2.749703e-05
## [35,] 2.870849e-07 -2.748910e-05 9.568122e-06 -0.07082861 -2.748910e-05
## [36,] 2.920529e-07 -2.747686e-05 9.926326e-06 -0.07244950 -2.747686e-05
## [37,] 2.968647e-07 -2.746155e-05 1.028134e-05 -0.07401208 -2.746155e-05
## [38,] 3.014082e-07 -2.743889e-05 1.062934e-05 -0.07547263 -2.743889e-05
## [39,] 3.057775e-07 -2.739486e-05 1.097291e-05 -0.07683912 -2.739486e-05
## [40,] 3.101796e-07 -2.731009e-05 1.131857e-05 -0.07815926 -2.731009e-05
## [41,] 3.148132e-07 -2.716945e-05 1.167268e-05 -0.07948902 -2.716945e-05
## [42,] 3.197839e-07 -2.697069e-05 1.203790e-05 -0.08086364 -2.697069e-05
## [43,] 3.250906e-07 -2.672327e-05 1.241246e-05 -0.08229388 -2.672327e-05
## [44,] 3.306851e-07 -2.643767e-05 1.279269e-05 -0.08378474 -2.643767e-05
## [45,] 3.365326e-07 -2.611861e-05 1.317570e-05 -0.08534773 -2.611861e-05
## [46,] 3.426131e-07 -2.576865e-05 1.355941e-05 -0.08699472 -2.576865e-05
## [47,] 3.488909e-07 -2.539439e-05 1.394061e-05 -0.08872961 -2.539439e-05
## [48,] 3.553162e-07 -2.500579e-05 1.431427e-05 -0.09055583 -2.500579e-05
## [49,] 3.618694e-07 -2.460724e-05 1.467475e-05 -0.09249866 -2.460724e-05
## [50,] 3.686047e-07 -2.418791e-05 1.501725e-05 -0.09462602 -2.418791e-05
## [51,] 3.756511e-07 -2.371826e-05 1.533745e-05 -0.09705170 -2.371826e-05
## [52,] 3.831609e-07 -2.315610e-05 1.562972e-05 -0.09990860 -2.315610e-05
## [53,] 3.912142e-07 -2.246543e-05 1.588537e-05 -0.10328101 -2.246543e-05
## [54,] 3.997252e-07 -2.164209e-05 1.609122e-05 -0.10712792 -2.164209e-05
## [55,] 4.084639e-07 -2.072295e-05 1.622852e-05 -0.11128295 -2.072295e-05
## [56,] 4.172135e-07 -1.977225e-05 1.627312e-05 -0.11553249 -1.977225e-05
## [57,] 4.259202e-07 -1.886786e-05 1.619623e-05 -0.11968808 -1.886786e-05
## [58,] 4.347632e-07 -1.809100e-05 1.596352e-05 -0.12364640 -1.809100e-05
## [59,] 4.441694e-07 -1.752008e-05 1.553291e-05 -0.12743379 -1.752008e-05
## [60,] 4.547289e-07 -1.725183e-05 1.485421e-05 -0.13118469 -1.725183e-05
## [61,] 4.670016e-07 -1.746349e-05 1.387401e-05 -0.13504765 -1.746349e-05
## [62,] 4.814440e-07 -1.849026e-05 1.254694e-05 -0.13908569 -1.849026e-05
## [63,] 4.987526e-07 -2.085709e-05 1.085664e-05 -0.14327045 -2.085709e-05
## [64,] 5.204011e-07 -2.518726e-05 8.864372e-06 -0.14762132 -2.518726e-05
## [65,] 5.484120e-07 -3.192786e-05 6.805690e-06 -0.15235193 -3.192786e-05
## [66,] 5.838113e-07 -4.093527e-05 5.153655e-06 -0.15779853 -4.093527e-05
## [67,] 6.245716e-07 -5.113168e-05 4.407305e-06 -0.16409808 -5.113168e-05
## [68,] 6.649320e-07 -6.073590e-05 4.632175e-06 -0.17070155 -6.073590e-05
## [69,] 6.979010e-07 -6.820199e-05 5.365730e-06 -0.17641568 -6.820199e-05
## [70,] 7.199238e-07 -7.303124e-05 6.087846e-06 -0.18039092 -7.303124e-05
## [71,] 7.326195e-07 -7.574779e-05 6.581494e-06 -0.18275980 -7.574779e-05
## [72,] 7.395312e-07 -7.717350e-05 6.866039e-06 -0.18409794 -7.717350e-05
## [73,] 7.432771e-07 -7.790827e-05 7.018697e-06 -0.18485074 -7.790827e-05
## [74,] 7.453120e-07 -7.828736e-05 7.098542e-06 -0.18527220 -7.828736e-05
## [75,] 7.464137e-07 -7.848381e-05 7.140014e-06 -0.18550518 -7.848381e-05
## [76,] 7.470063e-07 -7.858601e-05 7.161551e-06 -0.18563218 -7.858601e-05
## [77,] 7.473230e-07 -7.863932e-05 7.172755e-06 -0.18570061 -7.863932e-05
## [78,] 7.474913e-07 -7.866720e-05 7.178598e-06 -0.18573716 -7.866720e-05
## [79,] 7.475805e-07 -7.868179e-05 7.181650e-06 -0.18575658 -7.868179e-05
## [80,] 7.476277e-07 -7.868944e-05 7.183248e-06 -0.18576686 -7.868944e-05
## [81,] 7.476525e-07 -7.869346e-05 7.184085e-06 -0.18577228 -7.869346e-05
## [82,] 7.476656e-07 -7.869556e-05 7.184524e-06 -0.18577514 -7.869556e-05
## [83,] 7.476725e-07 -7.869667e-05 7.184755e-06 -0.18577664 -7.869667e-05
## [84,] 7.476762e-07 -7.869725e-05 7.184876e-06 -0.18577743 -7.869725e-05
## [85,] 7.476781e-07 -7.869756e-05 7.184939e-06 -0.18577785 -7.869756e-05
## [86,] 7.476791e-07 -7.869772e-05 7.184973e-06 -0.18577807 -7.869772e-05
## [87,] 7.476796e-07 -7.869780e-05 7.184990e-06 -0.18577819 -7.869780e-05
## [88,] 7.476799e-07 -7.869785e-05 7.185000e-06 -0.18577825 -7.869785e-05
## [89,] 7.476800e-07 -7.869787e-05 7.185004e-06 -0.18577828 -7.869787e-05
## [90,] 7.476801e-07 -7.869788e-05 7.185007e-06 -0.18577830 -7.869788e-05
## [91,] 7.476802e-07 -7.869789e-05 7.185008e-06 -0.18577830 -7.869789e-05
## [92,] 7.476802e-07 -7.869789e-05 7.185009e-06 -0.18577831 -7.869789e-05
## [93,] 7.476802e-07 -7.869790e-05 7.185009e-06 -0.18577831 -7.869790e-05
## [94,] 7.476802e-07 -7.869790e-05 7.185010e-06 -0.18577831 -7.869790e-05
## [95,] 7.476802e-07 -7.869790e-05 7.185010e-06 -0.18577831 -7.869790e-05
## [96,] 7.476802e-07 -7.869790e-05 7.185010e-06 -0.18577831 -7.869790e-05
## [97,] 7.476802e-07 -7.869790e-05 7.185010e-06 -0.18577831 -7.869790e-05
## [98,] 7.476802e-07 -7.869790e-05 7.185010e-06 -0.18577831 -7.869790e-05
## [99,] 7.476802e-07 -7.869790e-05 7.185010e-06 -0.18577831 -7.869790e-05
## [100,] 7.476802e-07 -7.869790e-05 7.185010e-06 -0.18577831 -7.869790e-05
## [,11] [,12] [,13] [,14] [,15]
## [1,] 1.69493752 0.0229414838 -0.81335512 2.403031e-04 0.0229414838
## [2,] 0.02038305 0.0001271585 -0.00983559 3.478401e-06 0.0001271585
## [3,] 0.02042402 0.0001434648 -0.01005653 3.630541e-06 0.0001434648
## [4,] 0.02048113 0.0001637805 -0.01031425 3.797484e-06 0.0001637805
## [5,] 0.02055270 0.0001866301 -0.01060006 3.975452e-06 0.0001866301
## [6,] 0.02063183 0.0002090064 -0.01089698 4.158741e-06 0.0002090064
## [7,] 0.02071067 0.0002279512 -0.01118849 4.341821e-06 0.0002279512
## [8,] 0.02078585 0.0002425391 -0.01147014 4.522444e-06 0.0002425391
## [9,] 0.02085904 0.0002537960 -0.01174986 4.702449e-06 0.0002537960
## [10,] 0.02093375 0.0002630108 -0.01203710 4.884889e-06 0.0002630108
## [11,] 0.02101240 0.0002704676 -0.01233289 5.070245e-06 0.0002704676
## [12,] 0.02109560 0.0002755487 -0.01262873 5.254993e-06 0.0002755487
## [13,] 0.02118318 0.0002777094 -0.01291329 5.433156e-06 0.0002777094
## [14,] 0.02127577 0.0002772642 -0.01317962 5.598984e-06 0.0002772642
## [15,] 0.02137525 0.0002754036 -0.01342771 5.748876e-06 0.0002754036
## [16,] 0.02148391 0.0002736246 -0.01366251 5.882208e-06 0.0002736246
## [17,] 0.02160233 0.0002730081 -0.01389041 6.001623e-06 0.0002730081
## [18,] 0.02172761 0.0002737174 -0.01411624 6.112928e-06 0.0002737174
## [19,] 0.02185354 0.0002750124 -0.01434197 6.223652e-06 0.0002750124
## [20,] 0.02197323 0.0002757769 -0.01456628 6.339475e-06 0.0002757769
## [21,] 0.02208223 0.0002751566 -0.01478495 6.460784e-06 0.0002751566
## [22,] 0.02217998 0.0002730297 -0.01499411 6.583951e-06 0.0002730297
## [23,] 0.02226953 0.0002702377 -0.01519522 6.706471e-06 0.0002702377
## [24,] 0.02235655 0.0002684583 -0.01539626 6.829284e-06 0.0002684583
## [25,] 0.02244845 0.0002698780 -0.01560967 6.955889e-06 0.0002698780
## [26,] 0.02255378 0.0002769241 -0.01585156 7.092695e-06 0.0002769241
## [27,] 0.02268136 0.0002919569 -0.01614093 7.249697e-06 0.0002919569
## [28,] 0.02283907 0.0003168162 -0.01649552 7.437781e-06 0.0003168162
## [29,] 0.02303265 0.0003524544 -0.01692584 7.663723e-06 0.0003524544
## [30,] 0.02326469 0.0003988108 -0.01743188 7.927060e-06 0.0003988108
## [31,] 0.02353449 0.0004549668 -0.01800452 8.221240e-06 0.0004549668
## [32,] 0.02383964 0.0005196742 -0.01863116 8.537735e-06 0.0005196742
## [33,] 0.02417760 0.0005918806 -0.01930110 8.869874e-06 0.0005918806
## [34,] 0.02454432 0.0006705534 -0.02000618 9.214176e-06 0.0006705534
## [35,] 0.02493012 0.0007538041 -0.02073590 9.568122e-06 0.0007538041
## [36,] 0.02531837 0.0008384151 -0.02147282 9.926326e-06 0.0008384151
## [37,] 0.02569096 0.0009207734 -0.02219670 1.028134e-05 0.0009207734
## [38,] 0.02603649 0.0009986381 -0.02289547 1.062934e-05 0.0009986381
## [39,] 0.02635455 0.0010721636 -0.02357128 1.097291e-05 0.0010721636
## [40,] 0.02665335 0.0011435193 -0.02423718 1.131857e-05 0.0011435193
## [41,] 0.02694360 0.0012156116 -0.02490911 1.167268e-05 0.0012156116
## [42,] 0.02723286 0.0012907635 -0.02559794 1.203790e-05 0.0012907635
## [43,] 0.02752478 0.0013701783 -0.02630673 1.241246e-05 0.0013701783
## [44,] 0.02782222 0.0014542981 -0.02703377 1.279269e-05 0.0014542981
## [45,] 0.02812914 0.0015432401 -0.02777641 1.317570e-05 0.0015432401
## [46,] 0.02844951 0.0016367590 -0.02853092 1.355941e-05 0.0016367590
## [47,] 0.02878628 0.0017340343 -0.02929039 1.394061e-05 0.0017340343
## [48,] 0.02914316 0.0018338016 -0.03004461 1.431427e-05 0.0018338016
## [49,] 0.02952889 0.0019348751 -0.03078299 1.467475e-05 0.0019348751
## [50,] 0.02996092 0.0020365785 -0.03149723 1.501725e-05 0.0020365785
## [51,] 0.03046574 0.0021386968 -0.03218073 1.533745e-05 0.0021386968
## [52,] 0.03107281 0.0022410584 -0.03282554 1.562972e-05 0.0022410584
## [53,] 0.03179941 0.0023430066 -0.03341891 1.588537e-05 0.0023430066
## [54,] 0.03263450 0.0024426497 -0.03393973 1.609122e-05 0.0024426497
## [55,] 0.03354123 0.0025359580 -0.03435524 1.622852e-05 0.0025359580
## [56,] 0.03447653 0.0026165100 -0.03462146 1.627312e-05 0.0026165100
## [57,] 0.03540847 0.0026758165 -0.03468619 1.619623e-05 0.0026758165
## [58,] 0.03633085 0.0027032667 -0.03448975 1.596352e-05 0.0027032667
## [59,] 0.03727450 0.0026856945 -0.03396329 1.553291e-05 0.0026856945
## [60,] 0.03830530 0.0026073242 -0.03303003 1.485421e-05 0.0026073242
## [61,] 0.03950890 0.0024510654 -0.03161654 1.387401e-05 0.0024510654
## [62,] 0.04097574 0.0022019611 -0.02968022 1.254694e-05 0.0022019611
## [63,] 0.04280259 0.0018524761 -0.02725453 1.085664e-05 0.0018524761
## [64,] 0.04511076 0.0014086001 -0.02451291 8.864372e-06 0.0014086001
## [65,] 0.04803065 0.0008981730 -0.02185054 6.805690e-06 0.0008981730
## [66,] 0.05159855 0.0003803998 -0.01990419 5.153655e-06 0.0003803998
## [67,] 0.05560123 -0.0000537485 -0.01930341 4.407305e-06 -0.0000537485
## [68,] 0.05950089 -0.0003205564 -0.02012179 4.632175e-06 -0.0003205564
## [69,] 0.06265682 -0.0004168764 -0.02164941 5.365730e-06 -0.0004168764
## [70,] 0.06475747 -0.0004146944 -0.02304308 6.087846e-06 -0.0004146944
## [71,] 0.06596819 -0.0003836240 -0.02397204 6.581494e-06 -0.0003836240
## [72,] 0.06662646 -0.0003552504 -0.02450769 6.866039e-06 -0.0003552504
## [73,] 0.06698166 -0.0003361819 -0.02479995 7.018697e-06 -0.0003361819
## [74,] 0.06717344 -0.0003247938 -0.02495630 7.098542e-06 -0.0003247938
## [75,] 0.06727661 -0.0003183608 -0.02503928 7.140014e-06 -0.0003183608
## [76,] 0.06733179 -0.0003148313 -0.02508313 7.161551e-06 -0.0003148313
## [77,] 0.06736115 -0.0003129260 -0.02510625 7.172755e-06 -0.0003129260
## [78,] 0.06737671 -0.0003119074 -0.02511843 7.178598e-06 -0.0003119074
## [79,] 0.06738493 -0.0003113660 -0.02512483 7.181650e-06 -0.0003113660
## [80,] 0.06738926 -0.0003110792 -0.02512820 7.183248e-06 -0.0003110792
## [81,] 0.06739154 -0.0003109277 -0.02512998 7.184085e-06 -0.0003109277
## [82,] 0.06739275 -0.0003108478 -0.02513091 7.184524e-06 -0.0003108478
## [83,] 0.06739338 -0.0003108057 -0.02513140 7.184755e-06 -0.0003108057
## [84,] 0.06739371 -0.0003107835 -0.02513166 7.184876e-06 -0.0003107835
## [85,] 0.06739388 -0.0003107718 -0.02513179 7.184939e-06 -0.0003107718
## [86,] 0.06739398 -0.0003107657 -0.02513186 7.184973e-06 -0.0003107657
## [87,] 0.06739402 -0.0003107625 -0.02513190 7.184990e-06 -0.0003107625
## [88,] 0.06739405 -0.0003107608 -0.02513192 7.185000e-06 -0.0003107608
## [89,] 0.06739406 -0.0003107599 -0.02513193 7.185004e-06 -0.0003107599
## [90,] 0.06739407 -0.0003107594 -0.02513194 7.185007e-06 -0.0003107594
## [91,] 0.06739407 -0.0003107591 -0.02513194 7.185008e-06 -0.0003107591
## [92,] 0.06739408 -0.0003107590 -0.02513194 7.185009e-06 -0.0003107590
## [93,] 0.06739408 -0.0003107590 -0.02513194 7.185009e-06 -0.0003107590
## [94,] 0.06739408 -0.0003107589 -0.02513194 7.185010e-06 -0.0003107589
## [95,] 0.06739408 -0.0003107589 -0.02513194 7.185010e-06 -0.0003107589
## [96,] 0.06739408 -0.0003107589 -0.02513194 7.185010e-06 -0.0003107589
## [97,] 0.06739408 -0.0003107589 -0.02513194 7.185010e-06 -0.0003107589
## [98,] 0.06739408 -0.0003107589 -0.02513194 7.185010e-06 -0.0003107589
## [99,] 0.06739408 -0.0003107589 -0.02513194 7.185010e-06 -0.0003107589
## [100,] 0.06739408 -0.0003107589 -0.02513194 7.185010e-06 -0.0003107589
## [,16]
## [1,] 0.341336070
## [2,] 0.002943093
## [3,] 0.002968335
## [4,] 0.002998336
## [5,] 0.003032607
## [6,] 0.003069995
## [7,] 0.003109424
## [8,] 0.003150890
## [9,] 0.003195489
## [10,] 0.003244307
## [11,] 0.003297273
## [12,] 0.003352967
## [13,] 0.003409432
## [14,] 0.003465146
## [15,] 0.003519469
## [16,] 0.003572475
## [17,] 0.003624530
## [18,] 0.003675902
## [19,] 0.003726584
## [20,] 0.003776280
## [21,] 0.003824495
## [22,] 0.003870812
## [23,] 0.003915248
## [24,] 0.003958330
## [25,] 0.004000950
## [26,] 0.004044333
## [27,] 0.004090003
## [28,] 0.004139461
## [29,] 0.004193678
## [30,] 0.004252771
## [31,] 0.004316082
## [32,] 0.004382586
## [33,] 0.004451294
## [34,] 0.004521396
## [35,] 0.004592109
## [36,] 0.004662482
## [37,] 0.004731521
## [38,] 0.004798597
## [39,] 0.004863688
## [40,] 0.004927241
## [41,] 0.004989823
## [42,] 0.005051775
## [43,] 0.005113064
## [44,] 0.005173374
## [45,] 0.005232274
## [46,] 0.005289249
## [47,] 0.005343634
## [48,] 0.005394592
## [49,] 0.005441188
## [50,] 0.005482465
## [51,] 0.005517404
## [52,] 0.005544778
## [53,] 0.005563015
## [54,] 0.005570229
## [55,] 0.005564455
## [56,] 0.005544068
## [57,] 0.005508333
## [58,] 0.005458033
## [59,] 0.005396183
## [60,] 0.005329041
## [61,] 0.005268078
## [62,] 0.005233990
## [63,] 0.005263199
## [64,] 0.005413782
## [65,] 0.005761703
## [66,] 0.006374739
## [67,] 0.007250778
## [68,] 0.008241222
## [69,] 0.009087222
## [70,] 0.009626272
## [71,] 0.009898778
## [72,] 0.010021880
## [73,] 0.010076928
## [74,] 0.010102416
## [75,] 0.010114693
## [76,] 0.010120789
## [77,] 0.010123879
## [78,] 0.010125465
## [79,] 0.010126287
## [80,] 0.010126715
## [81,] 0.010126938
## [82,] 0.010127055
## [83,] 0.010127117
## [84,] 0.010127149
## [85,] 0.010127166
## [86,] 0.010127174
## [87,] 0.010127179
## [88,] 0.010127182
## [89,] 0.010127183
## [90,] 0.010127184
## [91,] 0.010127184
## [92,] 0.010127184
## [93,] 0.010127184
## [94,] 0.010127184
## [95,] 0.010127184
## [96,] 0.010127184
## [97,] 0.010127184
## [98,] 0.010127184
## [99,] 0.010127184
## [100,] 0.010127184
result$beta[iter,] # Final Estimates
## [1] -3.44954840 0.00229396 0.77701357 -0.56003139
matrix(result$cov[iter,],4,4) # Final Covariance Estimates
## [,1] [,2] [,3] [,4]
## [1,] 0.8361184299 -2.035888e-04 -0.1857783138 -2.513194e-02
## [2,] -0.0002035888 7.476802e-07 -0.0000786979 7.185010e-06
## [3,] -0.1857783138 -7.869790e-05 0.0673940779 -3.107589e-04
## [4,] -0.0251319410 7.185010e-06 -0.0003107589 1.012718e-02
############### 3. EM Algorithm ###############
# Read the data
data = as.matrix(scan("mixture data.txt"))
class(data)
## [1] "matrix"
head(data,10)
## [,1]
## [1,] -1.6273898
## [2,] -1.4238937
## [3,] -0.6787816
## [4,] -2.2420957
## [5,] -1.8670102
## [6,] -1.9708137
## [7,] -2.0132615
## [8,] -0.1043120
## [9,] -2.4802637
## [10,] -1.2142840
n = length(data)
iter = 30; mu1 = -1.5; mu2 = 1.5; sigma = 0.5; # Parameters
delta = c(); delta[1] = 0.5 # Initial values
EU <- c()
# Expectation-Maximization Algorithm
for(k in 1:iter){
for(i in 1:n){
EU[i] <- delta[k] * dnorm(data[i],mu1,sigma)/((1-delta[k]) * dnorm(data[i],mu2,sigma) + delta[k] * dnorm(data[i],mu1,sigma))
}
delta[k+1] = sum(EU)/n
}
delta <- as.matrix(delta)
delta
## [,1]
## [1,] 0.5000000
## [2,] 0.7480462
## [3,] 0.7499138
## [4,] 0.7499282
## [5,] 0.7499283
## [6,] 0.7499283
## [7,] 0.7499283
## [8,] 0.7499283
## [9,] 0.7499283
## [10,] 0.7499283
## [11,] 0.7499283
## [12,] 0.7499283
## [13,] 0.7499283
## [14,] 0.7499283
## [15,] 0.7499283
## [16,] 0.7499283
## [17,] 0.7499283
## [18,] 0.7499283
## [19,] 0.7499283
## [20,] 0.7499283
## [21,] 0.7499283
## [22,] 0.7499283
## [23,] 0.7499283
## [24,] 0.7499283
## [25,] 0.7499283
## [26,] 0.7499283
## [27,] 0.7499283
## [28,] 0.7499283
## [29,] 0.7499283
## [30,] 0.7499283
## [31,] 0.7499283
plot(delta, type ="l", main="", xlab="EM Iteration for the Mixture Data")
############### 4. Metropolis-Hasting Algorithm ###############
n = 10500; burnin = 1:500; # Burn-in
s = 0.01; x.init1 = -3; x.init2 = -1.5; x.init3 = 2; # Starting points
f = function(x){0.6*dnorm(x,-1.5,0.5)+0.4*dnorm(x,1.5,0.5)}
x = seq(-4,4,by=0.01);
plot(x,f(x),type="l")
MCMC = function(x.init,s,n,f){
x = NULL; x[1] = x.init;
for(i in 1:(n-1)){
x[i+1] = rnorm(1,x[i],s)
MH.ratio = f(x[i+1])/f(x[i])
u = runif(1)
if(u > MH.ratio){x[i+1] = x[i]}
}
return(x)
}
x.set1 = MCMC(x.init1, s, n, f)
x.set2 = MCMC(x.init2, s, n, f)
x.set3 = MCMC(x.init3, s, n, f)
mean(x.set1[-burnin])
## [1] -2.15861
plot(x.set1[-burnin],type="l",ylab="x",xlab="t")
acf(x.set1)
par(mfrow=c(1,1))
z = seq(-4,4,by=0.01)
hist(x.set1[-burnin],breaks=50,freq=FALSE,xlab="x",main="Histogram of Mixture Data",ylab="Density")
points(z,f(z),type="l")
abline(v=-1.5);abline(v=1.5) # Highlight the approximate modes of the true distribution
mean(x.set2[-burnin])
## [1] -1.486066
plot(x.set2[-burnin],type="l",ylab="x",xlab="t")
acf(x.set2)
par(mfrow=c(1,1))
z = seq(-4,4,by=0.01)
hist(x.set2[-burnin],breaks=50,freq=FALSE,xlab="x",main="Histogram of Mixture Data",ylab="Density")
points(z,f(z),type="l")
abline(v=-1.5);abline(v=1.5) # Highlight the approximate modes of the true distribution
mean(x.set3[-burnin])
## [1] 2.33435
plot(x.set3[-burnin],type="l",ylab="x",xlab="t")
acf(x.set3)
par(mfrow=c(1,1))
z = seq(-4,4,by=0.01)
hist(x.set3[-burnin],breaks=50,freq=FALSE,xlab="x",main="Histogram of Mixture Data",ylab="Density")
points(z,f(z),type="l")
abline(v=-1.5);abline(v=1.5) # Highlight the approximate modes of the true distribution
# Reset the variance
s = 1.5
x.set1 = MCMC(x.init1, s, n, f)
x.set2 = MCMC(x.init2, s, n, f)
x.set3 = MCMC(x.init3, s, n, f)
mean(x.set1[-burnin])
## [1] -0.354797
plot(x.set1[-burnin],type="l",ylab="x",xlab="t")
acf(x.set1)
par(mfrow=c(1,1))
z = seq(-4,4,by=0.01)
hist(x.set1[-burnin],breaks=50,freq=FALSE,xlab="x",main="Histogram of Mixture Data",ylab="Density")
points(z,f(z),type="l")
abline(v=-1.5);abline(v=1.5) # Highlight the approximate modes of the true distribution
mean(x.set2[-burnin])
## [1] -0.3337031
plot(x.set2[-burnin],type="l",ylab="x",xlab="t")
acf(x.set2)
par(mfrow=c(1,1))
z = seq(-4,4,by=0.01)
hist(x.set2[-burnin],breaks=50,freq=FALSE,xlab="x",main="Histogram of Mixture Data",ylab="Density")
points(z,f(z),type="l")
abline(v=-1.5);abline(v=1.5) # Highlight the approximate modes of the true distribution
mean(x.set3[-burnin])
## [1] -0.26175
plot(x.set3[-burnin],type="l",ylab="x",xlab="t")
acf(x.set3)
par(mfrow=c(1,1))
z = seq(-4,4,by=0.01)
hist(x.set3[-burnin],breaks=50,freq=FALSE,xlab="x",main="Histogram of Mixture Data",ylab="Density")
points(z,f(z),type="l")
abline(v=-1.5);abline(v=1.5) # Highlight the approximate modes of the true distribution
############### 5. Gibbs Sampler ###############
# Read in data
data = read.table("Gibbs.txt", header = TRUE)
Y <- data$Y
x <- data$X
n = length(x)
# Initialize constants and parameters
set.seed(50)
N <- 12000 # length of chain
burn <- 2000 # burn-in length
taua = 0.001; taub = 0.001; alpha0 = 0.001; beta0 = 0.001;
alpha <- rep(0, N); beta <- rep(0, N); gamma <- rep(0, N);
alpha[1] <- rnorm(1,0,sqrt(1/taua))
beta[1] <- rnorm(1,0,sqrt(1/taub))
gamma[1] <- rgamma(1,shape=alpha0,scale=beta0)
for (i in 2:N) {
a1 = c(); a2 = c(); a3 = c(); a4 = c() # Compute the sums
for(j in 1:n){
a1[j] = Y[j]-beta[i]*x[j]
a2[j] = (x[j])^2
a3[j] = (Y[j]-alpha[i])*x[j]
a4[j] = (Y[j]-alpha[i]-beta[i]*x[j])^2
}
# sample the alpha direction conditional on beta and gamma
alpha[i] <- rnorm(1,gamma[i]*sum(a1)/(n*gamma[i]+taua),1/(n*gamma[i]+taua))
# sample the beta direction conditional on alpha and gamma
beta[i] <- rnorm(1,gamma[i]*sum(a3)/(gamma[i]*sum(a2)+taub),1/(gamma[i]*sum(a2)+taub))
# sample the gamma direction conditional on alpha and beta
gamma[i] <- rgamma(1, shape = alpha0+0.5*n, scale = beta0 + 0.5*sum(a4))
}
b <- burn + 1 # burn-in
parameters <- cbind(alpha,beta,gamma) # the chain, a bivariate sample
para <- parameters[b:N, ]
# compare sample statistics to parameters
colMeans(para)
## alpha beta gamma
## 3.3567376 -0.1949705 773.5494810
cov(para)
## alpha beta gamma
## alpha 1012463.951 -17087.6678 3934.3119
## beta -17087.668 1013277.5316 389.9635
## gamma 3934.312 389.9635 198844.7251
cor(para)
## alpha beta gamma
## alpha 1.000000000 -0.0168705330 0.0087684325
## beta -0.016870533 1.0000000000 0.0008687658
## gamma 0.008768432 0.0008687658 1.0000000000
# 3D Scatterplot
library(scatterplot3d)
## Warning: package 'scatterplot3d' was built under R version 3.2.3
scatterplot3d(para, main="3D Scatterplot")
