#################################################################################################################### #### Simulates from the posterior distribution of the parameters $(\mu,\sigma,\gamma)$ of the skew-Exponential #### power distribution (with fixed shape parameter) and the skew-Student's t distribution #### (with fixed degrees of freedom) using the proposal prior in Rubio and Steel (2011) such that AG ~ Beta(a,b) #### with the parameterization in Mudholkar and Hutson (2000). #### See http://www2.warwick.ac.uk/fac/sci/statistics/crism/research/2011/paper11-13 #### The posterior distributions are sampled using the command 'metrop' from the library 'mcmc'. #### Author F. Javier Rubio ##################################################################################################################### rm(list=ls()) # Required packages library(mcmc) library(normalp) #-------------------------------------------------------------------------------------------------------- # data : Vector of observations. # model : 1 for Exponential Power distribution; 2 for Student's t distribution. # a : Shape parameter in the proposal prior AG ~ Beta(a,b). # b : Shape parameter in the proposal prior AG ~ Beta(a,b). # ap : Additional parameter, if model = 1 this represents the shape parameter of the exponential # power distribution; if model = 2 this represents the degrees of freedom of the Student's t # distribution. # sc : Controls the proposal step size in the Metropolis sampling algorithm. # burnin : Burn-in period. # thin : Thinning period. # NMH : Length of the Markov Chain. # ini : Starting point for the MCMC sampler. #-------------------------------------------------------------------------------------------------------- post = function(data,model,a,b,ap,sc,burnin,thin,NMH,ini){ ns = length(data) if(model==1){ lp = function(par){ if(par[3]>-1 & par[3]<1 & par[2]>0) return( -(ns+1)*log(par[2]) + sum(log( dnormp((data-par[1])/(par[2]*( 1 - par[3]*sign(data-par[1]) )), p=ap) )) + (a-1)*log(1-par[3]) + (b-1)*log(1+par[3]) ) else return(-Inf) } out <- metrop(lp, scale = sc, initial = ini, nbatch = NMH) ind = seq(burnin,NMH,thin) sim = t(rbind( out$batch[ , 1][ind], out$batch[ , 2][ind],-out$batch[ , 3][ind])) colnames(sim) <- c("mu","sigma","gamma") acc = out$accept cat('Acceptance rate = ') cat(acc, fill = TRUE) return(sim) } if(model==2){ lp = function(par){ if(par[3]>-1 & par[3]<1 & par[2]>0) return( -(ns+1)*log(par[2]) + sum(log( dt((data-par[1])/(par[2]*( 1 - par[3]*sign(data-par[1]) )), df=ap) )) + (a-1)*log(1-par[3]) + (b-1)*log(1+par[3]) ) else return(-Inf) } out <- metrop(lp, scale = sc, initial = ini, nbatch = NMH) ind = seq(burnin,NMH,thin) sim = t(rbind( out$batch[ , 1][ind], out$batch[ , 2][ind],-out$batch[ , 3][ind])) colnames(sim) <- c("mu","sigma","gamma") acc = out$accept cat('Acceptance rate = ') cat(acc, fill = TRUE) return(sim) } } ###################################################################################################### # Example with simulated data from a Normal(0,1). AG ~ Uniform(-1,1) ###################################################################################################### p = post(rnorm(100),1,1,1,2,0.1,1000,10,26000,c(0,1,0)) hist(p[,1],probability=TRUE,xlab=~mu,main=~mu) hist(p[,2],probability=TRUE,xlab=~sigma,main=~sigma) hist(p[,3],probability=TRUE,xlab=~gamma,main=~gamma) ###################################################################################################### # Example with real data from Mudholkar and Hutson (2000). AG ~ Uniform(-1,1) ###################################################################################################### dat = 10*c(0.9,0.9,0.7,0.6,0.6,0.6,0.5,0.2, 1.9,1.9,1.7,1.6,1.6,1.3,1.1,1.0,1.0,1.0, 2.9,2.9,2.9,2.8,2.8,2.7,2.7,2.6,2.6,2.6,2.5,2.5,2.4,2.4,2.4,2.2,2.2,2.2,2.1,2.1,2.0,2.0, 3.9,3.9,3.8,3.7,3.6,3.6,3.6,3.5,3.5,3.5,3.5,3.4,3.4,3.2,3.2,3.1,3.1,3.0, 4.9,4.9,4.9,4.9,4.9,4.8,4.8,4.7,4.6,4.4,4.4,4.4,4.3,4.3,4.3,4.3,4.3,4.2,4.1,4.1,4.1,4.0, 5.9,5.9,5.7,5.7,5.7,5.6,5.6,5.6,5.6,5.6,5.6,5.5,5.5,5.4,5.4,5.3,5.2,5.2,5.2,5.1,5.1,5.0,5.0, 6.9,6.8,6.8,6.7,6.7,6.7,6.6,6.6,6.6,6.6,6.5,6.5,6.4,6.4,6.1,6.1,6.0,6.0, 7.9,7.8,7.8,7.8,7.7,7.6,7.5,7.5,7.5,7.4,7.3,7.3,7.2,7.1,7.1,7.1,7.0,7.0,7.0,7.0, 8.9,8.7,8.6,8.5,8.3,8.3,8.3,8.2,8.2,8.2,8.2,8.1, 9.9,9.7,9.7,9.6,9.5,9.5,9.4,9.4,9.3,9.3,9.2,9.1,9.0,9.0,9.0, 10.9,10.8,10.6,10.5,10.4,10.4,10.3,10.3,10.2,10.2,10.1,10.1,10.0, 11.9,11.6,11.6,11.4,11.3,11.3,11.2,11.1,11.1,11.0, 12.6,12.5,12.4,12.4,12.4,12.2,12.1,12.1, 13.8,13.7,13.4,13.3,13.0, 14.0,14.0, 15.7,15.6,15.6, 16.5,16.2, 17.9,17.2, 18.5, 19.9,19.7,19.3,19.3,19.0) hist(dat,probability=TRUE,xlab="Height", main="Volcanoes' data") p = post(dat,1,1,1,2,0.1,50000,100,550000,c(30,40,0.65)) hist(p[,1],probability=TRUE,xlab=~mu,main=~mu) hist(p[,2],probability=TRUE,xlab=~sigma,main=~sigma) hist(p[,3],probability=TRUE,xlab=~gamma,main=~gamma)