All of the algorithms run on R 3.1.0 on my Mac, but fail on R 3.1.2 on the Linux machines at Uni.

certifiedwaif commented 9 years ago

This is very odd.

Mac log:

R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details.

Natural language support but running in an English locale

R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R.

accuracy.R

source("zero_inflated_model.R") source("generate.R") source("mcmc.R") source("rwmh.R")

calculate_accuracies = function(mult, mcmc_samples, var_result, approximation, print_flag=FALSE, plot_flag=FALSE)

{

TODO: Add support for checking the accuracy over multiple dimensions

cubature$adaptIntegrate

if (plot_flag) pdf(paste0("accuracyplots", approximation, ".pdf"))

return(var_result)

vbeta accuracy

calculate_accuracy = function(mcmc_samples, dist_fn, param1, param2)

{

mcmc_density = density(mcmc_samples)

mcmc_fn = splinefun(mcmc_density$x, mcmc_density$y)

integrand <- function(x)

{

return(abs(mcmc_fn(x) - dist_fn(x, param1, param2)))

}

result = integrate(integrand, min(mcmc_density$x), max(mcmc_density$x),

subdivisions = length(mcmc_density$x))

accuracy = 1 - .5 * result$value

return(accuracy)

}

Compare MCMC distribution with variational approximation for each parameter

vnu[i] ~ Normal, dnorm

sigma2_u ~ IG, dgamma(1/x)

rho ~ Beta, dbeta

vr[i] ~ Bernoulli, dbinom

For each parameter of interest,

* estimate density of MCMC

* compare with q distribution using L_1 norm

Kernel density estimates of MCMC-estimated posteriors

Use L_1 distance to compare against variational approximations of posteriors

accuracy_plot = function(mcmc_samples, dist_fn, param1, param2)

{

mcmc_density = density(mcmc_samples)

plot(mcmc_density)

curve(dist_fn(x, param1, param2),

from=min(mcmc_density$x), to=max(mcmc_density$x),

add=TRUE, lty=2, col="blue")

}

vbeta_accuracy = rep(NA, ncol(mult$mX))

for (i in 1:ncol(mult$mX)) {

vbeta_accuracy[i] = calculate_accuracy(mcmc_samples$vbeta[,i], dnorm,

var_result$vmu[i], sqrt(var_result$mLambda[i,i]))

if (print_flag) cat("vbeta[", i, "]", approximation, "accuracy:", vbeta_accuracy[i], "\n")

param_name = sprintf("vbeta[%d]", i)

if (plot_flag) accuracy_plot(mcmc_samples$vbeta[,i], dnorm,

var_result$vmu[i], sqrt(var_result$mLambda[i,i]))

}

vu accuracy

FIXME: To check for random slopes accuracy, this section will have

to get more complex.

print(dim(mult$mZ))

print(dim(mcmc_samples$u))

vu_accuracy = rep(NA, ncol(mult$mZ))

B = mult$blocksize

b_idx = 1

for (i in 1:ncol(mult$mZ)) {

TODO: The B - (i %% B) expression only works for B=2. Rewrite this.

m_idx = ceiling(i/B)

vu_accuracy[i] = calculate_accuracy(mcmc_samples$vu[,m_idx,b_idx], dnorm,

var_result$vmu[i+mult$p], sqrt(var_result$mLambda[i+mult$p,i+mult$p]))

if (print_flag) cat("vu[", i, "]", approximation, "accuracy:", vu_accuracy[i], "\n")

if (plot_flag) accuracy_plot(mcmc_samples$vu[,m_idx,b_idx], dnorm,

var_result$vmu[i+mult$p], sqrt(var_result$mLambda[i+mult$p,i+mult$p]))

b_idx = b_idx + 1

if (b_idx > B)

b_idx=1

}

sigma2_u accuracy

FIXME - this may be wrong for blocksize != 1?

This is totally wrong for the Inverse Wishart model?

sigma2_u_accuracy = calculate_accuracy(1/mcmc_samples$sigma_u^2, dgamma,

var_result$a_sigma, var_result$b_sigma)

if (print_flag) cat("sigma2_u", approximation, "accuracy:", sigma2_u_accuracy, "\n")

if (plot_flag) accuracy_plot(1/mcmc_samples$sigma_u^2, dgamma,

var_result$a_sigma, var_result$b_sigma)

rho accuracy

rho_accuracy = calculate_accuracy(mcmc_samples$rho, dbeta,

var_result$a_rho, var_result$b_rho)

if (print_flag) cat("rho", approximation, "accuracy: ", rho_accuracy, "\n")

if (plot_flag) accuracy_plot(mcmc_samples$rho, dbeta,

var_result$a_rho, var_result$b_rho)

if (plot_flag) dev.off()

return(list(var_result=var_result,

vbeta_accuracy=vbeta_accuracy,

vu_accuracy=vu_accuracy,

sigma2_u_accuracy=sigma2_u_accuracy,

rho_accuracy=rho_accuracy))

}

test_accuracy = function(mult, mcmc_samples, approximation, plot=FALSE)

{

var_result = zero_infl_var(mult, method=approximation, verbose=TRUE)

return(calculate_accuracies(mult, mcmc_samples, var_result, approximation, plot_flag=plot))

}
Calculate accuracy ----

Approximate the L1 norm between the variational approximation and

the MCMC approximation

calculate_univariate_accuracy <- function(result_mcmc, var_result)

{

density_mcmc_rho = density(result_mcmc$vrho)

integrand <- function(x)

{

fn = splinefun(density_mcmc_rho$x, density_mcmc_rho$y)

return(abs(fn(x) - dbeta(x, var_result$a_rho, var_result$b_rho)))

}

result = integrate(integrand, min(density_mcmc_rho$x), max(density_mcmc_rho$x), subdivisions = length(density_mcmc_rho$x))

rho_accuracy = 1 - .5 * result$value

density_mcmc_lambda = density(result_mcmc$vlambda)

integrand <- function(x)

{

fn = splinefun(density_mcmc_lambda$x, density_mcmc_lambda$y)

return(abs(fn(x) - dgamma(x, var_result$a_lambda, var_result$b_lambda)))

}

result = integrate(integrand, min(density_mcmc_lambda$x), max(density_mcmc_lambda$x), subdivisions = length(density_mcmc_lambda$x))

lambda_accuracy = 1 - .5 * result$value

return(list(rho_accuracy=rho_accuracy, lambda_accuracy=lambda_accuracy))

}

test_accuracies = function()

{

Need to be able to compare the solution paths of each approximation

Generate data

for (i in 1:100) {

set.seed(i)

mult = generate_test_data(20, 100)

Monte Carlo Markov Chains approximation

mcmc_samples = mcmc_approximation(mult, iterations=1e6)

Save the results, because this takes such a long time to run.

}

save(mult, mcmc_samples, file="accuracy_good.RData")

set.seed(1)

mult = generate_test_data(10, 100)

Monte Carlo Markov Chains approximation

mcmc_samples = mcmc_approximation(mult, iterations=1e6, warmup = 1e4)

Save the results, because this takes such a long time to run.

save(mult, mcmc_samples, file="accuracy.RData")

save(mult, mcmc_samples, file="data/accuracy_int.RData")

load(file="data/accuracy_int_2015_02_17.RData")

mult$spline_dim = 0

load(file="accuracy.RData")

Test all other approximations against it

load(file="accuracy.RData")

Test multivariate approximation's accuracy

now = Sys.time()

var1 = test_accuracy(mult, mcmc_samples, "laplacian")

print(Sys.time() - now)

print(image(Matrix(var1$var_result$mLambda)))

print(var1)

now = Sys.time()

var2 = test_accuracy(mult, mcmc_samples, "gva")

print(Sys.time() - now)

print(image(Matrix(var2$var_result$mLambda)))

print(var2)

now = Sys.time()

var3 = test_accuracy(mult, mcmc_samples, "gva2", plot=TRUE)

print(Sys.time() - now)

print(image(Matrix(var3$mLambda)))

print(var3)

Rprof()

now = Sys.time()

var4 = test_accuracy(mult, mcmc_samples, "gva2new")

print(Sys.time() - now)

print(image(Matrix(var4$var_result$mLambda)))

print(var4)

Rprof(NULL)

print(summaryRprof())

print(image(Matrix(var3_new$mLambda)))

now = Sys.time()

var5 = test_accuracy(mult, mcmc_samples, "gva_nr")

print(Sys.time() - now)

print(image(Matrix(var4$result_var$mLambda)))

print(var5)

save(var1, var2, var3, var4, var5, file="accuracy_results_int.RData")

for (i in 1:100) {

set.seed(i)

mult = generate_test_data(20, 100)

mcmc_samples = mcmc_approximation(mult, iterations=1e4)

var1 = test_accuracy(mult, mcmc_samples, "laplacian")

var2 = test_accuracy(mult, mcmc_samples, "gva")

var3 = test_accuracy(mult, mcmc_samples, "gva2")

var4 = test_accuracy(mult, mcmc_samples, "gva_nr")

}

}
test_accuracies()

test_accuracies_slope = function()

{

Monte Carlo Markov Chains approximation

set.seed(1)

mult = generate_slope_test_data()

mcmc_samples = mcmc_approximation(mult, iterations=1e6, warmup=1e4, mc.cores = 1)

save(mult, mcmc_samples, file="data/accuracy_slope_2015_03_03.RData")

load(file="data/accuracy_slope_2015_03_03.RData")

now = Sys.time()

var1 = test_accuracy(mult, mcmc_samples, "laplacian", plot=TRUE)

print(Sys.time() - now)

print(var1)

now = Sys.time()

var2 = test_accuracy(mult, mcmc_samples, "gva", plot=TRUE)

print(Sys.time() - now)

print(image(Matrix(var2$var_result$mLambda)))

print(var2)

now = Sys.time()

var3 = test_accuracy(mult, mcmc_samples, "gva2", plot=TRUE)

print(image(Matrix(var3$var_result$mLambda)))

print(Sys.time() - now)

print(var3)

print("gva2new")

now = Sys.time()

Rprof(line.profiling=TRUE)

var4 = test_accuracy(mult, mcmc_samples, "gva2new", plot=TRUE)

Rprof(NULL)

print(summaryRprof(lines = "both"))

print(image(Matrix(var4$var_result$mLambda)))

var4 = test_accuracy(mult, mcmc_samples, "gva2new")

print(Sys.time() - now)

print(var4)

now = Sys.time()

GVA NR is unstable, and sometimes fails with an error

var5 = test_accuracy(mult, mcmc_samples, "gva_nr")

var5 = test_accuracy(mult, mcmc_samples, "gva_nr", plot=TRUE)

print(Sys.time() - now)

print(image(Matrix(var5$var_result$mLambda)))

print(var5)

save(var1, var2, var3, var4, var5, file="accuracy_results_slope.RData")

print(image(Matrix(var3_new$mLambda)))

now = Sys.time()

var4 = test_accuracy(mult, mcmc_samples, "gva_nr")

Sys.time() - now

print(image(Matrix(var4$mLambda)))

mult = generate_spline_test_data()

mcmc_samples = mcmc_approximation(mult, iterations=1e4)

save(mult, mcmc_samples, file="accuracy_spline.RData")

now = Sys.time()

var1 = test_accuracy(mult, mcmc_samples, "laplacian")

Sys.time() - now

print(image(Matrix(var1$var_result$mLambda)))

var1

now = Sys.time()

var2 = test_accuracy(mult, mcmc_samples, "gva")

Sys.time() - now

print(image(Matrix(var2$var_result$mLambda)))

var2

now = Sys.time()

var3 = test_accuracy(mult, mcmc_samples, "gva2")

Sys.time() - now

print(image(Matrix(var3$var_result$mLambda)))

var3

now = Sys.time()

var4 = test_accuracy(mult, mcmc_samples, "gva2new")

Sys.time() - now

print(image(Matrix(var4$var_result$mLambda)))

var4

Rprof()

now = Sys.time()

var5 = test_accuracy(mult, mcmc_samples, "gva_nr")

print(Sys.time() - now)

print(image(Matrix(var5$var_result$mLambda)))

var5

} test_accuracies_slope() N 100 p 2 m 5 blocksize 2 spline_dim 0 calculate_lower_bound: T1 -4.10617e+17 T2 -105.5615 T3 -38.79287 a_sigma b_sigma calculate_lower_bound: T1 -4254.533 T2 -48.13332 T3 -47.93269 Iteration 2 : lower bound -4350.599 difference 4.10617e+17 parameters vmu 1.970248 2.387512 1.101093 -2.379948 1.825163 -2.837605 2.43698 -2.691774 -0.1999721 -1.972204 diag(mLambda) 0.008005845 0.00468602 0.01032378 0.007912696 0.009163322 0.006823373 0.008631102 0.005365371 0.01762683 0.01577659 a_rho 84.0028 b_rho 17.9972 a_sigma b_sigma calculate_lower_bound: T1 -4248.745 T2 -48.09777 T3 -47.92874 Iteration 3 : lower bound -4344.772 difference 5.827025 parameters vmu 2.093207 2.303359 0.9793739 -2.292995 1.701905 -2.752487 2.39622 -2.570463 -0.3002183 -1.871611 diag(mLambda) 0.007165469 0.004302308 0.009434513 0.007534457 0.008300133 0.006445907 0.007727148 0.004922677 0.01642191 0.01512188 a_rho 84.00333 b_rho 17.99667 a_sigma b_sigma calculate_lower_bound: T1 -4248.74 T2 -48.10542 T3 -47.92929 Iteration 4 : lower bound -4344.774 difference -0.002537746 parameters vmu 2.09571 2.301589 0.9766196 -2.291273 1.699333 -2.750677 2.393756 -2.568625 -0.3045866 -1.870705 diag(mLambda) 0.007154014 0.00429596 0.009425739 0.007527698 0.008289927 0.006438819 0.007716174 0.004915965 0.01642994 0.01513047 a_rho 84.00326 b_rho 17.99674 a_sigma b_sigma [1] 100 8 NULL Time difference of 2.466712 secs $var_result $var_result$vmu ,1 2.0957098 x 2.3015888 groups2 0.9766196 x:groups2 -2.2912729 groups3 1.6993328 x:groups3 -2.7506774 groups4 2.3937563 x:groups4 -2.5686248 groups5 -0.3045866 x:groups5 -1.8707046

$var_result$mLambda (Intercept) x groups2 x:groups2 groups3 (Intercept) 0.007154014 -0.005154323 -0.007136626 0.005163922 -0.007144532 x -0.005154323 0.004295960 0.005143068 -0.004301498 0.005148240 groups2 -0.007136626 0.005143068 0.009425739 -0.005162749 0.007127168 x:groups2 0.005163922 -0.004301498 -0.005162749 0.007527698 -0.005157824 groups3 -0.007144532 0.005148240 0.007127168 -0.005157824 0.008289927 x:groups3 0.005162642 -0.004300897 -0.005151365 0.004306450 -0.004903491 groups4 -0.007149408 0.005151363 0.007132031 -0.005160954 0.007139932 x:groups4 0.005156909 -0.004297504 -0.005145647 0.004303047 -0.005150822 groups5 -0.007097451 0.005116721 0.007080207 -0.005126237 0.007088048 x:groups5 0.005150843 -0.004291441 -0.005139592 0.004296980 -0.005144762 x:groups3 groups4 x:groups4 groups5 x:groups5 (Intercept) 0.005162642 -0.007149408 0.005156909 -0.007097451 0.005150843 x -0.004300897 0.005151363 -0.004297504 0.005116721 -0.004291441 groups2 -0.005151365 0.007132031 -0.005145647 0.007080207 -0.005139592 x:groups2 0.004306450 -0.005160954 0.004303047 -0.005126237 0.004296980 groups3 -0.004903491 0.007139932 -0.005150822 0.007088048 -0.005144762 x:groups3 0.006438819 -0.005159675 0.004302446 -0.005124969 0.004296379 groups4 -0.005159675 0.007716174 -0.005056923 0.007092883 -0.005147884 x:groups4 0.004302446 -0.005056923 0.004915965 -0.005119285 0.004292986 groups5 -0.005124969 0.007092883 -0.005119285 0.016429939 -0.009887394 x:groups5 0.004296379 -0.005147884 0.004292986 -0.009887394 0.015130470

$var_result$a_rho [1] 84.00326

$var_result$b_rho [1] 17.99674

$var_result$a_sigma numeric(0)

$var_result$b_sigma numeric(0)

$var_result$vlower_bound [1] -4.106170e+17 -4.350599e+03 -4.344772e+03 -4.344774e+03

$vbeta_accuracy [1] 0.9070440 0.8987978

$vu_accuracy [1] 0.9472439 0.8896154 0.9140297 0.9021092 0.9046075 0.9117498 0.8924308 [8] 0.8997999

$rho_accuracy [1] 0.9036479

N 100 p 2 m 5 blocksize 2 spline_dim 0 calculate_lower_bound: T1 -4290.5 T2 -49.87337 T3 -47.95629 a_sigma b_sigma calculate_lower_bound: T1 -4248.556 T2 -48.19411 T3 -47.92917 Iteration 2 : lower bound -4344.679 difference 43.65032 parameters vmu 2.09413 2.302459 0.9748742 -2.291327 1.699557 -2.752339 2.394875 -2.569562 -0.3210693 -1.869889 diag(mLambda) 0.008482581 0.004920235 0.01120775 0.00882215 0.009749764 0.007525693 0.009149858 0.005617782 0.02105552 0.02150821 a_rho 84.00327 b_rho 17.99673 a_sigma b_sigma calculate_lower_bound: T1 -4248.571 T2 -48.15914 T3 -47.92901 Iteration 3 : lower bound -4344.66 difference 0.01938413 parameters vmu 2.093147 2.302883 0.9762678 -2.291192 1.700575 -2.752579 2.395646 -2.570119 -0.315083 -1.866814 diag(mLambda) 0.008453674 0.004908523 0.01116801 0.008813112 0.009715358 0.007520294 0.009118412 0.005608128 0.02092709 0.0214335 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15807 T3 -47.92899 Iteration 4 : lower bound -4344.659 difference 0.000556743 parameters vmu 2.093039 2.302952 0.976395 -2.291246 1.700687 -2.752644 2.395749 -2.570193 -0.3148 -1.866767 diag(mLambda) 0.008453144 0.00490834 0.01116712 0.008812989 0.009714664 0.007520192 0.009117801 0.005607992 0.020923 0.02143089 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 5 : lower bound -4344.659 difference 2.094874e-05 parameters vmu 2.093033 2.302956 0.9764014 -2.291249 1.700693 -2.752648 2.395754 -2.570197 -0.314788 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812984 0.009714661 0.007520192 0.0091178 0.005607994 0.0209229 0.02143087 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 6 : lower bound -4344.659 difference 1.130324e-06 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395755 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453147 0.004908342 0.01116711 0.008812985 0.009714661 0.007520193 0.009117801 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 7 : lower bound -4344.659 difference -8.684765e-08 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395754 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453147 0.004908343 0.01116711 0.008812991 0.009714661 0.007520196 0.0091178 0.005607996 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma [1] 100 8 NULL Time difference of 1.995439 secs $var_result $var_result$vmu [1] 2.0930326 2.3029559 0.9764019 -2.2912491 1.7006935 -2.7526478 [7] 2.3957545 -2.5701973 -0.3147871 -1.8667661

$var_result$mLambda [,1] [,2] [,3] [,4] [,5] [1,] 0.008453147 -0.006044301 -0.008430261 0.006054040 -0.008441129 [2,] -0.006044301 0.004908343 0.006029263 -0.004913887 0.006036490 [3,] -0.008430261 0.006029263 0.011167114 -0.006484650 0.008418278 [4,] 0.006054040 -0.004913887 -0.006484650 0.008812991 -0.006046214 [5,] -0.008441129 0.006036490 0.008418278 -0.006046214 0.009714661 [6,] 0.006055018 -0.004914797 -0.006039948 0.004920357 -0.005868733 [7,] -0.008446603 0.006040060 0.008423736 -0.006049790 0.008434595 [8,] 0.006048103 -0.004910680 -0.006033053 0.004916230 -0.006040286 [9,] -0.008360108 0.005982407 0.008337485 -0.005992026 0.008348228 [10,] 0.006054226 -0.004911368 -0.006039150 0.004916932 -0.006046394 [,6] [,7] [,8] [,9] [,10] [1,] 0.006055018 -0.008446603 0.006048103 -0.008360108 0.006054226 [2,] -0.004914797 0.006040060 -0.004910680 0.005982407 -0.004911368 [3,] -0.006039948 0.008423736 -0.006033053 0.008337485 -0.006039150 [4,] 0.004920357 -0.006049790 0.004916230 -0.005992026 0.004916932 [5,] -0.005868733 0.008434595 -0.006040286 0.008348228 -0.006046394 [6,] 0.007520196 -0.006050767 0.004917140 -0.005992995 0.004917842 [7,] -0.006050767 0.009117800 -0.005890535 0.008353639 -0.006049972 [8,] 0.004917140 -0.005890535 0.005607996 -0.005986162 0.004913711 [9,] -0.005992995 0.008353639 -0.005986162 0.020922893 -0.011423462 [10,] 0.004917842 -0.006049972 0.004913711 -0.011423462 0.021430839

$var_result$a_rho [1] 84.0033

$var_result$b_rho [1] 17.9967

$var_result$a_sigma numeric(0)

$var_result$b_sigma numeric(0)

$var_result$vlower_bound [1] -4388.329 -4344.679 -4344.660 -4344.659 -4344.659 -4344.659 -4344.659

$vbeta_accuracy [1] 0.9218658 0.9196421

$vu_accuracy [1] 0.9555139 0.9244812 0.9300057 0.9253950 0.9257944 0.9294784 0.9325449 [8] 0.9585504

$rho_accuracy [1] 0.9036515

N 100 p 2 m 5 blocksize 2 spline_dim 0 calculate_lower_bound: T1 -4.10617e+17 T2 -105.5615 T3 -38.79287 a_sigma b_sigma calculate_lower_bound: T1 -4254.581 T2 -48.17991 T3 -47.93239 Iteration 2 : lower bound -4350.693 difference 4.10617e+17 parameters vmu 1.967162 2.38913 1.101495 -2.380752 1.826973 -2.839919 2.439368 -2.693567 -0.2090716 -1.967071 diag(mLambda) 0.009626389 0.005505333 0.0123472 0.009262223 0.01084605 0.008085332 0.01031284 0.006211642 0.02329809 0.02103267 a_rho 84.00284 b_rho 17.99716 a_sigma b_sigma calculate_lower_bound: T1 -4248.666 T2 -48.1396 T3 -47.92847 Iteration 3 : lower bound -4344.734 difference 5.959453 parameters vmu 2.089955 2.305214 0.9802474 -2.2943 1.70389 -2.755052 2.398955 -2.572443 -0.3080842 -1.8672 diag(mLambda) 0.009357424 0.005479383 0.01210666 0.009380454 0.01067038 0.007993289 0.009852098 0.006041833 0.02253293 0.0211349 a_rho 84.00337 b_rho 17.99663 a_sigma b_sigma calculate_lower_bound: T1 -4248.657 T2 -48.15205 T3 -47.92845 Iteration 4 : lower bound -4344.738 difference -0.004210441 parameters vmu 2.089665 2.305241 0.9794486 -2.294553 1.704081 -2.754486 2.399184 -2.57229 -0.3092569 -1.868841 diag(mLambda) 0.009506722 0.005478045 0.01233618 0.009400417 0.01077984 0.008051188 0.009967946 0.006027864 0.02250331 0.02121321 a_rho 84.00337 b_rho 17.99663 a_sigma b_sigma [1] 100 8 NULL Time difference of 2.084503 secs $var_result $var_result$vmu [1] 2.0896655 2.3052413 0.9794486 -2.2945533 1.7040812 -2.7544864 [7] 2.3991839 -2.5722897 -0.3092569 -1.8688408

$var_result$mLambda (Intercept) x groups2 x:groups2 groups3 (Intercept) 0.009506722 -0.006810024 -0.009693908 0.007077483 -0.009543433 x -0.006810024 0.005478045 0.006940870 -0.005700327 0.006831581 groups2 -0.009693908 0.006940870 0.012336184 -0.007224387 0.009722892 x:groups2 0.007077483 -0.005700327 -0.007224387 0.009400417 -0.007097980 groups3 -0.009543433 0.006831581 0.009722892 -0.007097980 0.010779844 x:groups3 0.006991532 -0.005632224 -0.007123440 0.005849564 -0.006738496 groups4 -0.009445196 0.006760559 0.009622777 -0.007024161 0.009474876 x:groups4 0.006747138 -0.005439216 -0.006875162 0.005649940 -0.006767274 groups5 -0.010435842 0.007494819 0.010619410 -0.007773163 0.010468587 x:groups5 0.008099759 -0.006516700 -0.008236823 0.006750237 -0.008124073 x:groups3 groups4 x:groups4 groups5 x:groups5 (Intercept) 0.006991532 -0.009445196 0.006747138 -0.010435842 0.008099759 x -0.005632224 0.006760559 -0.005439216 0.007494819 -0.006516700 groups2 -0.007123440 0.009622777 -0.006875162 0.010619410 -0.008236823 x:groups2 0.005849564 -0.007024161 0.005649940 -0.007773163 0.006750237 groups3 -0.006738496 0.009474876 -0.006767274 0.010468587 -0.008124073 x:groups3 0.008051188 -0.006940005 0.005583984 -0.007687403 0.006681336 groups4 -0.006940005 0.009967946 -0.006592944 0.010364304 -0.008048975 x:groups4 0.005583984 -0.006592944 0.006027864 -0.007419028 0.006467086 groups5 -0.007687403 0.010364304 -0.007419028 0.022503309 -0.015027962 x:groups5 0.006681336 -0.008048975 0.006467086 -0.015027962 0.021213209

$var_result$a_rho [1] 84.00337

$var_result$b_rho [1] 17.99663

$var_result$a_sigma numeric(0)

$var_result$b_sigma numeric(0)

$var_result$vlower_bound [1] -4.106170e+17 -4.350693e+03 -4.344734e+03 -4.344738e+03

$vbeta_accuracy [1] 0.9223178 0.9321255

$vu_accuracy [1] 0.9530949 0.9346686 0.9326032 0.9323161 0.9308158 0.9364944 0.9189237 [8] 0.9602546

$rho_accuracy [1] 0.9036579

[1] "gva2new" N 100 p 2 m 5 blocksize 2 spline_dim 0 calculate_lower_bound: T1 -4292.737 T2 -49.86207 T3 -47.95629 a_sigma b_sigma calculate_lower_bound: T1 -4250.696 T2 -48.16968 T3 -47.9283 Iteration 2 : lower bound -4346.794 difference 43.76168 parameters vmu 2.097054 2.300623 0.9717529 -2.289799 1.696478 -2.750359 2.391713 -2.567674 -0.3238487 -1.869233 diag(mLambda) 0.0002772039 0.0002585945 0.002714801 0.003825567 0.001274922 0.002529919 0.0006452087 0.0006648273 0.01206992 0.01604501 a_rho 84.00339 b_rho 17.99661 a_sigma b_sigma calculate_lower_bound: T1 -4250.713 T2 -48.13132 T3 -47.92848 Iteration 3 : lower bound -4346.773 difference 0.02088818 parameters vmu 2.097566 2.299943 0.9716072 -2.288637 1.695939 -2.749491 2.390994 -2.567088 -0.3199915 -1.86376 diag(mLambda) 0.0002773868 0.0002585225 0.002711821 0.003823748 0.001274627 0.002565101 0.0006452489 0.0006651087 0.01199677 0.01599983 a_rho 84.00337 b_rho 17.99663 a_sigma b_sigma calculate_lower_bound: T1 -4250.709 T2 -48.13713 T3 -47.92777 Iteration 4 : lower bound -4346.774 difference -0.001010009 parameters vmu 2.094479 2.302267 0.9748108 -2.291137 1.6991 -2.751906 2.394069 -2.569463 -0.3155559 -1.868097 diag(mLambda) 0.0002774194 0.0002585368 0.00271153 0.003823654 0.001274604 0.002568743 0.0006452575 0.0006651432 0.01199424 0.01599953 a_rho 84.00347 b_rho 17.99653 a_sigma b_sigma [1] 100 8 NULL Time difference of 2.075891 secs $var_result $var_result$vmu [1] 2.0944786 2.3022674 0.9748108 -2.2911367 1.6990997 -2.7519064 [7] 2.3940692 -2.5694626 -0.3155559 -1.8680974

$var_result$mLambda 1 2 1 2 3 1 2.774194e-04 1.865033e-06 -2.463851e-05 -5.330677e-06 -1.859036e-05 2 1.865033e-06 2.585368e-04 -3.162288e-06 -2.158638e-05 -7.083789e-07 1 -2.463851e-05 -3.162288e-06 2.711530e-03 7.232348e-07 1.657832e-06 2 -5.330677e-06 -2.158638e-05 7.232348e-07 3.823654e-03 4.058501e-07 3 -1.859036e-05 -7.083789e-07 1.657832e-06 4.058501e-07 1.274604e-03 4 -2.425131e-06 -1.951697e-05 4.414230e-07 1.672172e-06 2.065185e-07 5 -1.106122e-05 7.271966e-06 8.972287e-07 -3.998449e-07 7.246534e-07 6 7.312536e-06 -3.111613e-08 -6.485195e-07 -1.338202e-07 -4.898446e-07 7 -4.127052e-05 9.158416e-07 3.651536e-06 6.935497e-07 2.762916e-06 8 -2.072022e-05 -4.610888e-05 2.373080e-06 4.230161e-06 1.492236e-06 4 5 6 7 8 1 -2.425131e-06 -1.106122e-05 7.312536e-06 -4.127052e-05 -2.072022e-05 2 -1.951697e-05 7.271966e-06 -3.111613e-08 9.158416e-07 -4.610888e-05 1 4.414230e-07 8.972287e-07 -6.485195e-07 3.651536e-06 2.373080e-06 2 1.672172e-06 -3.998449e-07 -1.338202e-07 6.935497e-07 4.230161e-06 3 2.065185e-07 7.246534e-07 -4.898446e-07 2.762916e-06 1.492236e-06 4 2.568743e-03 -4.574443e-07 -5.786903e-08 2.707655e-07 3.648649e-06 5 -4.574443e-07 6.452575e-04 2.125136e-05 1.679440e-06 -4.801365e-07 6 -5.786903e-08 2.125136e-05 6.651432e-04 -1.088226e-06 -5.318928e-07 7 2.707655e-07 1.679440e-06 -1.088226e-06 1.199424e-02 -4.195596e-03 8 3.648649e-06 -4.801365e-07 -5.318928e-07 -4.195596e-03 1.599953e-02

$var_result$a_rho [1] 84.00347

$var_result$b_rho [1] 17.99653

$var_result$a_sigma numeric(0)

$var_result$b_sigma numeric(0)

$var_result$vlower_bound [1] -4390.555 -4346.794 -4346.773 -4346.774

$vbeta_accuracy [1] 0.3031139 0.3435177

$vu_accuracy [1] 0.6712743 0.7325994 0.5152712 0.6898212 0.4117733 0.4804914 0.8437027 [8] 0.9131477

$rho_accuracy [1] 0.9036661

N 100 p 2 m 5 blocksize 2 spline_dim 0 calculate_lower_bound: T1 -4290.483 T2 -49.87728 T3 -47.95629 a_sigma b_sigma calculate_lower_bound: T1 -4248.556 T2 -48.19414 T3 -47.92916 Iteration 2 : lower bound -4344.679 difference 43.63771 parameters vmu 2.094092 2.30249 0.9749076 -2.291343 1.699598 -2.752365 2.394914 -2.569594 -0.3210242 -1.869919 diag(mLambda) 0.008481814 0.004919978 0.01120775 0.008823185 0.009748932 0.007530014 0.009149182 0.005618411 0.02106534 0.02151484 a_rho 84.00328 b_rho 17.99672 a_sigma b_sigma calculate_lower_bound: T1 -4248.571 T2 -48.15913 T3 -47.92901 Iteration 3 : lower bound -4344.66 difference 0.01938018 parameters vmu 2.093146 2.302883 0.9762684 -2.291192 1.700575 -2.752579 2.395646 -2.570119 -0.3150817 -1.866813 diag(mLambda) 0.008453264 0.004908316 0.01116759 0.008812906 0.009714952 0.007520098 0.009118001 0.005607927 0.02092649 0.02143308 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15807 T3 -47.92899 Iteration 4 : lower bound -4344.659 difference 0.00055306 parameters vmu 2.093039 2.302952 0.9763949 -2.291246 1.700687 -2.752644 2.395749 -2.570193 -0.3148001 -1.866767 diag(mLambda) 0.008453131 0.004908331 0.01116711 0.008812973 0.009714652 0.007520181 0.009117788 0.005607983 0.02092299 0.02143091 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 5 : lower bound -4344.659 difference 2.0825e-05 parameters vmu 2.093033 2.302956 0.9764015 -2.291249 1.700693 -2.752648 2.395754 -2.570197 -0.3147877 -1.866766 diag(mLambda) 0.008453145 0.004908341 0.01116711 0.008812985 0.00971466 0.007520193 0.009117799 0.005607994 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 6 : lower bound -4344.659 difference 8.548432e-07 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395754 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812986 0.009714661 0.007520194 0.0091178 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 7 : lower bound -4344.659 difference 3.64098e-08 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395755 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812986 0.009714661 0.007520194 0.0091178 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 8 : lower bound -4344.659 difference 1.578883e-09 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395755 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812986 0.009714661 0.007520194 0.0091178 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 9 : lower bound -4344.659 difference 6.548362e-11 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395755 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812986 0.009714661 0.007520194 0.0091178 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 10 : lower bound -4344.659 difference 2.728484e-12 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395755 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812986 0.009714661 0.007520194 0.0091178 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 11 : lower bound -4344.659 difference 1.818989e-12 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395755 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812986 0.009714661 0.007520194 0.0091178 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 12 : lower bound -4344.659 difference 0 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395755 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812986 0.009714661 0.007520194 0.0091178 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma [1] 100 8 NULL Time difference of 1.688941 secs $var_result $var_result$vmu [,1] [1,] 2.0930326 [2,] 2.3029559 [3,] 0.9764019 [4,] -2.2912491 [5,] 1.7006935 [6,] -2.7526478 [7,] 2.3957545 [8,] -2.5701973 [9,] -0.3147871 [10,] -1.8667662

$var_result$mLambda [,1] [,2] [,3] [,4] [,5] [1,] 0.008453146 -0.006044300 -0.008430260 0.006054035 -0.008441129 [2,] -0.006044300 0.004908342 0.006029262 -0.004913884 0.006036490 [3,] -0.008430260 0.006029262 0.011167113 -0.006484645 0.008418277 [4,] 0.006054035 -0.004913884 -0.006484645 0.008812986 -0.006046209 [5,] -0.008441129 0.006036490 0.008418277 -0.006046209 0.009714661 [6,] 0.006055016 -0.004914795 -0.006039945 0.004920353 -0.005868731 [7,] -0.008446602 0.006040059 0.008423735 -0.006049785 0.008434595 [8,] 0.006048102 -0.004910678 -0.006033051 0.004916226 -0.006040285 [9,] -0.008360106 0.005982404 0.008337483 -0.005992021 0.008348227 [10,] 0.006054223 -0.004911365 -0.006039149 0.004916929 -0.006046393 [,6] [,7] [,8] [,9] [,10] [1,] 0.006055016 -0.008446602 0.006048102 -0.008360106 0.006054223 [2,] -0.004914795 0.006040059 -0.004910678 0.005982404 -0.004911365 [3,] -0.006039945 0.008423735 -0.006033051 0.008337483 -0.006039149 [4,] 0.004920353 -0.006049785 0.004916226 -0.005992021 0.004916929 [5,] -0.005868731 0.008434595 -0.006040285 0.008348227 -0.006046393 [6,] 0.007520194 -0.006050765 0.004917138 -0.005992992 0.004917840 [7,] -0.006050765 0.009117800 -0.005890534 0.008353638 -0.006049971 [8,] 0.004917138 -0.005890534 0.005607995 -0.005986160 0.004913709 [9,] -0.005992992 0.008353638 -0.005986160 0.020922889 -0.011423464 [10,] 0.004917840 -0.006049971 0.004913709 -0.011423464 0.021430839

$var_result$a_rho [1] 84.0033

$var_result$b_rho [1] 17.9967

$var_result$a_sigma numeric(0)

$var_result$b_sigma numeric(0)

$var_result$vlower_bound [1] -4388.317 -4344.679 -4344.660 -4344.659 -4344.659 -4344.659 -4344.659 [8] -4344.659 -4344.659 -4344.659 -4344.659 -4344.659

$vbeta_accuracy [1] 0.9218658 0.9196421

$vu_accuracy [1] 0.9555139 0.9244811 0.9300057 0.9253950 0.9257944 0.9294784 0.9325448 [8] 0.9585505

$rho_accuracy [1] 0.9036515

Linux log:

R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R.

accuracy.R

source("zero_inflated_model.R") source("generate.R") source("mcmc.R") source("rwmh.R")

calculate_accuracy = function(mcmc_samples, dist_fn, param1, param2)

{

mcmc_density = density(mcmc_samples)

mcmc_fn = splinefun(mcmc_density$x, mcmc_density$y)

integrand <- function(x)

{

return(abs(mcmc_fn(x) - dist_fn(x, param1, param2)))

}

result = integrate(integrand, min(mcmc_density$x), max(mcmc_density$x),

subdivisions = length(mcmc_density$x))

accuracy = 1 - .5 * result$value

return(accuracy)

}

accuracy_plot = function(mcmc_samples, dist_fn, param1, param2)

{

mcmc_density = density(mcmc_samples)

plot(mcmc_density)

curve(dist_fn(x, param1, param2),

from=min(mcmc_density$x), to=max(mcmc_density$x),

add=TRUE, lty=2, col="blue")

}

calculate_accuracies = function(mult, mcmc_samples, var_result, approximation, print_flag=FALSE, plot_flag=FALSE)

{

TODO: Add support for checking the accuracy over multiple dimensions

cubature$adaptIntegrate

if (plot_flag) pdf(paste0("accuracyplots", approximation, ".pdf"))

return(var_result)

vbeta accuracy

Compare MCMC distribution with variational approximation for each parameter

vnu[i] ~ Normal, dnorm

sigma2_u ~ IG, dgamma(1/x)

rho ~ Beta, dbeta

vr[i] ~ Bernoulli, dbinom

For each parameter of interest,

* estimate density of MCMC

* compare with q distribution using L_1 norm

Kernel density estimates of MCMC-estimated posteriors

Use L_1 distance to compare against variational approximations of posteriors

vbeta_accuracy = rep(NA, ncol(mult$mX))

for (i in 1:ncol(mult$mX)) {

vbeta_accuracy[i] = calculate_accuracy(mcmc_samples$vbeta[,i], dnorm,

var_result$vmu[i], sqrt(var_result$mLambda[i,i]))

if (print_flag) cat("vbeta[", i, "]", approximation, "accuracy:", vbeta_accuracy[i], "\n")

param_name = sprintf("vbeta[%d]", i)

if (plot_flag) accuracy_plot(mcmc_samples$vbeta[,i], dnorm,

var_result$vmu[i], sqrt(var_result$mLambda[i,i]))

}

vu accuracy

FIXME: To check for random slopes accuracy, this section will have

to get more complex.

print(dim(mult$mZ))

print(dim(mcmc_samples$u))

vu_accuracy = rep(NA, ncol(mult$mZ))

B = mult$blocksize

b_idx = 1

for (i in 1:ncol(mult$mZ)) {

m_idx = ceiling(i/B)

vu_accuracy[i] = calculate_accuracy(mcmc_samples$vu[,m_idx,b_idx], dnorm,

var_result$vmu[i+mult$p], sqrt(var_result$mLambda[i+mult$p,i+mult$p]))

if (print_flag) cat("vu[", i, "]", approximation, "accuracy:", vu_accuracy[i], "\n")

if (plot_flag) accuracy_plot(mcmc_samples$vu[,m_idx,b_idx], dnorm,

var_result$vmu[i+mult$p], sqrt(var_result$mLambda[i+mult$p,i+mult$p]))

b_idx = b_idx + 1

if (b_idx > B)

b_idx=1

}

sigma2_u accuracy

FIXME - this may be wrong for blocksize != 1?

This is totally wrong for the Inverse Wishart model?

sigma2_u_accuracy = calculate_accuracy(1/mcmc_samples$sigma_u^2, dgamma,

var_result$a_sigma, var_result$b_sigma)

if (print_flag) cat("sigma2_u", approximation, "accuracy:", sigma2_u_accuracy, "\n")

if (plot_flag) accuracy_plot(1/mcmc_samples$sigma_u^2, dgamma,

var_result$a_sigma, var_result$b_sigma)

rho accuracy

rho_accuracy = calculate_accuracy(mcmc_samples$rho, dbeta,

var_result$a_rho, var_result$b_rho)

if (print_flag) cat("rho", approximation, "accuracy: ", rho_accuracy, "\n")

if (plot_flag) accuracy_plot(mcmc_samples$rho, dbeta,

var_result$a_rho, var_result$b_rho)

if (plot_flag) dev.off()

return(list(var_result=var_result,

vbeta_accuracy=vbeta_accuracy,

vu_accuracy=vu_accuracy,

sigma2_u_accuracy=sigma2_u_accuracy,

rho_accuracy=rho_accuracy))

}

test_accuracy = function(mult, mcmc_samples, approximation, plot=FALSE)

{

var_result = zero_infl_var(mult, method=approximation, verbose=TRUE)

return(calculate_accuracies(mult, mcmc_samples, var_result, approximation, plot_flag=plot))

}
Calculate accuracy ----

Approximate the L1 norm between the variational approximation and

the MCMC approximation

calculate_univariate_accuracy <- function(result_mcmc, var_result)

{

density_mcmc_rho = density(result_mcmc$vrho)

integrand <- function(x)

{

fn = splinefun(density_mcmc_rho$x, density_mcmc_rho$y)

return(abs(fn(x) - dbeta(x, var_result$a_rho, var_result$b_rho)))

}

result = integrate(integrand, min(density_mcmc_rho$x), max(density_mcmc_rho$x), subdivisions = length(density_mcmc_rho$x))

rho_accuracy = 1 - .5 * result$value

density_mcmc_lambda = density(result_mcmc$vlambda)

integrand <- function(x)

{

fn = splinefun(density_mcmc_lambda$x, density_mcmc_lambda$y)

return(abs(fn(x) - dgamma(x, var_result$a_lambda, var_result$b_lambda)))

}

result = integrate(integrand, min(density_mcmc_lambda$x), max(density_mcmc_lambda$x), subdivisions = length(density_mcmc_lambda$x))

lambda_accuracy = 1 - .5 * result$value

return(list(rho_accuracy=rho_accuracy, lambda_accuracy=lambda_accuracy))

}

test_accuracies = function()

{

Need to be able to compare the solution paths of each approximation

Generate data

for (i in 1:100) {

set.seed(i)

mult = generate_test_data(20, 100)

Monte Carlo Markov Chains approximation

mcmc_samples = mcmc_approximation(mult, iterations=1e6)

Save the results, because this takes such a long time to run.

}

save(mult, mcmc_samples, file="accuracy_good.RData")

set.seed(1)

mult = generate_test_data(10, 100)

Monte Carlo Markov Chains approximation

mcmc_samples = mcmc_approximation(mult, iterations=1e6, warmup = 1e4)

Save the results, because this takes such a long time to run.

save(mult, mcmc_samples, file="accuracy.RData")

save(mult, mcmc_samples, file="data/accuracy_int.RData")

load(file="data/accuracy_int_2015_02_17.RData")

mult$spline_dim = 0

load(file="accuracy.RData")

Test all other approximations against it

load(file="accuracy.RData")

Test multivariate approximation's accuracy

now = Sys.time()

var1 = test_accuracy(mult, mcmc_samples, "laplacian")

print(Sys.time() - now)

print(image(Matrix(var1$var_result$mLambda)))

print(var1)

now = Sys.time()

var2 = test_accuracy(mult, mcmc_samples, "gva")

print(Sys.time() - now)

print(image(Matrix(var2$var_result$mLambda)))

print(var2)

now = Sys.time()

var3 = test_accuracy(mult, mcmc_samples, "gva2", plot=TRUE)

print(Sys.time() - now)

print(image(Matrix(var3$mLambda)))

print(var3)

Rprof()

now = Sys.time()

var4 = test_accuracy(mult, mcmc_samples, "gva2new")

print(Sys.time() - now)

print(image(Matrix(var4$var_result$mLambda)))

print(var4)

Rprof(NULL)

print(summaryRprof())

print(image(Matrix(var3_new$mLambda)))

now = Sys.time()

var5 = test_accuracy(mult, mcmc_samples, "gva_nr")

print(Sys.time() - now)

print(image(Matrix(var4$result_var$mLambda)))

print(var5)

save(var1, var2, var3, var4, var5, file="accuracy_results_int.RData")

for (i in 1:100) {

set.seed(i)

mult = generate_test_data(20, 100)

mcmc_samples = mcmc_approximation(mult, iterations=1e4)

var1 = test_accuracy(mult, mcmc_samples, "laplacian")

var2 = test_accuracy(mult, mcmc_samples, "gva")

var3 = test_accuracy(mult, mcmc_samples, "gva2")

var4 = test_accuracy(mult, mcmc_samples, "gva_nr")

}

}
test_accuracies()

test_accuracies_slope = function()

{

Monte Carlo Markov Chains approximation

seed = 1

set.seed(seed)

mult = generate_slope_test_data()

mcmc_samples = mcmc_approximation(mult, seed=seed, iterations=1e6, warmup=1e5)

save(mult, mcmc_samples, file="data/accuracy_slope_2015_03_03.RData")

load(file="data/accuracy_slope_2015_03_03.RData")

now = Sys.time()

var1 = test_accuracy(mult, mcmc_samples, "laplacian", plot=TRUE)

print(Sys.time() - now)

print(var1)

now = Sys.time()

var2 = test_accuracy(mult, mcmc_samples, "gva", plot=TRUE)

print(Sys.time() - now)

print(image(Matrix(var2$var_result$mLambda)))

print(var2)

now = Sys.time()

var3 = test_accuracy(mult, mcmc_samples, "gva2", plot=TRUE)

print(image(Matrix(var3$var_result$mLambda)))

print(Sys.time() - now)

print(var3)

print("gva2new")

now = Sys.time()

Rprof(line.profiling=TRUE)

var4 = test_accuracy(mult, mcmc_samples, "gva2new", plot=TRUE)

Rprof(NULL)

print(summaryRprof(lines = "both"))

print(image(Matrix(var4$var_result$mLambda)))

var4 = test_accuracy(mult, mcmc_samples, "gva2new")

print(Sys.time() - now)

print(var4)

now = Sys.time()

GVA NR is unstable, and sometimes fails with an error

var5 = test_accuracy(mult, mcmc_samples, "gva_nr")

var5 = test_accuracy(mult, mcmc_samples, "gva_nr", plot=TRUE)

print(Sys.time() - now)

print(image(Matrix(var5$var_result$mLambda)))

print(var5)

save(var1, var2, var3, var4, var5, file="accuracy_results_slope.RData")

print(image(Matrix(var3_new$mLambda)))

now = Sys.time()

var4 = test_accuracy(mult, mcmc_samples, "gva_nr")

Sys.time() - now

print(image(Matrix(var4$mLambda)))

mult = generate_spline_test_data()

mcmc_samples = mcmc_approximation(mult, iterations=1e4)

save(mult, mcmc_samples, file="accuracy_spline.RData")

now = Sys.time()

var1 = test_accuracy(mult, mcmc_samples, "laplacian")

Sys.time() - now

print(image(Matrix(var1$var_result$mLambda)))

var1

now = Sys.time()

var2 = test_accuracy(mult, mcmc_samples, "gva")

Sys.time() - now

print(image(Matrix(var2$var_result$mLambda)))

var2

now = Sys.time()

var3 = test_accuracy(mult, mcmc_samples, "gva2")

Sys.time() - now

print(image(Matrix(var3$var_result$mLambda)))

var3

now = Sys.time()

var4 = test_accuracy(mult, mcmc_samples, "gva2new")

Sys.time() - now

print(image(Matrix(var4$var_result$mLambda)))

var4

Rprof()

now = Sys.time()

var5 = test_accuracy(mult, mcmc_samples, "gva_nr")

print(Sys.time() - now)

print(image(Matrix(var5$var_result$mLambda)))

var5

} test_accuracies_slope() N 100 p 2 m 5 blocksize 2 spline_dim 0 calculate_lower_bound: T1 -4.10617e+17 T2 -105.5609 T3 -38.79288 a_sigma b_sigma calculate_lower_bound: T1 -4254.533 T2 -48.13332 T3 -47.93269 Iteration 2 : lower bound -4350.599 difference 4.10617e+17 parameters vmu 1.970248 2.387512 1.101093 -2.379948 1.825163 -2.837605 2.436981 -2.691774 -0.1999721 -1.972204 diag(mLambda) 0.008005845 0.00468602 0.01032378 0.007912696 0.009163322 0.006823373 0.008631102 0.005365371 0.01762683 0.01577659 a_rho 84.0028 b_rho 17.9972 a_sigma b_sigma calculate_lower_bound: T1 -4248.745 T2 -48.09777 T3 -47.92874 Iteration 3 : lower bound -4344.772 difference 5.826987 parameters vmu 2.093207 2.303359 0.9793739 -2.292995 1.701905 -2.752487 2.39622 -2.570463 -0.3002183 -1.871611 diag(mLambda) 0.007165469 0.004302308 0.009434513 0.007534457 0.008300133 0.006445907 0.007727148 0.004922677 0.01642191 0.01512188 a_rho 84.00333 b_rho 17.99667 a_sigma b_sigma calculate_lower_bound: T1 -4248.74 T2 -48.10542 T3 -47.92929 Iteration 4 : lower bound -4344.774 difference -0.002537742 parameters vmu 2.09571 2.301589 0.9766196 -2.291273 1.699333 -2.750677 2.393756 -2.568625 -0.3045866 -1.870705 diag(mLambda) 0.007154014 0.00429596 0.009425739 0.007527698 0.008289927 0.006438819 0.007716174 0.004915965 0.01642994 0.01513047 a_rho 84.00326 b_rho 17.99674 a_sigma b_sigma [1] 100 8 NULL Time difference of 3.77699 secs $var_result $var_result$vmu ,1 2.0957098 x 2.3015888 groups2 0.9766196 x:groups2 -2.2912729 groups3 1.6993328 x:groups3 -2.7506774 groups4 2.3937563 x:groups4 -2.5686248 groups5 -0.3045866 x:groups5 -1.8707046

$var_result$mLambda (Intercept) x groups2 x:groups2 groups3 (Intercept) 0.007154014 -0.005154323 -0.007136626 0.005163922 -0.007144532 x -0.005154323 0.004295960 0.005143068 -0.004301498 0.005148240 groups2 -0.007136626 0.005143068 0.009425739 -0.005162749 0.007127168 x:groups2 0.005163922 -0.004301498 -0.005162749 0.007527698 -0.005157824 groups3 -0.007144532 0.005148240 0.007127168 -0.005157824 0.008289927 x:groups3 0.005162642 -0.004300897 -0.005151365 0.004306450 -0.004903491 groups4 -0.007149408 0.005151363 0.007132031 -0.005160954 0.007139932 x:groups4 0.005156909 -0.004297504 -0.005145647 0.004303047 -0.005150822 groups5 -0.007097451 0.005116721 0.007080207 -0.005126237 0.007088048 x:groups5 0.005150843 -0.004291441 -0.005139592 0.004296980 -0.005144762 x:groups3 groups4 x:groups4 groups5 x:groups5 (Intercept) 0.005162642 -0.007149408 0.005156909 -0.007097451 0.005150843 x -0.004300897 0.005151363 -0.004297504 0.005116721 -0.004291441 groups2 -0.005151365 0.007132031 -0.005145647 0.007080207 -0.005139592 x:groups2 0.004306450 -0.005160954 0.004303047 -0.005126237 0.004296980 groups3 -0.004903491 0.007139932 -0.005150822 0.007088048 -0.005144762 x:groups3 0.006438819 -0.005159675 0.004302446 -0.005124969 0.004296379 groups4 -0.005159675 0.007716174 -0.005056923 0.007092883 -0.005147884 x:groups4 0.004302446 -0.005056923 0.004915965 -0.005119285 0.004292986 groups5 -0.005124969 0.007092883 -0.005119285 0.016429939 -0.009887394 x:groups5 0.004296379 -0.005147884 0.004292986 -0.009887394 0.015130470

$var_result$a_rho [1] 84.00326

$var_result$b_rho [1] 17.99674

$var_result$a_sigma numeric(0)

$var_result$b_sigma numeric(0)

$var_result$vlower_bound [1] -4.106170e+17 -4.350599e+03 -4.344772e+03 -4.344774e+03

$vbeta_accuracy [1] 0.9035163 0.9057978

$vu_accuracy [1] 0.9168333 0.9106049 0.9179837 0.9240022 0.9067588 0.9171403 0.9006352 [8] 0.8843505

$rho_accuracy [1] 0.9294158

N 100 p 2 m 5 blocksize 2 spline_dim 0 calculate_lower_bound: T1 -4290.507 T2 -49.87344 T3 -47.95629 a_sigma b_sigma calculate_lower_bound: T1 -4248.556 T2 -48.19419 T3 -47.92916 Iteration 2 : lower bound -4344.679 difference 43.6573 parameters vmu 2.09409 2.302492 0.974904 -2.291347 1.699599 -2.752365 2.394914 -2.569596 -0.3210042 -1.869959 diag(mLambda) 0.008481679 0.004919913 0.0112074 0.008822642 0.009748994 0.007528633 0.009148985 0.005618053 0.0210578 0.02150873 a_rho 84.00328 b_rho 17.99672 a_sigma b_sigma calculate_lower_bound: T1 -4248.571 T2 -48.15913 T3 -47.92901 Iteration 3 : lower bound -4344.66 difference 0.01940456 parameters vmu 2.093146 2.302883 0.9762685 -2.291192 1.700576 -2.752579 2.395646 -2.570119 -0.3150815 -1.866813 diag(mLambda) 0.008453228 0.004908299 0.01116755 0.008812885 0.009714915 0.00752008 0.009117964 0.005607909 0.02092631 0.02143282 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15807 T3 -47.92899 Iteration 4 : lower bound -4344.659 difference 0.0005529186 parameters vmu 2.093039 2.302952 0.976395 -2.291246 1.700687 -2.752644 2.395749 -2.570193 -0.3148003 -1.866767 diag(mLambda) 0.008453131 0.00490833 0.01116711 0.008812966 0.009714652 0.007520178 0.009117788 0.005607981 0.02092304 0.02143098 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 5 : lower bound -4344.659 difference 2.113186e-05 parameters vmu 2.093033 2.302956 0.9764016 -2.291249 1.700693 -2.752648 2.395754 -2.570197 -0.3147875 -1.866766 diag(mLambda) 0.008453145 0.004908341 0.01116711 0.008812988 0.009714659 0.007520194 0.009117799 0.005607995 0.02092289 0.02143083 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 6 : lower bound -4344.659 difference 5.468974e-07 parameters vmu 2.093032 2.302956 0.9764023 -2.291249 1.700694 -2.752648 2.395755 -2.570198 -0.3147866 -1.866767 diag(mLambda) 0.008453146 0.004908345 0.01116712 0.008812998 0.009714662 0.007520199 0.009117801 0.005607998 0.02092289 0.02143081 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 7 : lower bound -4344.659 difference 3.403011e-07 parameters vmu 2.093033 2.302956 0.9764017 -2.291249 1.700693 -2.752648 2.395754 -2.570197 -0.3147873 -1.866766 diag(mLambda) 0.008453147 0.004908342 0.01116711 0.008812986 0.009714661 0.007520193 0.009117799 0.005607994 0.02092288 0.02143082 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 8 : lower bound -4344.659 difference -1.964881e-07 parameters vmu 2.093032 2.302956 0.9764023 -2.291249 1.700694 -2.752648 2.395755 -2.570198 -0.3147865 -1.866767 diag(mLambda) 0.008453147 0.004908342 0.01116711 0.00881299 0.00971466 0.007520198 0.009117801 0.005607997 0.02092289 0.02143083 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma [1] 100 8 NULL Time difference of 3.73424 secs $var_result $var_result$vmu [1] 2.0930321 2.3029562 0.9764023 -2.2912495 1.7006940 -2.7526481 [7] 2.3957550 -2.5701976 -0.3147865 -1.8667666

$var_result$mLambda [,1] [,2] [,3] [,4] [,5] [1,] 0.008453147 -0.006044300 -0.008430261 0.006054039 -0.008441128 [2,] -0.006044300 0.004908342 0.006029261 -0.004913886 0.006036489 [3,] -0.008430261 0.006029261 0.011167115 -0.006484648 0.008418278 [4,] 0.006054039 -0.004913886 -0.006484648 0.008812990 -0.006046212 [5,] -0.008441128 0.006036489 0.008418278 -0.006046212 0.009714660 [6,] 0.006055017 -0.004914797 -0.006039947 0.004920356 -0.005868732 [7,] -0.008446602 0.006040059 0.008423736 -0.006049789 0.008434593 [8,] 0.006048102 -0.004910679 -0.006033052 0.004916228 -0.006040285 [9,] -0.008360108 0.005982405 0.008337484 -0.005992025 0.008348227 [10,] 0.006054225 -0.004911368 -0.006039150 0.004916931 -0.006046394 [,6] [,7] [,8] [,9] [,10] [1,] 0.006055017 -0.008446602 0.006048102 -0.008360108 0.006054225 [2,] -0.004914797 0.006040059 -0.004910679 0.005982405 -0.004911368 [3,] -0.006039947 0.008423736 -0.006033052 0.008337484 -0.006039150 [4,] 0.004920356 -0.006049789 0.004916228 -0.005992025 0.004916931 [5,] -0.005868732 0.008434593 -0.006040285 0.008348227 -0.006046394 [6,] 0.007520198 -0.006050766 0.004917140 -0.005992994 0.004917842 [7,] -0.006050766 0.009117801 -0.005890535 0.008353638 -0.006049972 [8,] 0.004917140 -0.005890535 0.005607997 -0.005986161 0.004913711 [9,] -0.005992994 0.008353638 -0.005986161 0.020922894 -0.011423447 [10,] 0.004917842 -0.006049972 0.004913711 -0.011423447 0.021430829

$var_result$a_rho [1] 84.0033

$var_result$b_rho [1] 17.9967

$var_result$a_sigma numeric(0)

$var_result$b_sigma numeric(0)

$var_result$vlower_bound [1] -4388.336 -4344.679 -4344.660 -4344.659 -4344.659 -4344.659 -4344.659 [8] -4344.659

$vbeta_accuracy [1] 0.9308112 0.9251300

$vu_accuracy [1] 0.9332386 0.9255522 0.9373222 0.9434991 0.9286187 0.9388787 0.9601703 [8] 0.9324045

$rho_accuracy [1] 0.9294196

N 100 p 2 m 5 blocksize 2 spline_dim 0 Error in optim(par = vmu, fn = f.GVA_new, gr = vg.GVA_new, method = method, : non-finite value supplied by optim Calls: test_accuracies_slope ... zero_infl_var -> zero_infl_var.multivariate -> fit.GVA_new -> optim Execution halted

JohnOrmerod commented 9 years ago

argh!

From: Mark Greenaway [notifications@github.com] Sent: Wednesday, March 04, 2015 3:36 PM To: certifiedwaif/phd Subject: [phd] All of the algorithms run on R 3.1.0 on my Mac, but fail on R 3.1.2 on the Linux machines at Uni. (#8)

This is very odd.

Mac log:

R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details.

Natural language support but running in an English locale

R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R.

accuracy.R

source("zero_inflated_model.R") source("generate.R") source("mcmc.R") source("rwmh.R")

calculate_accuracies = function(mult, mcmc_samples, var_result, approximation, print_flag=FALSE, plot_flag=FALSE)

{
TODO: Add support for checking the accuracy over multiple dimensions
cubature$adaptIntegrate

*
if (plot_flag) pdf(paste0("accuracyplots", approximation, ".pdf"))
return(var_result)
vbeta accuracy
calculate_accuracy = function(mcmc_samples, dist_fn, param1, param2)
{
mcmc_density = density(mcmc_samples)
mcmc_fn = splinefun(mcmc_density$x, mcmc_density$y) *
integrand <- function(x)
{
return(abs(mcmc_fn(x) - dist_fn(x, param1, param2)))
}
result = integrate(integrand, min(mcmc_density$x), max(mcmc_density$x),
subdivisions = length(mcmc_density$x))
accuracy = 1 - .5 * result$value
return(accuracy)
} *
Compare MCMC distribution with variational approximation for each parameter
vnu[i] ~ Normal, dnorm
sigma2_u ~ IG, dgamma(1/x)
rho ~ Beta, dbeta
vr[i] ~ Bernoulli, dbinom
For each parameter of interest,
* estimate density of MCMC
* compare with q distribution using L_1 norm

*
Kernel density estimates of MCMC-estimated posteriors
Use L_1 distance to compare against variational approximations of posteriors

*
accuracy_plot = function(mcmc_samples, dist_fn, param1, param2)
{
mcmc_density = density(mcmc_samples)
plot(mcmc_density)
curve(dist_fn(x, param1, param2),
from=min(mcmc_density$x), to=max(mcmc_density$x),
add=TRUE, lty=2, col="blue")
} *
vbeta_accuracy = rep(NA, ncol(mult$mX))
for (i in 1:ncol(mult$mX)) {
vbeta_accuracy[i] = calculate_accuracy(mcmc_samples$vbeta[,i], dnorm,
var_result$vmu[i], sqrt(var_result$mLambda[i,i]))
if (print_flag) cat("vbeta[", i, "]", approximation, "accuracy:", vbeta_accuracy[i], "\n") *
param_name = sprintf("vbeta[%d]", i)
if (plot_flag) accuracy_plot(mcmc_samples$vbeta[,i], dnorm,
var_result$vmu[i], sqrt(var_result$mLambda[i,i]))
} *
vu accuracy
FIXME: To check for random slopes accuracy, this section will have
to get more complex.
print(dim(mult$mZ))
print(dim(mcmc_samples$u))
vu_accuracy = rep(NA, ncol(mult$mZ))
B = mult$blocksize
b_idx = 1
for (i in 1:ncol(mult$mZ)) {
TODO: The B - (i %% B) expression only works for B=2. Rewrite this.
m_idx = ceiling(i/B)
vu_accuracy[i] = calculate_accuracy(mcmc_samples$vu[,m_idx,b_idx], dnorm,
var_result$vmu[i+mult$p], sqrt(var_result$mLambda[i+mult$p,i+mult$p]))
if (print_flag) cat("vu[", i, "]", approximation, "accuracy:", vu_accuracy[i], "\n")
if (plot_flag) accuracy_plot(mcmc_samples$vu[,m_idx,b_idx], dnorm,
var_result$vmu[i+mult$p], sqrt(var_result$mLambda[i+mult$p,i+mult$p])) *
b_idx = b_idx + 1
if (b_idx > B)
b_idx=1
} *
sigma2_u accuracy
FIXME - this may be wrong for blocksize != 1?
This is totally wrong for the Inverse Wishart model?
sigma2_u_accuracy = calculate_accuracy(1/mcmc_samples$sigma_u^2, dgamma,
var_result$a_sigma, var_result$b_sigma)
if (print_flag) cat("sigma2_u", approximation, "accuracy:", sigma2_u_accuracy, "\n")
if (plot_flag) accuracy_plot(1/mcmc_samples$sigma_u^2, dgamma,
var_result$a_sigma, var_result$b_sigma)

*
rho accuracy
rho_accuracy = calculate_accuracy(mcmc_samples$rho, dbeta,
var_result$a_rho, var_result$b_rho)
if (print_flag) cat("rho", approximation, "accuracy: ", rho_accuracy, "\n")
if (plot_flag) accuracy_plot(mcmc_samples$rho, dbeta,
var_result$a_rho, var_result$b_rho)
if (plot_flag) dev.off()
return(list(var_result=var_result,
vbeta_accuracy=vbeta_accuracy,
vu_accuracy=vu_accuracy,
sigma2_u_accuracy=sigma2_u_accuracy,
rho_accuracy=rho_accuracy))
}

test_accuracy = function(mult, mcmc_samples, approximation, plot=FALSE)

{
var_result = zero_infl_var(mult, method=approximation, verbose=TRUE)
return(calculate_accuracies(mult, mcmc_samples, var_result, approximation, plot_flag=plot))
}

Calculate accuracy ---- Approximate the L1 norm between the variational approximation and the MCMC approximation

calculate_univariate_accuracy <- function(result_mcmc, var_result)

{
density_mcmc_rho = density(result_mcmc$vrho)
integrand <- function(x)
{
fn = splinefun(density_mcmc_rho$x, density_mcmc_rho$y)
return(abs(fn(x) - dbeta(x, var_result$a_rho, var_result$b_rho)))
}
result = integrate(integrand, min(density_mcmc_rho$x), max(density_mcmc_rho$x), subdivisions = length(density_mcmc_rho$x))
rho_accuracy = 1 - .5 * result$value *
density_mcmc_lambda = density(result_mcmc$vlambda)
integrand <- function(x)
{
fn = splinefun(density_mcmc_lambda$x, density_mcmc_lambda$y)
return(abs(fn(x) - dgamma(x, var_result$a_lambda, var_result$b_lambda)))
}
result = integrate(integrand, min(density_mcmc_lambda$x), max(density_mcmc_lambda$x), subdivisions = length(density_mcmc_lambda$x))
lambda_accuracy = 1 - .5 * result$value
return(list(rho_accuracy=rho_accuracy, lambda_accuracy=lambda_accuracy))
}

test_accuracies = function()

{
Need to be able to compare the solution paths of each approximation

*
Generate data
for (i in 1:100) {
set.seed(i)
mult = generate_test_data(20, 100)
Monte Carlo Markov Chains approximation
mcmc_samples = mcmc_approximation(mult, iterations=1e6)
Save the results, because this takes such a long time to run.
}
save(mult, mcmc_samples, file="accuracy_good.RData")
set.seed(1)
mult = generate_test_data(10, 100)
Monte Carlo Markov Chains approximation
mcmc_samples = mcmc_approximation(mult, iterations=1e6, warmup = 1e4)
Save the results, because this takes such a long time to run.
save(mult, mcmc_samples, file="accuracy.RData")
save(mult, mcmc_samples, file="data/accuracy_int.RData")
load(file="data/accuracy_int_2015_02_17.RData")
mult$spline_dim = 0
load(file="accuracy.RData")
Test all other approximations against it
load(file="accuracy.RData")

*
Test multivariate approximation's accuracy
now = Sys.time()
var1 = test_accuracy(mult, mcmc_samples, "laplacian")
print(Sys.time() - now)
print(image(Matrix(var1$var_result$mLambda)))
print(var1) *
now = Sys.time()
var2 = test_accuracy(mult, mcmc_samples, "gva")
print(Sys.time() - now)
print(image(Matrix(var2$var_result$mLambda)))
print(var2) *
now = Sys.time()
var3 = test_accuracy(mult, mcmc_samples, "gva2", plot=TRUE)
print(Sys.time() - now)
print(image(Matrix(var3$mLambda)))
print(var3) *
Rprof()
now = Sys.time()
var4 = test_accuracy(mult, mcmc_samples, "gva2new")
print(Sys.time() - now)
print(image(Matrix(var4$var_result$mLambda)))
print(var4) *
Rprof(NULL)
print(summaryRprof())
print(image(Matrix(var3_new$mLambda)))

*
now = Sys.time()
var5 = test_accuracy(mult, mcmc_samples, "gva_nr")
print(Sys.time() - now)
print(image(Matrix(var4$result_var$mLambda)))
print(var5) *
save(var1, var2, var3, var4, var5, file="accuracy_results_int.RData")
for (i in 1:100) {
set.seed(i)
mult = generate_test_data(20, 100)
mcmc_samples = mcmc_approximation(mult, iterations=1e4)
#
var1 = test_accuracy(mult, mcmc_samples, "laplacian")
var2 = test_accuracy(mult, mcmc_samples, "gva")
var3 = test_accuracy(mult, mcmc_samples, "gva2")
var4 = test_accuracy(mult, mcmc_samples, "gva_nr")
}

*
} #test_accuracies()

test_accuracies_slope = function()

{
Monte Carlo Markov Chains approximation
set.seed(1)
mult = generate_slope_test_data()
mcmc_samples = mcmc_approximation(mult, iterations=1e6, warmup=1e4, mc.cores = 1)
save(mult, mcmc_samples, file="data/accuracy_slope_2015_03_03.RData")
load(file="data/accuracy_slope_2015_03_03.RData") *
now = Sys.time()
var1 = test_accuracy(mult, mcmc_samples, "laplacian", plot=TRUE)
print(Sys.time() - now)
print(var1) *
now = Sys.time()
var2 = test_accuracy(mult, mcmc_samples, "gva", plot=TRUE)
print(Sys.time() - now)
print(image(Matrix(var2$var_result$mLambda)))
print(var2) *
now = Sys.time()
var3 = test_accuracy(mult, mcmc_samples, "gva2", plot=TRUE)
print(image(Matrix(var3$var_result$mLambda)))
print(Sys.time() - now)
print(var3) *
print("gva2new")
now = Sys.time()
Rprof(line.profiling=TRUE)
var4 = test_accuracy(mult, mcmc_samples, "gva2new", plot=TRUE)
Rprof(NULL)
print(summaryRprof(lines = "both"))
print(image(Matrix(var4$var_result$mLambda)))
var4 = test_accuracy(mult, mcmc_samples, "gva2new")
print(Sys.time() - now)
print(var4) *
now = Sys.time()
GVA NR is unstable, and sometimes fails with an error
var5 = test_accuracy(mult, mcmc_samples, "gva_nr")
var5 = test_accuracy(mult, mcmc_samples, "gva_nr", plot=TRUE)
print(Sys.time() - now)
print(image(Matrix(var5$var_result$mLambda)))
print(var5) *
save(var1, var2, var3, var4, var5, file="accuracy_results_slope.RData")

*
print(image(Matrix(var3_new$mLambda)))

*
now = Sys.time()
var4 = test_accuracy(mult, mcmc_samples, "gva_nr")
Sys.time() - now
print(image(Matrix(var4$mLambda)))

*
mult = generate_spline_test_data()
mcmc_samples = mcmc_approximation(mult, iterations=1e4)
save(mult, mcmc_samples, file="accuracy_spline.RData")

*
now = Sys.time()
var1 = test_accuracy(mult, mcmc_samples, "laplacian")
Sys.time() - now
print(image(Matrix(var1$var_result$mLambda)))
var1

*
now = Sys.time()
var2 = test_accuracy(mult, mcmc_samples, "gva")
Sys.time() - now
print(image(Matrix(var2$var_result$mLambda)))
var2

*
now = Sys.time()
var3 = test_accuracy(mult, mcmc_samples, "gva2")
Sys.time() - now
print(image(Matrix(var3$var_result$mLambda)))
var3

*
now = Sys.time()
var4 = test_accuracy(mult, mcmc_samples, "gva2new")
Sys.time() - now
print(image(Matrix(var4$var_result$mLambda)))
var4

*
Rprof()
now = Sys.time()
var5 = test_accuracy(mult, mcmc_samples, "gva_nr")
print(Sys.time() - now)
print(image(Matrix(var5$var_result$mLambda)))
var5
} test_accuracies_slope() N 100 p 2 m 5 blocksize 2 spline_dim 0 calculate_lower_bound: T1 -4.10617e+17 T2 -105.5615 T3 -38.79287 a_sigma b_sigma calculate_lower_bound: T1 -4254.533 T2 -48.13332 T3 -47.93269 Iteration 2 : lower bound -4350.599 difference 4.10617e+17 parameters vmu 1.970248 2.387512 1.101093 -2.379948 1.825163 -2.837605 2.43698 -2.691774 -0.1999721 -1.972204 diag(mLambda) 0.008005845 0.00468602 0.01032378 0.007912696 0.009163322 0.006823373 0.008631102 0.005365371 0.01762683 0.01577659 a_rho 84.0028 b_rho 17.9972 a_sigma b_sigma calculate_lower_bound: T1 -4248.745 T2 -48.09777 T3 -47.92874 Iteration 3 : lower bound -4344.772 difference 5.827025 parameters vmu 2.093207 2.303359 0.9793739 -2.292995 1.701905 -2.752487 2.39622 -2.570463 -0.3002183 -1.871611 diag(mLambda) 0.007165469 0.004302308 0.009434513 0.007534457 0.008300133 0.006445907 0.007727148 0.004922677 0.01642191 0.01512188 a_rho 84.00333 b_rho 17.99667 a_sigma b_sigma calculate_lower_bound: T1 -4248.74 T2 -48.10542 T3 -47.92929 Iteration 4 : lower bound -4344.774 difference -0.002537746 parameters vmu 2.09571 2.301589 0.9766196 -2.291273 1.699333 -2.750677 2.393756 -2.568625 -0.3045866 -1.870705 diag(mLambda) 0.007154014 0.00429596 0.009425739 0.007527698 0.008289927 0.006438819 0.007716174 0.004915965 0.01642994 0.01513047 a_rho 84.00326 b_rho 17.99674 a_sigma b_sigma [1] 100 8 NULL Time difference of 2.466712 secs $var_result $var_result$vmu ,1 2.0957098 x 2.3015888 groups2 0.9766196 x:groups2 -2.2912729 groups3 1.6993328 x:groups3 -2.7506774 groups4 2.3937563 x:groups4 -2.5686248 groups5 -0.3045866 x:groups5 -1.8707046

$var_result$mLambda (Intercept) x groups2 x:groups2 groups3 (Intercept) 0.007154014 -0.005154323 -0.007136626 0.005163922 -0.007144532 x -0.005154323 0.004295960 0.005143068 -0.004301498 0.005148240 groups2 -0.007136626 0.005143068 0.009425739 -0.005162749 0.007127168 x:groups2 0.005163922 -0.004301498 -0.005162749 0.007527698 -0.005157824 groups3 -0.007144532 0.005148240 0.007127168 -0.005157824 0.008289927 x:groups3 0.005162642 -0.004300897 -0.005151365 0.004306450 -0.004903491 groups4 -0.007149408 0.005151363 0.007132031 -0.005160954 0.007139932 x:groups4 0.005156909 -0.004297504 -0.005145647 0.004303047 -0.005150822 groups5 -0.007097451 0.005116721 0.007080207 -0.005126237 0.007088048 x:groups5 0.005150843 -0.004291441 -0.005139592 0.004296980 -0.005144762 x:groups3 groups4 x:groups4 groups5 x:groups5 (Intercept) 0.005162642 -0.007149408 0.005156909 -0.007097451 0.005150843 x -0.004300897 0.005151363 -0.004297504 0.005116721 -0.004291441 groups2 -0.005151365 0.007132031 -0.005145647 0.007080207 -0.005139592 x:groups2 0.004306450 -0.005160954 0.004303047 -0.005126237 0.004296980 groups3 -0.004903491 0.007139932 -0.005150822 0.007088048 -0.005144762 x:groups3 0.006438819 -0.005159675 0.004302446 -0.005124969 0.004296379 groups4 -0.005159675 0.007716174 -0.005056923 0.007092883 -0.005147884 x:groups4 0.004302446 -0.005056923 0.004915965 -0.005119285 0.004292986 groups5 -0.005124969 0.007092883 -0.005119285 0.016429939 -0.009887394 x:groups5 0.004296379 -0.005147884 0.004292986 -0.009887394 0.015130470

$var_result$a_rho [1] 84.00326

$var_result$b_rho [1] 17.99674

$var_result$a_sigma numeric(0)

$var_result$b_sigma numeric(0)

$var_result$vlower_bound [1] -4.106170e+17 -4.350599e+03 -4.344772e+03 -4.344774e+03

$vbeta_accuracy [1] 0.9070440 0.8987978

$vu_accuracy [1] 0.9472439 0.8896154 0.9140297 0.9021092 0.9046075 0.9117498 0.8924308 [8] 0.8997999

$rho_accuracy [1] 0.9036479

N 100 p 2 m 5 blocksize 2 spline_dim 0 calculate_lower_bound: T1 -4290.5 T2 -49.87337 T3 -47.95629 a_sigma b_sigma calculate_lower_bound: T1 -4248.556 T2 -48.19411 T3 -47.92917 Iteration 2 : lower bound -4344.679 difference 43.65032 parameters vmu 2.09413 2.302459 0.9748742 -2.291327 1.699557 -2.752339 2.394875 -2.569562 -0.3210693 -1.869889 diag(mLambda) 0.008482581 0.004920235 0.01120775 0.00882215 0.009749764 0.007525693 0.009149858 0.005617782 0.02105552 0.02150821 a_rho 84.00327 b_rho 17.99673 a_sigma b_sigma calculate_lower_bound: T1 -4248.571 T2 -48.15914 T3 -47.92901 Iteration 3 : lower bound -4344.66 difference 0.01938413 parameters vmu 2.093147 2.302883 0.9762678 -2.291192 1.700575 -2.752579 2.395646 -2.570119 -0.315083 -1.866814 diag(mLambda) 0.008453674 0.004908523 0.01116801 0.008813112 0.009715358 0.007520294 0.009118412 0.005608128 0.02092709 0.0214335 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15807 T3 -47.92899 Iteration 4 : lower bound -4344.659 difference 0.000556743 parameters vmu 2.093039 2.302952 0.976395 -2.291246 1.700687 -2.752644 2.395749 -2.570193 -0.3148 -1.866767 diag(mLambda) 0.008453144 0.00490834 0.01116712 0.008812989 0.009714664 0.007520192 0.009117801 0.005607992 0.020923 0.02143089 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 5 : lower bound -4344.659 difference 2.094874e-05 parameters vmu 2.093033 2.302956 0.9764014 -2.291249 1.700693 -2.752648 2.395754 -2.570197 -0.314788 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812984 0.009714661 0.007520192 0.0091178 0.005607994 0.0209229 0.02143087 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 6 : lower bound -4344.659 difference 1.130324e-06 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395755 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453147 0.004908342 0.01116711 0.008812985 0.009714661 0.007520193 0.009117801 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 7 : lower bound -4344.659 difference -8.684765e-08 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395754 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453147 0.004908343 0.01116711 0.008812991 0.009714661 0.007520196 0.0091178 0.005607996 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma [1] 100 8 NULL Time difference of 1.995439 secs $var_result $var_result$vmu [1] 2.0930326 2.3029559 0.9764019 -2.2912491 1.7006935 -2.7526478 [7] 2.3957545 -2.5701973 -0.3147871 -1.8667661

$var_result$mLambda [,1] [,2] [,3] [,4] [,5] [1,] 0.008453147 -0.006044301 -0.008430261 0.006054040 -0.008441129 [2,] -0.006044301 0.004908343 0.006029263 -0.004913887 0.006036490 [3,] -0.008430261 0.006029263 0.011167114 -0.006484650 0.008418278 [4,] 0.006054040 -0.004913887 -0.006484650 0.008812991 -0.006046214 [5,] -0.008441129 0.006036490 0.008418278 -0.006046214 0.009714661 [6,] 0.006055018 -0.004914797 -0.006039948 0.004920357 -0.005868733 [7,] -0.008446603 0.006040060 0.008423736 -0.006049790 0.008434595 [8,] 0.006048103 -0.004910680 -0.006033053 0.004916230 -0.006040286 [9,] -0.008360108 0.005982407 0.008337485 -0.005992026 0.008348228 [10,] 0.006054226 -0.004911368 -0.006039150 0.004916932 -0.006046394 [,6] [,7] [,8] [,9] [,10] [1,] 0.006055018 -0.008446603 0.006048103 -0.008360108 0.006054226 [2,] -0.004914797 0.006040060 -0.004910680 0.005982407 -0.004911368 [3,] -0.006039948 0.008423736 -0.006033053 0.008337485 -0.006039150 [4,] 0.004920357 -0.006049790 0.004916230 -0.005992026 0.004916932 [5,] -0.005868733 0.008434595 -0.006040286 0.008348228 -0.006046394 [6,] 0.007520196 -0.006050767 0.004917140 -0.005992995 0.004917842 [7,] -0.006050767 0.009117800 -0.005890535 0.008353639 -0.006049972 [8,] 0.004917140 -0.005890535 0.005607996 -0.005986162 0.004913711 [9,] -0.005992995 0.008353639 -0.005986162 0.020922893 -0.011423462 [10,] 0.004917842 -0.006049972 0.004913711 -0.011423462 0.021430839

$var_result$a_rho [1] 84.0033

$var_result$b_rho [1] 17.9967

$var_result$a_sigma numeric(0)

$var_result$b_sigma numeric(0)

$var_result$vlower_bound [1] -4388.329 -4344.679 -4344.660 -4344.659 -4344.659 -4344.659 -4344.659

$vbeta_accuracy [1] 0.9218658 0.9196421

$vu_accuracy [1] 0.9555139 0.9244812 0.9300057 0.9253950 0.9257944 0.9294784 0.9325449 [8] 0.9585504

$rho_accuracy [1] 0.9036515

N 100 p 2 m 5 blocksize 2 spline_dim 0 calculate_lower_bound: T1 -4.10617e+17 T2 -105.5615 T3 -38.79287 a_sigma b_sigma calculate_lower_bound: T1 -4254.581 T2 -48.17991 T3 -47.93239 Iteration 2 : lower bound -4350.693 difference 4.10617e+17 parameters vmu 1.967162 2.38913 1.101495 -2.380752 1.826973 -2.839919 2.439368 -2.693567 -0.2090716 -1.967071 diag(mLambda) 0.009626389 0.005505333 0.0123472 0.009262223 0.01084605 0.008085332 0.01031284 0.006211642 0.02329809 0.02103267 a_rho 84.00284 b_rho 17.99716 a_sigma b_sigma calculate_lower_bound: T1 -4248.666 T2 -48.1396 T3 -47.92847 Iteration 3 : lower bound -4344.734 difference 5.959453 parameters vmu 2.089955 2.305214 0.9802474 -2.2943 1.70389 -2.755052 2.398955 -2.572443 -0.3080842 -1.8672 diag(mLambda) 0.009357424 0.005479383 0.01210666 0.009380454 0.01067038 0.007993289 0.009852098 0.006041833 0.02253293 0.0211349 a_rho 84.00337 b_rho 17.99663 a_sigma b_sigma calculate_lower_bound: T1 -4248.657 T2 -48.15205 T3 -47.92845 Iteration 4 : lower bound -4344.738 difference -0.004210441 parameters vmu 2.089665 2.305241 0.9794486 -2.294553 1.704081 -2.754486 2.399184 -2.57229 -0.3092569 -1.868841 diag(mLambda) 0.009506722 0.005478045 0.01233618 0.009400417 0.01077984 0.008051188 0.009967946 0.006027864 0.02250331 0.02121321 a_rho 84.00337 b_rho 17.99663 a_sigma b_sigma [1] 100 8 NULL Time difference of 2.084503 secs $var_result $var_result$vmu [1] 2.0896655 2.3052413 0.9794486 -2.2945533 1.7040812 -2.7544864 [7] 2.3991839 -2.5722897 -0.3092569 -1.8688408

$var_result$mLambda (Intercept) x groups2 x:groups2 groups3 (Intercept) 0.009506722 -0.006810024 -0.009693908 0.007077483 -0.009543433 x -0.006810024 0.005478045 0.006940870 -0.005700327 0.006831581 groups2 -0.009693908 0.006940870 0.012336184 -0.007224387 0.009722892 x:groups2 0.007077483 -0.005700327 -0.007224387 0.009400417 -0.007097980 groups3 -0.009543433 0.006831581 0.009722892 -0.007097980 0.010779844 x:groups3 0.006991532 -0.005632224 -0.007123440 0.005849564 -0.006738496 groups4 -0.009445196 0.006760559 0.009622777 -0.007024161 0.009474876 x:groups4 0.006747138 -0.005439216 -0.006875162 0.005649940 -0.006767274 groups5 -0.010435842 0.007494819 0.010619410 -0.007773163 0.010468587 x:groups5 0.008099759 -0.006516700 -0.008236823 0.006750237 -0.008124073 x:groups3 groups4 x:groups4 groups5 x:groups5 (Intercept) 0.006991532 -0.009445196 0.006747138 -0.010435842 0.008099759 x -0.005632224 0.006760559 -0.005439216 0.007494819 -0.006516700 groups2 -0.007123440 0.009622777 -0.006875162 0.010619410 -0.008236823 x:groups2 0.005849564 -0.007024161 0.005649940 -0.007773163 0.006750237 groups3 -0.006738496 0.009474876 -0.006767274 0.010468587 -0.008124073 x:groups3 0.008051188 -0.006940005 0.005583984 -0.007687403 0.006681336 groups4 -0.006940005 0.009967946 -0.006592944 0.010364304 -0.008048975 x:groups4 0.005583984 -0.006592944 0.006027864 -0.007419028 0.006467086 groups5 -0.007687403 0.010364304 -0.007419028 0.022503309 -0.015027962 x:groups5 0.006681336 -0.008048975 0.006467086 -0.015027962 0.021213209

$var_result$a_rho [1] 84.00337

$var_result$b_rho [1] 17.99663

$var_result$a_sigma numeric(0)

$var_result$b_sigma numeric(0)

$var_result$vlower_bound [1] -4.106170e+17 -4.350693e+03 -4.344734e+03 -4.344738e+03

$vbeta_accuracy [1] 0.9223178 0.9321255

$vu_accuracy [1] 0.9530949 0.9346686 0.9326032 0.9323161 0.9308158 0.9364944 0.9189237 [8] 0.9602546

$rho_accuracy [1] 0.9036579

[1] "gva2new" N 100 p 2 m 5 blocksize 2 spline_dim 0 calculate_lower_bound: T1 -4292.737 T2 -49.86207 T3 -47.95629 a_sigma b_sigma calculate_lower_bound: T1 -4250.696 T2 -48.16968 T3 -47.9283 Iteration 2 : lower bound -4346.794 difference 43.76168 parameters vmu 2.097054 2.300623 0.9717529 -2.289799 1.696478 -2.750359 2.391713 -2.567674 -0.3238487 -1.869233 diag(mLambda) 0.0002772039 0.0002585945 0.002714801 0.003825567 0.001274922 0.002529919 0.0006452087 0.0006648273 0.01206992 0.01604501 a_rho 84.00339 b_rho 17.99661 a_sigma b_sigma calculate_lower_bound: T1 -4250.713 T2 -48.13132 T3 -47.92848 Iteration 3 : lower bound -4346.773 difference 0.02088818 parameters vmu 2.097566 2.299943 0.9716072 -2.288637 1.695939 -2.749491 2.390994 -2.567088 -0.3199915 -1.86376 diag(mLambda) 0.0002773868 0.0002585225 0.002711821 0.003823748 0.001274627 0.002565101 0.0006452489 0.0006651087 0.01199677 0.01599983 a_rho 84.00337 b_rho 17.99663 a_sigma b_sigma calculate_lower_bound: T1 -4250.709 T2 -48.13713 T3 -47.92777 Iteration 4 : lower bound -4346.774 difference -0.001010009 parameters vmu 2.094479 2.302267 0.9748108 -2.291137 1.6991 -2.751906 2.394069 -2.569463 -0.3155559 -1.868097 diag(mLambda) 0.0002774194 0.0002585368 0.00271153 0.003823654 0.001274604 0.002568743 0.0006452575 0.0006651432 0.01199424 0.01599953 a_rho 84.00347 b_rho 17.99653 a_sigma b_sigma [1] 100 8 NULL Time difference of 2.075891 secs $var_result $var_result$vmu [1] 2.0944786 2.3022674 0.9748108 -2.2911367 1.6990997 -2.7519064 [7] 2.3940692 -2.5694626 -0.3155559 -1.8680974

$var_result$mLambda 1 2 1 2 3 1 2.774194e-04 1.865033e-06 -2.463851e-05 -5.330677e-06 -1.859036e-05 2 1.865033e-06 2.585368e-04 -3.162288e-06 -2.158638e-05 -7.083789e-07 1 -2.463851e-05 -3.162288e-06 2.711530e-03 7.232348e-07 1.657832e-06 2 -5.330677e-06 -2.158638e-05 7.232348e-07 3.823654e-03 4.058501e-07 3 -1.859036e-05 -7.083789e-07 1.657832e-06 4.058501e-07 1.274604e-03 4 -2.425131e-06 -1.951697e-05 4.414230e-07 1.672172e-06 2.065185e-07 5 -1.106122e-05 7.271966e-06 8.972287e-07 -3.998449e-07 7.246534e-07 6 7.312536e-06 -3.111613e-08 -6.485195e-07 -1.338202e-07 -4.898446e-07 7 -4.127052e-05 9.158416e-07 3.651536e-06 6.935497e-07 2.762916e-06 8 -2.072022e-05 -4.610888e-05 2.373080e-06 4.230161e-06 1.492236e-06 4 5 6 7 8 1 -2.425131e-06 -1.106122e-05 7.312536e-06 -4.127052e-05 -2.072022e-05 2 -1.951697e-05 7.271966e-06 -3.111613e-08 9.158416e-07 -4.610888e-05 1 4.414230e-07 8.972287e-07 -6.485195e-07 3.651536e-06 2.373080e-06 2 1.672172e-06 -3.998449e-07 -1.338202e-07 6.935497e-07 4.230161e-06 3 2.065185e-07 7.246534e-07 -4.898446e-07 2.762916e-06 1.492236e-06 4 2.568743e-03 -4.574443e-07 -5.786903e-08 2.707655e-07 3.648649e-06 5 -4.574443e-07 6.452575e-04 2.125136e-05 1.679440e-06 -4.801365e-07 6 -5.786903e-08 2.125136e-05 6.651432e-04 -1.088226e-06 -5.318928e-07 7 2.707655e-07 1.679440e-06 -1.088226e-06 1.199424e-02 -4.195596e-03 8 3.648649e-06 -4.801365e-07 -5.318928e-07 -4.195596e-03 1.599953e-02

$var_result$a_rho [1] 84.00347

$var_result$b_rho [1] 17.99653

$var_result$a_sigma numeric(0)

$var_result$b_sigma numeric(0)

$var_result$vlower_bound [1] -4390.555 -4346.794 -4346.773 -4346.774

$vbeta_accuracy [1] 0.3031139 0.3435177

$vu_accuracy [1] 0.6712743 0.7325994 0.5152712 0.6898212 0.4117733 0.4804914 0.8437027 [8] 0.9131477

$rho_accuracy [1] 0.9036661

N 100 p 2 m 5 blocksize 2 spline_dim 0 calculate_lower_bound: T1 -4290.483 T2 -49.87728 T3 -47.95629 a_sigma b_sigma calculate_lower_bound: T1 -4248.556 T2 -48.19414 T3 -47.92916 Iteration 2 : lower bound -4344.679 difference 43.63771 parameters vmu 2.094092 2.30249 0.9749076 -2.291343 1.699598 -2.752365 2.394914 -2.569594 -0.3210242 -1.869919 diag(mLambda) 0.008481814 0.004919978 0.01120775 0.008823185 0.009748932 0.007530014 0.009149182 0.005618411 0.02106534 0.02151484 a_rho 84.00328 b_rho 17.99672 a_sigma b_sigma calculate_lower_bound: T1 -4248.571 T2 -48.15913 T3 -47.92901 Iteration 3 : lower bound -4344.66 difference 0.01938018 parameters vmu 2.093146 2.302883 0.9762684 -2.291192 1.700575 -2.752579 2.395646 -2.570119 -0.3150817 -1.866813 diag(mLambda) 0.008453264 0.004908316 0.01116759 0.008812906 0.009714952 0.007520098 0.009118001 0.005607927 0.02092649 0.02143308 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15807 T3 -47.92899 Iteration 4 : lower bound -4344.659 difference 0.00055306 parameters vmu 2.093039 2.302952 0.9763949 -2.291246 1.700687 -2.752644 2.395749 -2.570193 -0.3148001 -1.866767 diag(mLambda) 0.008453131 0.004908331 0.01116711 0.008812973 0.009714652 0.007520181 0.009117788 0.005607983 0.02092299 0.02143091 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 5 : lower bound -4344.659 difference 2.0825e-05 parameters vmu 2.093033 2.302956 0.9764015 -2.291249 1.700693 -2.752648 2.395754 -2.570197 -0.3147877 -1.866766 diag(mLambda) 0.008453145 0.004908341 0.01116711 0.008812985 0.00971466 0.007520193 0.009117799 0.005607994 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 6 : lower bound -4344.659 difference 8.548432e-07 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395754 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812986 0.009714661 0.007520194 0.0091178 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 7 : lower bound -4344.659 difference 3.64098e-08 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395755 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812986 0.009714661 0.007520194 0.0091178 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 8 : lower bound -4344.659 difference 1.578883e-09 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395755 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812986 0.009714661 0.007520194 0.0091178 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 9 : lower bound -4344.659 difference 6.548362e-11 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395755 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812986 0.009714661 0.007520194 0.0091178 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 10 : lower bound -4344.659 difference 2.728484e-12 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395755 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812986 0.009714661 0.007520194 0.0091178 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 11 : lower bound -4344.659 difference 1.818989e-12 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395755 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812986 0.009714661 0.007520194 0.0091178 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 12 : lower bound -4344.659 difference 0 parameters vmu 2.093033 2.302956 0.9764019 -2.291249 1.700694 -2.752648 2.395755 -2.570197 -0.3147871 -1.866766 diag(mLambda) 0.008453146 0.004908342 0.01116711 0.008812986 0.009714661 0.007520194 0.0091178 0.005607995 0.02092289 0.02143084 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma [1] 100 8 NULL Time difference of 1.688941 secs $var_result $var_result$vmu [,1] [1,] 2.0930326 [2,] 2.3029559 [3,] 0.9764019 [4,] -2.2912491 [5,] 1.7006935 [6,] -2.7526478 [7,] 2.3957545 [8,] -2.5701973 [9,] -0.3147871 [10,] -1.8667662

$var_result$mLambda [,1] [,2] [,3] [,4] [,5] [1,] 0.008453146 -0.006044300 -0.008430260 0.006054035 -0.008441129 [2,] -0.006044300 0.004908342 0.006029262 -0.004913884 0.006036490 [3,] -0.008430260 0.006029262 0.011167113 -0.006484645 0.008418277 [4,] 0.006054035 -0.004913884 -0.006484645 0.008812986 -0.006046209 [5,] -0.008441129 0.006036490 0.008418277 -0.006046209 0.009714661 [6,] 0.006055016 -0.004914795 -0.006039945 0.004920353 -0.005868731 [7,] -0.008446602 0.006040059 0.008423735 -0.006049785 0.008434595 [8,] 0.006048102 -0.004910678 -0.006033051 0.004916226 -0.006040285 [9,] -0.008360106 0.005982404 0.008337483 -0.005992021 0.008348227 [10,] 0.006054223 -0.004911365 -0.006039149 0.004916929 -0.006046393 [,6] [,7] [,8] [,9] [,10] [1,] 0.006055016 -0.008446602 0.006048102 -0.008360106 0.006054223 [2,] -0.004914795 0.006040059 -0.004910678 0.005982404 -0.004911365 [3,] -0.006039945 0.008423735 -0.006033051 0.008337483 -0.006039149 [4,] 0.004920353 -0.006049785 0.004916226 -0.005992021 0.004916929 [5,] -0.005868731 0.008434595 -0.006040285 0.008348227 -0.006046393 [6,] 0.007520194 -0.006050765 0.004917138 -0.005992992 0.004917840 [7,] -0.006050765 0.009117800 -0.005890534 0.008353638 -0.006049971 [8,] 0.004917138 -0.005890534 0.005607995 -0.005986160 0.004913709 [9,] -0.005992992 0.008353638 -0.005986160 0.020922889 -0.011423464 [10,] 0.004917840 -0.006049971 0.004913709 -0.011423464 0.021430839

$var_result$a_rho [1] 84.0033

$var_result$b_rho [1] 17.9967

$var_result$a_sigma numeric(0)

$var_result$b_sigma numeric(0)

$var_result$vlower_bound [1] -4388.317 -4344.679 -4344.660 -4344.659 -4344.659 -4344.659 -4344.659 [8] -4344.659 -4344.659 -4344.659 -4344.659 -4344.659

$vbeta_accuracy [1] 0.9218658 0.9196421

$vu_accuracy [1] 0.9555139 0.9244811 0.9300057 0.9253950 0.9257944 0.9294784 0.9325448 [8] 0.9585505

$rho_accuracy [1] 0.9036515

Linux log:

R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R.

accuracy.R

source("zero_inflated_model.R") source("generate.R") source("mcmc.R") source("rwmh.R")

calculate_accuracy = function(mcmc_samples, dist_fn, param1, param2)

{
mcmc_density = density(mcmc_samples)
mcmc_fn = splinefun(mcmc_density$x, mcmc_density$y) *
integrand <- function(x)
{
return(abs(mcmc_fn(x) - dist_fn(x, param1, param2)))
}
result = integrate(integrand, min(mcmc_density$x), max(mcmc_density$x),
subdivisions = length(mcmc_density$x))
accuracy = 1 - .5 * result$value
return(accuracy)
}

accuracy_plot = function(mcmc_samples, dist_fn, param1, param2)

{
mcmc_density = density(mcmc_samples)
plot(mcmc_density)
curve(dist_fn(x, param1, param2),
from=min(mcmc_density$x), to=max(mcmc_density$x),
add=TRUE, lty=2, col="blue")
}

calculate_accuracies = function(mult, mcmc_samples, var_result, approximation, print_flag=FALSE, plot_flag=FALSE)

{
TODO: Add support for checking the accuracy over multiple dimensions
cubature$adaptIntegrate

*
if (plot_flag) pdf(paste0("accuracyplots", approximation, ".pdf"))
return(var_result)
vbeta accuracy

*
Compare MCMC distribution with variational approximation for each parameter
vnu[i] ~ Normal, dnorm
sigma2_u ~ IG, dgamma(1/x)
rho ~ Beta, dbeta
vr[i] ~ Bernoulli, dbinom
For each parameter of interest,
* estimate density of MCMC
* compare with q distribution using L_1 norm

*
Kernel density estimates of MCMC-estimated posteriors
Use L_1 distance to compare against variational approximations of posteriors

*
vbeta_accuracy = rep(NA, ncol(mult$mX))
for (i in 1:ncol(mult$mX)) {
vbeta_accuracy[i] = calculate_accuracy(mcmc_samples$vbeta[,i], dnorm,
var_result$vmu[i], sqrt(var_result$mLambda[i,i]))
if (print_flag) cat("vbeta[", i, "]", approximation, "accuracy:", vbeta_accuracy[i], "\n") *
param_name = sprintf("vbeta[%d]", i)
if (plot_flag) accuracy_plot(mcmc_samples$vbeta[,i], dnorm,
var_result$vmu[i], sqrt(var_result$mLambda[i,i]))
} *
vu accuracy
FIXME: To check for random slopes accuracy, this section will have
to get more complex.
print(dim(mult$mZ))
print(dim(mcmc_samples$u))
vu_accuracy = rep(NA, ncol(mult$mZ))
B = mult$blocksize
b_idx = 1
for (i in 1:ncol(mult$mZ)) {
m_idx = ceiling(i/B)
vu_accuracy[i] = calculate_accuracy(mcmc_samples$vu[,m_idx,b_idx], dnorm,
var_result$vmu[i+mult$p], sqrt(var_result$mLambda[i+mult$p,i+mult$p]))
if (print_flag) cat("vu[", i, "]", approximation, "accuracy:", vu_accuracy[i], "\n")
if (plot_flag) accuracy_plot(mcmc_samples$vu[,m_idx,b_idx], dnorm,
var_result$vmu[i+mult$p], sqrt(var_result$mLambda[i+mult$p,i+mult$p])) *
b_idx = b_idx + 1
if (b_idx > B)
b_idx=1
} *
sigma2_u accuracy
FIXME - this may be wrong for blocksize != 1?
This is totally wrong for the Inverse Wishart model?
sigma2_u_accuracy = calculate_accuracy(1/mcmc_samples$sigma_u^2, dgamma,
var_result$a_sigma, var_result$b_sigma)
if (print_flag) cat("sigma2_u", approximation, "accuracy:", sigma2_u_accuracy, "\n")
if (plot_flag) accuracy_plot(1/mcmc_samples$sigma_u^2, dgamma,
var_result$a_sigma, var_result$b_sigma)

*
rho accuracy
rho_accuracy = calculate_accuracy(mcmc_samples$rho, dbeta,
var_result$a_rho, var_result$b_rho)
if (print_flag) cat("rho", approximation, "accuracy: ", rho_accuracy, "\n")
if (plot_flag) accuracy_plot(mcmc_samples$rho, dbeta,
var_result$a_rho, var_result$b_rho)
if (plot_flag) dev.off()
return(list(var_result=var_result,
vbeta_accuracy=vbeta_accuracy,
vu_accuracy=vu_accuracy,
sigma2_u_accuracy=sigma2_u_accuracy,
rho_accuracy=rho_accuracy))
}

test_accuracy = function(mult, mcmc_samples, approximation, plot=FALSE)

{
var_result = zero_infl_var(mult, method=approximation, verbose=TRUE)
return(calculate_accuracies(mult, mcmc_samples, var_result, approximation, plot_flag=plot))
}

Calculate accuracy ---- Approximate the L1 norm between the variational approximation and the MCMC approximation

calculate_univariate_accuracy <- function(result_mcmc, var_result)

{
density_mcmc_rho = density(result_mcmc$vrho)
integrand <- function(x)
{
fn = splinefun(density_mcmc_rho$x, density_mcmc_rho$y)
return(abs(fn(x) - dbeta(x, var_result$a_rho, var_result$b_rho)))
}
result = integrate(integrand, min(density_mcmc_rho$x), max(density_mcmc_rho$x), subdivisions = length(density_mcmc_rho$x))
rho_accuracy = 1 - .5 * result$value *
density_mcmc_lambda = density(result_mcmc$vlambda)
integrand <- function(x)
{
fn = splinefun(density_mcmc_lambda$x, density_mcmc_lambda$y)
return(abs(fn(x) - dgamma(x, var_result$a_lambda, var_result$b_lambda)))
}
result = integrate(integrand, min(density_mcmc_lambda$x), max(density_mcmc_lambda$x), subdivisions = length(density_mcmc_lambda$x))
lambda_accuracy = 1 - .5 * result$value
return(list(rho_accuracy=rho_accuracy, lambda_accuracy=lambda_accuracy))
}

test_accuracies = function()

{
Need to be able to compare the solution paths of each approximation

*
Generate data
for (i in 1:100) {
set.seed(i)
mult = generate_test_data(20, 100)
Monte Carlo Markov Chains approximation
mcmc_samples = mcmc_approximation(mult, iterations=1e6)
Save the results, because this takes such a long time to run.
}
save(mult, mcmc_samples, file="accuracy_good.RData")
set.seed(1)
mult = generate_test_data(10, 100)
Monte Carlo Markov Chains approximation
mcmc_samples = mcmc_approximation(mult, iterations=1e6, warmup = 1e4)
Save the results, because this takes such a long time to run.
save(mult, mcmc_samples, file="accuracy.RData")
save(mult, mcmc_samples, file="data/accuracy_int.RData")
load(file="data/accuracy_int_2015_02_17.RData")
mult$spline_dim = 0
load(file="accuracy.RData")
Test all other approximations against it
load(file="accuracy.RData")

*
Test multivariate approximation's accuracy
now = Sys.time()
var1 = test_accuracy(mult, mcmc_samples, "laplacian")
print(Sys.time() - now)
print(image(Matrix(var1$var_result$mLambda)))
print(var1) *
now = Sys.time()
var2 = test_accuracy(mult, mcmc_samples, "gva")
print(Sys.time() - now)
print(image(Matrix(var2$var_result$mLambda)))
print(var2) *
now = Sys.time()
var3 = test_accuracy(mult, mcmc_samples, "gva2", plot=TRUE)
print(Sys.time() - now)
print(image(Matrix(var3$mLambda)))
print(var3) *
Rprof()
now = Sys.time()
var4 = test_accuracy(mult, mcmc_samples, "gva2new")
print(Sys.time() - now)
print(image(Matrix(var4$var_result$mLambda)))
print(var4) *
Rprof(NULL)
print(summaryRprof())
print(image(Matrix(var3_new$mLambda)))

*
now = Sys.time()
var5 = test_accuracy(mult, mcmc_samples, "gva_nr")
print(Sys.time() - now)
print(image(Matrix(var4$result_var$mLambda)))
print(var5) *
save(var1, var2, var3, var4, var5, file="accuracy_results_int.RData")
for (i in 1:100) {
set.seed(i)
mult = generate_test_data(20, 100)
mcmc_samples = mcmc_approximation(mult, iterations=1e4)
#
var1 = test_accuracy(mult, mcmc_samples, "laplacian")
var2 = test_accuracy(mult, mcmc_samples, "gva")
var3 = test_accuracy(mult, mcmc_samples, "gva2")
var4 = test_accuracy(mult, mcmc_samples, "gva_nr")
}

*
} #test_accuracies()

test_accuracies_slope = function()

{
Monte Carlo Markov Chains approximation
seed = 1
set.seed(seed)
mult = generate_slope_test_data()
mcmc_samples = mcmc_approximation(mult, seed=seed, iterations=1e6, warmup=1e5)
save(mult, mcmc_samples, file="data/accuracy_slope_2015_03_03.RData")
load(file="data/accuracy_slope_2015_03_03.RData")

*
now = Sys.time()
var1 = test_accuracy(mult, mcmc_samples, "laplacian", plot=TRUE)
print(Sys.time() - now)
print(var1) *
now = Sys.time()
var2 = test_accuracy(mult, mcmc_samples, "gva", plot=TRUE)
print(Sys.time() - now)
print(image(Matrix(var2$var_result$mLambda)))
print(var2) *
now = Sys.time()
var3 = test_accuracy(mult, mcmc_samples, "gva2", plot=TRUE)
print(image(Matrix(var3$var_result$mLambda)))
print(Sys.time() - now)
print(var3) *
print("gva2new")
now = Sys.time()
Rprof(line.profiling=TRUE)
var4 = test_accuracy(mult, mcmc_samples, "gva2new", plot=TRUE)
Rprof(NULL)
print(summaryRprof(lines = "both"))
print(image(Matrix(var4$var_result$mLambda)))
var4 = test_accuracy(mult, mcmc_samples, "gva2new")
print(Sys.time() - now)
print(var4) *
now = Sys.time()
GVA NR is unstable, and sometimes fails with an error
var5 = test_accuracy(mult, mcmc_samples, "gva_nr")
var5 = test_accuracy(mult, mcmc_samples, "gva_nr", plot=TRUE)
print(Sys.time() - now)
print(image(Matrix(var5$var_result$mLambda)))
print(var5) *
save(var1, var2, var3, var4, var5, file="accuracy_results_slope.RData")

*
print(image(Matrix(var3_new$mLambda)))

*
now = Sys.time()
var4 = test_accuracy(mult, mcmc_samples, "gva_nr")
Sys.time() - now
print(image(Matrix(var4$mLambda)))

*
mult = generate_spline_test_data()
mcmc_samples = mcmc_approximation(mult, iterations=1e4)
save(mult, mcmc_samples, file="accuracy_spline.RData")

*
now = Sys.time()
var1 = test_accuracy(mult, mcmc_samples, "laplacian")
Sys.time() - now
print(image(Matrix(var1$var_result$mLambda)))
var1

*
now = Sys.time()
var2 = test_accuracy(mult, mcmc_samples, "gva")
Sys.time() - now
print(image(Matrix(var2$var_result$mLambda)))
var2

*
now = Sys.time()
var3 = test_accuracy(mult, mcmc_samples, "gva2")
Sys.time() - now
print(image(Matrix(var3$var_result$mLambda)))
var3

*
now = Sys.time()
var4 = test_accuracy(mult, mcmc_samples, "gva2new")
Sys.time() - now
print(image(Matrix(var4$var_result$mLambda)))
var4

*
Rprof()
now = Sys.time()
var5 = test_accuracy(mult, mcmc_samples, "gva_nr")
print(Sys.time() - now)
print(image(Matrix(var5$var_result$mLambda)))
var5
} test_accuracies_slope() N 100 p 2 m 5 blocksize 2 spline_dim 0 calculate_lower_bound: T1 -4.10617e+17 T2 -105.5609 T3 -38.79288 a_sigma b_sigma calculate_lower_bound: T1 -4254.533 T2 -48.13332 T3 -47.93269 Iteration 2 : lower bound -4350.599 difference 4.10617e+17 parameters vmu 1.970248 2.387512 1.101093 -2.379948 1.825163 -2.837605 2.436981 -2.691774 -0.1999721 -1.972204 diag(mLambda) 0.008005845 0.00468602 0.01032378 0.007912696 0.009163322 0.006823373 0.008631102 0.005365371 0.01762683 0.01577659 a_rho 84.0028 b_rho 17.9972 a_sigma b_sigma calculate_lower_bound: T1 -4248.745 T2 -48.09777 T3 -47.92874 Iteration 3 : lower bound -4344.772 difference 5.826987 parameters vmu 2.093207 2.303359 0.9793739 -2.292995 1.701905 -2.752487 2.39622 -2.570463 -0.3002183 -1.871611 diag(mLambda) 0.007165469 0.004302308 0.009434513 0.007534457 0.008300133 0.006445907 0.007727148 0.004922677 0.01642191 0.01512188 a_rho 84.00333 b_rho 17.99667 a_sigma b_sigma calculate_lower_bound: T1 -4248.74 T2 -48.10542 T3 -47.92929 Iteration 4 : lower bound -4344.774 difference -0.002537742 parameters vmu 2.09571 2.301589 0.9766196 -2.291273 1.699333 -2.750677 2.393756 -2.568625 -0.3045866 -1.870705 diag(mLambda) 0.007154014 0.00429596 0.009425739 0.007527698 0.008289927 0.006438819 0.007716174 0.004915965 0.01642994 0.01513047 a_rho 84.00326 b_rho 17.99674 a_sigma b_sigma [1] 100 8 NULL Time difference of 3.77699 secs $var_result $var_result$vmu ,1 2.0957098 x 2.3015888 groups2 0.9766196 x:groups2 -2.2912729 groups3 1.6993328 x:groups3 -2.7506774 groups4 2.3937563 x:groups4 -2.5686248 groups5 -0.3045866 x:groups5 -1.8707046

$var_result$mLambda (Intercept) x groups2 x:groups2 groups3 (Intercept) 0.007154014 -0.005154323 -0.007136626 0.005163922 -0.007144532 x -0.005154323 0.004295960 0.005143068 -0.004301498 0.005148240 groups2 -0.007136626 0.005143068 0.009425739 -0.005162749 0.007127168 x:groups2 0.005163922 -0.004301498 -0.005162749 0.007527698 -0.005157824 groups3 -0.007144532 0.005148240 0.007127168 -0.005157824 0.008289927 x:groups3 0.005162642 -0.004300897 -0.005151365 0.004306450 -0.004903491 groups4 -0.007149408 0.005151363 0.007132031 -0.005160954 0.007139932 x:groups4 0.005156909 -0.004297504 -0.005145647 0.004303047 -0.005150822 groups5 -0.007097451 0.005116721 0.007080207 -0.005126237 0.007088048 x:groups5 0.005150843 -0.004291441 -0.005139592 0.004296980 -0.005144762 x:groups3 groups4 x:groups4 groups5 x:groups5 (Intercept) 0.005162642 -0.007149408 0.005156909 -0.007097451 0.005150843 x -0.004300897 0.005151363 -0.004297504 0.005116721 -0.004291441 groups2 -0.005151365 0.007132031 -0.005145647 0.007080207 -0.005139592 x:groups2 0.004306450 -0.005160954 0.004303047 -0.005126237 0.004296980 groups3 -0.004903491 0.007139932 -0.005150822 0.007088048 -0.005144762 x:groups3 0.006438819 -0.005159675 0.004302446 -0.005124969 0.004296379 groups4 -0.005159675 0.007716174 -0.005056923 0.007092883 -0.005147884 x:groups4 0.004302446 -0.005056923 0.004915965 -0.005119285 0.004292986 groups5 -0.005124969 0.007092883 -0.005119285 0.016429939 -0.009887394 x:groups5 0.004296379 -0.005147884 0.004292986 -0.009887394 0.015130470

$var_result$a_rho [1] 84.00326

$var_result$b_rho [1] 17.99674

$var_result$a_sigma numeric(0)

$var_result$b_sigma numeric(0)

$var_result$vlower_bound [1] -4.106170e+17 -4.350599e+03 -4.344772e+03 -4.344774e+03

$vbeta_accuracy [1] 0.9035163 0.9057978

$vu_accuracy [1] 0.9168333 0.9106049 0.9179837 0.9240022 0.9067588 0.9171403 0.9006352 [8] 0.8843505

$rho_accuracy [1] 0.9294158

N 100 p 2 m 5 blocksize 2 spline_dim 0 calculate_lower_bound: T1 -4290.507 T2 -49.87344 T3 -47.95629 a_sigma b_sigma calculate_lower_bound: T1 -4248.556 T2 -48.19419 T3 -47.92916 Iteration 2 : lower bound -4344.679 difference 43.6573 parameters vmu 2.09409 2.302492 0.974904 -2.291347 1.699599 -2.752365 2.394914 -2.569596 -0.3210042 -1.869959 diag(mLambda) 0.008481679 0.004919913 0.0112074 0.008822642 0.009748994 0.007528633 0.009148985 0.005618053 0.0210578 0.02150873 a_rho 84.00328 b_rho 17.99672 a_sigma b_sigma calculate_lower_bound: T1 -4248.571 T2 -48.15913 T3 -47.92901 Iteration 3 : lower bound -4344.66 difference 0.01940456 parameters vmu 2.093146 2.302883 0.9762685 -2.291192 1.700576 -2.752579 2.395646 -2.570119 -0.3150815 -1.866813 diag(mLambda) 0.008453228 0.004908299 0.01116755 0.008812885 0.009714915 0.00752008 0.009117964 0.005607909 0.02092631 0.02143282 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15807 T3 -47.92899 Iteration 4 : lower bound -4344.659 difference 0.0005529186 parameters vmu 2.093039 2.302952 0.976395 -2.291246 1.700687 -2.752644 2.395749 -2.570193 -0.3148003 -1.866767 diag(mLambda) 0.008453131 0.00490833 0.01116711 0.008812966 0.009714652 0.007520178 0.009117788 0.005607981 0.02092304 0.02143098 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 5 : lower bound -4344.659 difference 2.113186e-05 parameters vmu 2.093033 2.302956 0.9764016 -2.291249 1.700693 -2.752648 2.395754 -2.570197 -0.3147875 -1.866766 diag(mLambda) 0.008453145 0.004908341 0.01116711 0.008812988 0.009714659 0.007520194 0.009117799 0.005607995 0.02092289 0.02143083 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 6 : lower bound -4344.659 difference 5.468974e-07 parameters vmu 2.093032 2.302956 0.9764023 -2.291249 1.700694 -2.752648 2.395755 -2.570198 -0.3147866 -1.866767 diag(mLambda) 0.008453146 0.004908345 0.01116712 0.008812998 0.009714662 0.007520199 0.009117801 0.005607998 0.02092289 0.02143081 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 7 : lower bound -4344.659 difference 3.403011e-07 parameters vmu 2.093033 2.302956 0.9764017 -2.291249 1.700693 -2.752648 2.395754 -2.570197 -0.3147873 -1.866766 diag(mLambda) 0.008453147 0.004908342 0.01116711 0.008812986 0.009714661 0.007520193 0.009117799 0.005607994 0.02092288 0.02143082 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma calculate_lower_bound: T1 -4248.572 T2 -48.15803 T3 -47.92899 Iteration 8 : lower bound -4344.659 difference -1.964881e-07 parameters vmu 2.093032 2.302956 0.9764023 -2.291249 1.700694 -2.752648 2.395755 -2.570198 -0.3147865 -1.866767 diag(mLambda) 0.008453147 0.004908342 0.01116711 0.00881299 0.00971466 0.007520198 0.009117801 0.005607997 0.02092289 0.02143083 a_rho 84.0033 b_rho 17.9967 a_sigma b_sigma [1] 100 8 NULL Time difference of 3.73424 secs $var_result $var_result$vmu [1] 2.0930321 2.3029562 0.9764023 -2.2912495 1.7006940 -2.7526481 [7] 2.3957550 -2.5701976 -0.3147865 -1.8667666

$var_result$mLambda [,1] [,2] [,3] [,4] [,5] [1,] 0.008453147 -0.006044300 -0.008430261 0.006054039 -0.008441128 [2,] -0.006044300 0.004908342 0.006029261 -0.004913886 0.006036489 [3,] -0.008430261 0.006029261 0.011167115 -0.006484648 0.008418278 [4,] 0.006054039 -0.004913886 -0.006484648 0.008812990 -0.006046212 [5,] -0.008441128 0.006036489 0.008418278 -0.006046212 0.009714660 [6,] 0.006055017 -0.004914797 -0.006039947 0.004920356 -0.005868732 [7,] -0.008446602 0.006040059 0.008423736 -0.006049789 0.008434593 [8,] 0.006048102 -0.004910679 -0.006033052 0.004916228 -0.006040285 [9,] -0.008360108 0.005982405 0.008337484 -0.005992025 0.008348227 [10,] 0.006054225 -0.004911368 -0.006039150 0.004916931 -0.006046394 [,6] [,7] [,8] [,9] [,10] [1,] 0.006055017 -0.008446602 0.006048102 -0.008360108 0.006054225 [2,] -0.004914797 0.006040059 -0.004910679 0.005982405 -0.004911368 [3,] -0.006039947 0.008423736 -0.006033052 0.008337484 -0.006039150 [4,] 0.004920356 -0.006049789 0.004916228 -0.005992025 0.004916931 [5,] -0.005868732 0.008434593 -0.006040285 0.008348227 -0.006046394 [6,] 0.007520198 -0.006050766 0.004917140 -0.005992994 0.004917842 [7,] -0.006050766 0.009117801 -0.005890535 0.008353638 -0.006049972 [8,] 0.004917140 -0.005890535 0.005607997 -0.005986161 0.004913711 [9,] -0.005992994 0.008353638 -0.005986161 0.020922894 -0.011423447 [10,] 0.004917842 -0.006049972 0.004913711 -0.011423447 0.021430829

$var_result$a_rho [1] 84.0033

$var_result$b_rho [1] 17.9967

$var_result$a_sigma numeric(0)

$var_result$b_sigma numeric(0)

$var_result$vlower_bound [1] -4388.336 -4344.679 -4344.660 -4344.659 -4344.659 -4344.659 -4344.659 [8] -4344.659

$vbeta_accuracy [1] 0.9308112 0.9251300

$vu_accuracy [1] 0.9332386 0.9255522 0.9373222 0.9434991 0.9286187 0.9388787 0.9601703 [8] 0.9324045

$rho_accuracy [1] 0.9294196

N 100 p 2 m 5 blocksize 2 spline_dim 0 Error in optim(par = vmu, fn = f.GVA_new, gr = vg.GVA_new, method = method, : non-finite value supplied by optim Calls: test_accuracies_slope ... zero_infl_var -> zero_infl_var.multivariate -> fit.GVA_new -> optim Execution halted

Reply to this email directly or view it on GitHubhttps://github.com/certifiedwaif/phd/issues/8.

certifiedwaif commented 9 years ago

Indeed.

Hey, I bet this isn't really helping ...

Browse[2]> mR[1:2, 1:2]
            (Intercept)       x
(Intercept)  1056490771       0
x            1685014370 5848370
Browse[2]> log(mR[1:2, 1:2])
            (Intercept)        x
(Intercept)    20.77822     -Inf
x              21.24504 15.58167

So my guess would be it's the inverse Cholesky parameterisation combined with the fact that we're only taking logs and exponentials of the diagonal of mLambda. If we initialised differently or took logs of the entire lower triangular matrix, we could probably side-step this.

But yes, it really unsettles me when the same code with the same random number seed behaves differently on different computers. R and all the other tools that we're using are supposed to be portable, and minor revisions shouldn't break things.

certifiedwaif commented 9 years ago

Think I've caught it. We were initialising the parameters for the GVA2 fit using the Laplacian approximation first, and that was causing the crash. So I'm initialising the model fit parameters normally now.

https://github.com/certifiedwaif/phd/commit/656c83537449492ae9338229e639ab9abef160f0

certifiedwaif commented 9 years ago

Resolved by frequent R version upgrades.

certifiedwaif / phd

All of the algorithms run on R 3.1.0 on my Mac, but fail on R 3.1.2 on the Linux machines at Uni. #8

accuracy.R

TODO: Add support for checking the accuracy over multiple dimensions

cubature$adaptIntegrate

return(var_result)

vbeta accuracy

Compare MCMC distribution with variational approximation for each parameter

vnu[i] ~ Normal, dnorm

sigma2_u ~ IG, dgamma(1/x)

rho ~ Beta, dbeta

vr[i] ~ Bernoulli, dbinom

For each parameter of interest,

* estimate density of MCMC

* compare with q distribution using L_1 norm

Kernel density estimates of MCMC-estimated posteriors

Use L_1 distance to compare against variational approximations of posteriors

vu accuracy

FIXME: To check for random slopes accuracy, this section will have

to get more complex.

TODO: The B - (i %% B) expression only works for B=2. Rewrite this.

sigma2_u accuracy

FIXME - this may be wrong for blocksize != 1?

This is totally wrong for the Inverse Wishart model?

sigma2_u_accuracy = calculate_accuracy(1/mcmc_samples$sigma_u^2, dgamma,

var_result$a_sigma, var_result$b_sigma)

if (print_flag) cat("sigma2_u", approximation, "accuracy:", sigma2_u_accuracy, "\n")

if (plot_flag) accuracy_plot(1/mcmc_samples$sigma_u^2, dgamma,

var_result$a_sigma, var_result$b_sigma)

rho accuracy

sigma2_u_accuracy=sigma2_u_accuracy,

Calculate accuracy ----

Approximate the L1 norm between the variational approximation and

the MCMC approximation

Need to be able to compare the solution paths of each approximation

Generate data

for (i in 1:100) {

set.seed(i)

mult = generate_test_data(20, 100)

Monte Carlo Markov Chains approximation

mcmc_samples = mcmc_approximation(mult, iterations=1e6)

Save the results, because this takes such a long time to run.

}

save(mult, mcmc_samples, file="accuracy_good.RData")

set.seed(1)

mult = generate_test_data(10, 100)

Monte Carlo Markov Chains approximation

mcmc_samples = mcmc_approximation(mult, iterations=1e6, warmup = 1e4)

Save the results, because this takes such a long time to run.

save(mult, mcmc_samples, file="accuracy.RData")

save(mult, mcmc_samples, file="data/accuracy_int.RData")

mult$spline_dim = 0

load(file="accuracy.RData")

Test all other approximations against it

load(file="accuracy.RData")

Test multivariate approximation's accuracy

print(image(Matrix(var1$var_result$mLambda)))

print(image(Matrix(var2$var_result$mLambda)))

print(image(Matrix(var3$mLambda)))

Rprof()

print(image(Matrix(var4$var_result$mLambda)))

Rprof(NULL)

print(summaryRprof())

print(image(Matrix(var3_new$mLambda)))

print(image(Matrix(var4$result_var$mLambda)))

save(var1, var2, var3, var4, var5, file="accuracy_results_int.RData")

for (i in 1:100) {

set.seed(i)

mult = generate_test_data(20, 100)

mcmc_samples = mcmc_approximation(mult, iterations=1e4)

var1 = test_accuracy(mult, mcmc_samples, "laplacian")

var2 = test_accuracy(mult, mcmc_samples, "gva")

var3 = test_accuracy(mult, mcmc_samples, "gva2")

var4 = test_accuracy(mult, mcmc_samples, "gva_nr")

}

test_accuracies()

Monte Carlo Markov Chains approximation

set.seed(1)

mult = generate_slope_test_data()

mcmc_samples = mcmc_approximation(mult, iterations=1e6, warmup=1e4, mc.cores = 1)