# change directory in R to the place where the following files are stored: # nb10-data.txt # nb10-model-1-mu.txt # nb10-model-1-sigma.txt # nb10-model-2-mu.txt # nb10-model-2-sigma.txt # nb10-model-2-nu.txt y <- dget( "nb10-data.txt" ) y <- sort( y$y ) mu.G <- matrix( scan( "nb10-model-1-mu.txt" ), 100000, 2, byrow = T )[ , 2 ] sigma.G <- matrix( scan( "nb10-model-1-sigma.txt" ), 100000, 2, byrow = T )[ , 2 ] mu.t <- matrix( scan( "nb10-model-2-mu.txt" ), 100000, 2, byrow = T )[ , 2 ] sigma.t <- matrix( scan( "nb10-model-2-sigma.txt" ), 100000, 2, byrow = T )[ , 2 ] nu.t <- matrix( scan( "nb10-model-2-nu.txt" ), 100000, 2, byrow = T )[ , 2 ] dt.s <- function( y, mu, sigma, nu ) { exp( lgamma( ( nu + 1 ) / 2 ) - ( ( nu + 1 ) / 2 ) * log( 1 + ( y - mu )^2 / ( nu * sigma^2 ) ) - lgamma( nu / 2 ) - log( nu * pi ) / 2 - log( sigma ) ) } LS.contributions <- matrix( 0, 100, 2 ) for ( j in 1:100 ) { LS.contributions[ j, 1 ] <- log( mean( dt.s( y[ j ], mu.t, sigma.t, nu.t ) ) ) LS.contributions[ j, 2 ] <- log( mean( dnorm( y[ j ], mu.G, sigma.G ) ) ) } cbind( y, LS.contributions, 0 + LS.contributions[ , 1 ] > LS.contributions[ , 2 ] ) sum( LS.contributions[ , 1 ] > LS.contributions[ , 2 ] ) / length( y ) LS.t <- mean( LS.contributions[ , 1 ] ) LS.G <- mean( LS.contributions[ , 2 ] ) c( LS.t, LS.G ) plot( y, LS.contributions[ , 1 ], ylim = c( min( LS.contributions ), max( LS.contributions ) ), ylab = 'Log Score Contributions' ) lines( y, LS.contributions[ , 1 ], lty = 1 ) points( y, LS.contributions[ , 2 ], pch = 2 ) lines( y, LS.contributions[ , 2 ], lty = 2 ) legend( 397.5, -10, c( "t", "Gaussian" ), pch = c( 1, 2 ) )