SSVS of George and McCulloch after Francesca's SSVSsim.b % data generation (matlab ) % the model includes 3rd, 5th and 7th variable. % beta_3 = -2.2; beta_5=4; beta_7 = 2.46 % tau = 1; clear all; rand('seed',1) randn('seed',1) x=floor(10 * rand(10,7) ); Y = -2.2 * x(:,3) + 4.0 * x(:,5) + 2.46 * x(:,7) + randn(10,1); #Model Begin model { for ( i in 1:nn ) { #nn is sample size Y[i] ~ dnorm(mu[i],tau) mu[i] <- beta[1]*x[i,1]+beta[2]*x[i,2]+beta[3]*x[i,3]+ beta[4]*x[i,4]+beta[5]*x[i,5]+beta[6]*x[i,6] +beta[7]*x[i,7] } #------------ for (j in 1:7) { #j is the index of variables gamma[j] <- g[j,k] #k is model number beta[j] ~ dnorm(0,tau.beta[j]) tau.beta[j] <- equals(gamma[j],1)*tau.1+equals(gamma[j],0)*tau.2 } #Priors: tau ~ dgamma(0.001,0.001) # discrete prior on set of all N models for (n in 1:N) { prior[n] <- 1/N} k ~ dcat(prior[]); # alternating precision parameters in tau.beta tau.1 <- 0.1 tau.2 <- 10 # Specification of models by gamma. The decision maker believes that variable #1 is NOT # in the model, however the rest 6 variables give 2^6=64 plausible models that are apriori # assumed equiprobable g[1,1]<-0; g[2,1]<-0; g[3,1]<-0; g[4,1]<-0; g[5,1]<-0; g[6,1]<-0; g[7,1]<-0; #0000000 g[1,2]<-0; g[2,2]<-0; g[3,2]<-0; g[4,2]<-0; g[5,2]<-0; g[6,2]<-0; g[7,2]<-1; #0000001 g[1,3]<-0; g[2,3]<-0; g[3,3]<-0; g[4,3]<-0; g[5,3]<-0; g[6,3]<-1; g[7,3]<-0; #0000010 g[1,4]<-0; g[2,4]<-0; g[3,4]<-0; g[4,4]<-0; g[5,4]<-0; g[6,4]<-1; g[7,4]<-1; #0000011 g[1,5]<-0; g[2,5]<-0; g[3,5]<-0; g[4,5]<-0; g[5,5]<-1; g[6,5]<-0; g[7,5]<-0; #0000100 g[1,6]<-0; g[2,6]<-0; g[3,6]<-0; g[4,6]<-0; g[5,6]<-1; g[6,6]<-0; g[7,6]<-1; #0000101 g[1,7]<-0; g[2,7]<-0; g[3,7]<-0; g[4,7]<-0; g[5,7]<-1; g[6,7]<-1; g[7,7]<-0; #0000110 g[1,8]<-0; g[2,8]<-0; g[3,8]<-0; g[4,8]<-0; g[5,8]<-1; g[6,8]<-1; g[7,8]<-1; #0000111 g[1,9]<-0; g[2,9]<-0; g[3,9]<-0; g[4,9]<-1; g[5,9]<-0; g[6,9]<-0; g[7,9]<-0; #0001000 g[1,10]<-0; g[2,10]<-0; g[3,10]<-0; g[4,10]<-1; g[5,10]<-0; g[6,10]<-0; g[7,10]<-1; #0001001 g[1,11]<-0; g[2,11]<-0; g[3,11]<-0; g[4,11]<-1; g[5,11]<-0; g[6,11]<-1; g[7,11]<-0; #0001010 g[1,12]<-0; g[2,12]<-0; g[3,12]<-0; g[4,12]<-1; g[5,12]<-0; g[6,12]<-1; g[7,12]<-1; #0001011 g[1,13]<-0; g[2,13]<-0; g[3,13]<-0; g[4,13]<-1; g[5,13]<-1; g[6,13]<-0; g[7,13]<-0; #0001100 g[1,14]<-0; g[2,14]<-0; g[3,14]<-0; g[4,14]<-1; g[5,14]<-1; g[6,14]<-0; g[7,14]<-1; #0001101 g[1,15]<-0; g[2,15]<-0; g[3,15]<-0; g[4,15]<-1; g[5,15]<-1; g[6,15]<-1; g[7,15]<-0; #0001110 g[1,16]<-0; g[2,16]<-0; g[3,16]<-0; g[4,16]<-1; g[5,16]<-1; g[6,16]<-1; g[7,16]<-1; #0001111 g[1,17]<-0; g[2,17]<-0; g[3,17]<-1; g[4,17]<-0; g[5,17]<-0; g[6,17]<-0; g[7,17]<-0; #0010000 g[1,18]<-0; g[2,18]<-0; g[3,18]<-1; g[4,18]<-0; g[5,18]<-0; g[6,18]<-0; g[7,18]<-1; #0010001 g[1,19]<-0; g[2,19]<-0; g[3,19]<-1; g[4,19]<-0; g[5,19]<-0; g[6,19]<-1; g[7,19]<-0; #0010010 g[1,20]<-0; g[2,20]<-0; g[3,20]<-1; g[4,20]<-0; g[5,20]<-0; g[6,20]<-1; g[7,20]<-1; #0010011 g[1,21]<-0; g[2,21]<-0; g[3,21]<-1; g[4,21]<-0; g[5,21]<-1; g[6,21]<-0; g[7,21]<-0; #0010100 g[1,22]<-0; g[2,22]<-0; g[3,22]<-1; g[4,22]<-0; g[5,22]<-1; g[6,22]<-0; g[7,22]<-1; #0010101***** g[1,23]<-0; g[2,23]<-0; g[3,23]<-1; g[4,23]<-0; g[5,23]<-1; g[6,23]<-1; g[7,23]<-0; #0010110 g[1,24]<-0; g[2,24]<-0; g[3,24]<-1; g[4,24]<-0; g[5,24]<-1; g[6,24]<-1; g[7,24]<-1; #0010111 g[1,25]<-0; g[2,25]<-0; g[3,25]<-1; g[4,25]<-1; g[5,25]<-0; g[6,25]<-0; g[7,25]<-0; #0011000 g[1,26]<-0; g[2,26]<-0; g[3,26]<-1; g[4,26]<-1; g[5,26]<-0; g[6,26]<-0; g[7,26]<-1; #0011001 g[1,27]<-0; g[2,27]<-0; g[3,27]<-1; g[4,27]<-1; g[5,27]<-0; g[6,27]<-1; g[7,27]<-0; #0011010 g[1,28]<-0; g[2,28]<-0; g[3,28]<-1; g[4,28]<-1; g[5,28]<-0; g[6,28]<-1; g[7,28]<-1; #0011011 g[1,29]<-0; g[2,29]<-0; g[3,29]<-1; g[4,29]<-1; g[5,29]<-1; g[6,29]<-0; g[7,29]<-0; #0011100 g[1,30]<-0; g[2,30]<-0; g[3,30]<-1; g[4,30]<-1; g[5,30]<-1; g[6,30]<-0; g[7,30]<-1; #0011101 g[1,31]<-0; g[2,31]<-0; g[3,31]<-1; g[4,31]<-1; g[5,31]<-1; g[6,31]<-1; g[7,31]<-0; #0011110 g[1,32]<-0; g[2,32]<-0; g[3,32]<-1; g[4,32]<-1; g[5,32]<-1; g[6,32]<-1; g[7,32]<-1; #0011111 g[1,33]<-0; g[2,33]<-1; g[3,33]<-0; g[4,33]<-0; g[5,33]<-0; g[6,33]<-0; g[7,33]<-0; #0100000 g[1,34]<-0; g[2,34]<-1; g[3,34]<-0; g[4,34]<-0; g[5,34]<-0; g[6,34]<-0; g[7,34]<-1; #0100001 g[1,35]<-0; g[2,35]<-1; g[3,35]<-0; g[4,35]<-0; g[5,35]<-0; g[6,35]<-1; g[7,35]<-0; #0100010 g[1,36]<-0; g[2,36]<-1; g[3,36]<-0; g[4,36]<-0; g[5,36]<-0; g[6,36]<-1; g[7,36]<-1; #0100011 g[1,37]<-0; g[2,37]<-1; g[3,37]<-0; g[4,37]<-0; g[5,37]<-1; g[6,37]<-0; g[7,37]<-0; #0100100 g[1,38]<-0; g[2,38]<-1; g[3,38]<-0; g[4,38]<-0; g[5,38]<-1; g[6,38]<-0; g[7,38]<-1; #0100101 g[1,39]<-0; g[2,39]<-1; g[3,39]<-0; g[4,39]<-0; g[5,39]<-1; g[6,39]<-1; g[7,39]<-0; #0100110 g[1,40]<-0; g[2,40]<-1; g[3,40]<-0; g[4,40]<-0; g[5,40]<-1; g[6,40]<-1; g[7,40]<-1; #0100111 g[1,41]<-0; g[2,41]<-1; g[3,41]<-0; g[4,41]<-1; g[5,41]<-0; g[6,41]<-0; g[7,41]<-0; #0101000 g[1,42]<-0; g[2,42]<-1; g[3,42]<-0; g[4,42]<-1; g[5,42]<-0; g[6,42]<-0; g[7,42]<-1; #0101001 g[1,43]<-0; g[2,43]<-1; g[3,43]<-0; g[4,43]<-1; g[5,43]<-0; g[6,43]<-1; g[7,43]<-0; #0101010 g[1,44]<-0; g[2,44]<-1; g[3,44]<-0; g[4,44]<-1; g[5,44]<-0; g[6,44]<-1; g[7,44]<-1; #0101011 g[1,45]<-0; g[2,45]<-1; g[3,45]<-0; g[4,45]<-1; g[5,45]<-1; g[6,45]<-0; g[7,45]<-0; #0101100 g[1,46]<-0; g[2,46]<-1; g[3,46]<-0; g[4,46]<-1; g[5,46]<-1; g[6,46]<-0; g[7,46]<-1; #0101101 g[1,47]<-0; g[2,47]<-1; g[3,47]<-0; g[4,47]<-1; g[5,47]<-1; g[6,47]<-1; g[7,47]<-0; #0101110 g[1,48]<-0; g[2,48]<-1; g[3,48]<-0; g[4,48]<-1; g[5,48]<-1; g[6,48]<-1; g[7,48]<-1; #0101111 g[1,49]<-0; g[2,49]<-1; g[3,49]<-1; g[4,49]<-0; g[5,49]<-0; g[6,49]<-0; g[7,49]<-0; #0110000 g[1,50]<-0; g[2,50]<-1; g[3,50]<-1; g[4,50]<-0; g[5,50]<-0; g[6,50]<-0; g[7,50]<-1; #0110001 g[1,51]<-0; g[2,51]<-1; g[3,51]<-1; g[4,51]<-0; g[5,51]<-0; g[6,51]<-1; g[7,51]<-0; #0110010 g[1,52]<-0; g[2,52]<-1; g[3,52]<-1; g[4,52]<-0; g[5,52]<-0; g[6,52]<-1; g[7,52]<-1; #0110011 g[1,53]<-0; g[2,53]<-1; g[3,53]<-1; g[4,53]<-0; g[5,53]<-1; g[6,53]<-0; g[7,53]<-0; #0110100 g[1,54]<-0; g[2,54]<-1; g[3,54]<-1; g[4,54]<-0; g[5,54]<-1; g[6,54]<-0; g[7,54]<-1; #0110101 g[1,55]<-0; g[2,55]<-1; g[3,55]<-1; g[4,55]<-0; g[5,55]<-1; g[6,55]<-1; g[7,55]<-0; #0110110 g[1,56]<-0; g[2,56]<-1; g[3,56]<-1; g[4,56]<-0; g[5,56]<-1; g[6,56]<-1; g[7,56]<-1; #0110111 g[1,57]<-0; g[2,57]<-1; g[3,57]<-1; g[4,57]<-1; g[5,57]<-0; g[6,57]<-0; g[7,57]<-0; #0111000 g[1,58]<-0; g[2,58]<-1; g[3,58]<-1; g[4,58]<-1; g[5,58]<-0; g[6,58]<-0; g[7,58]<-1; #0111001 g[1,59]<-0; g[2,59]<-1; g[3,59]<-1; g[4,59]<-1; g[5,59]<-0; g[6,59]<-1; g[7,59]<-0; #0111010 g[1,60]<-0; g[2,60]<-1; g[3,60]<-1; g[4,60]<-1; g[5,60]<-0; g[6,60]<-1; g[7,60]<-1; #0111011 g[1,61]<-0; g[2,61]<-1; g[3,61]<-1; g[4,61]<-1; g[5,61]<-1; g[6,61]<-0; g[7,61]<-0; #0111100 g[1,62]<-0; g[2,62]<-1; g[3,62]<-1; g[4,62]<-1; g[5,62]<-1; g[6,62]<-0; g[7,62]<-1; #0111101 g[1,63]<-0; g[2,63]<-1; g[3,63]<-1; g[4,63]<-1; g[5,63]<-1; g[6,63]<-1; g[7,63]<-0; #0111110 g[1,64]<-0; g[2,64]<-1; g[3,64]<-1; g[4,64]<-1; g[5,64]<-1; g[6,64]<-1; g[7,64]<-1; #0111111 } #Model End #Data Begin list( nn=10, #sample size N=64, #number of apriori plausible models Y=c( 37.7594, 18.9744, 22.4716, 38.3637, 40.3193, 6.5916, 33.8860, 1.5786, 7.9400, 23.3282), x=structure(.Data=c( 5, 1, 3, 7, 9, 0, 3, 4, 4, 3, 4, 4, 6, 4, 3, 1, 1, 1, 2, 7, 7, 0, 7, 0, 6, 9, 2, 1, 4 , 6, 3, 8, 9, 7, 5, 7, 4 , 8, 1, 5, 2, 2, 1, 4 , 0, 6, 8, 4, 0, 0, 5, 2, 1, 1, 8, 2, 3, 5 , 9, 6, 3, 8, 6, 0, 5, 2, 8, 6, 1, 2), .Dim=c(10,7))) #Data End #Inits Begin list( tau=0.4, k=4) #Inits End