library("mvtnorm")
dimensions = 4
num_samples = 1000
burn_in_threshold = 500
# Create a multivariate normal distribution f with some parameters
f_means = runif(dimensions)
f_cov = diag(dimensions)
f = function (x) {
return (dmvnorm(x, f_means, f_cov))
}
# Create a multivariate normal proposal distribution
means = runif(dimensions)
covariance_matrix = diag(dimensions)
# Initialize previous with a random value
previous = rmvnorm(1, mean = means, sigma = covariance_matrix)
# Create an empty data frame to store the accepted values in
df = data.frame(matrix(ncol = dimensions, nrow = 0))
# Iterate till we have data points equal to the number of samples needed
for(i in 1 : num_samples + burn_in_threshold){
# Continue this iteration till we accept a value from the proposal distribution
while(TRUE) {
# Generate a random value from proposal distribution
current = rmvnorm(1, mean = means, sigma = covariance_matrix)
# Calculate alpha using the current and previous values
numerator = f(current) * dmvnorm(previous, means, covariance_matrix)
denominator = f(previous) * dmvnorm(current, means, covariance_matrix)
alpha = numerator / denominator
# Check if the current value is feasible
if (runif(1) < min(alpha, 1)) {
# Append the current value in data frame
df = rbind(df, current)
# Update the means and covariance matrix once we have sufficient accepted values
if(i > burn_in_threshold + dimensions) {
covariance_matrix = cov(df)
means = colMeans(df)
}
# Update the value of previous to current if current is feasible
previous = current
# Move onto the next iteration
break
}
}
}
# Remove burn-in values from the list of accepted values
df = df[-c(1 : burn_in_threshold), ]
# Print the accepted values
print(head(df, 25))
## V1 V2 V3 V4
## 501 0.937182731 0.49291674 0.69840445 0.69461840
## 502 1.174342909 0.02803134 0.31533100 2.05081010
## 503 1.420808724 0.14502141 0.71766710 0.35914315
## 504 1.085062079 -0.07808030 0.05531863 -0.60799106
## 505 -0.338012252 -0.02448536 0.68344554 -0.79689119
## 506 1.494324706 0.83458106 -0.07826402 -0.42163765
## 507 1.073477863 1.22548696 -0.44932056 0.49317400
## 508 0.382647056 0.04853373 -0.17663388 0.94917131
## 509 0.455470819 0.14399613 -1.06026204 0.99344847
## 510 0.985347053 0.26057732 0.65196692 0.21740649
## 511 0.838843270 0.11633613 0.41619646 0.13040587
## 512 -0.672983627 0.01672900 -0.28694120 0.77741949
## 513 0.276158971 1.13218248 -0.33577398 0.58164508
## 514 0.931521858 0.58909366 -0.77476430 1.72889045
## 515 0.005282312 0.81851489 -0.97635552 1.46899673
## 516 0.190307122 0.22009193 0.96367606 -1.03715961
## 517 1.349760368 0.64495642 -1.01122329 0.06045482
## 518 0.412088761 1.29650808 -0.60387789 0.32251642
## 519 0.637034885 -0.45533294 0.99652160 0.81607992
## 520 0.428370736 -0.19250950 0.76709567 -0.78915791
## 521 0.388275335 0.24103710 -0.41677775 1.43941175
## 522 -0.302189517 3.23550660 -2.22052120 0.17375933
## 523 -0.750068523 2.26119969 -0.61469898 1.10502421
## 524 -0.044209594 0.93768727 -0.05111132 -0.26290527
## 525 -0.445806477 0.32199292 0.43202813 0.21099462
# Plot a graph of the accepted values
plot(
df,
col="orange",
main="Point cloud of accepted values"
)
# Estimated means of accepted samples
print(colMeans(df))
## V1 V2 V3 V4
## 0.6788509 0.4416534 -0.2034737 0.5851596
# Original means used for f
print(f_means)
## [1] 0.94177544 0.55743523 0.06759519 0.50682604
# Estimated covariance matrix of accepted samples
print(cov(df))
## V1 V2 V3 V4
## V1 0.81886837 -0.1786098 -0.20848177 0.03079301
## V2 -0.17860982 0.9208532 -0.14277306 0.14798477
## V3 -0.20848177 -0.1427731 0.80460843 0.09152212
## V4 0.03079301 0.1479848 0.09152212 0.86231092
# Original covariance matrix used for f
print(f_cov)
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1