Chain Temperatures

[FrankD]

@BennMacdonald : I came across an issue with how the chain temperatures are calculated while preparing the package for CRAN. I changed the chainNum in our example in the agm() function documentation from 5 to 4 (because it was taking more than 10 seconds to run on one of the test servers, which is frowned upon by CRAN, and I thought fewer chains might reduce the runtime). This threw an error, that I tracked down to the chain temperatures. You calculate them in this code snippet in doMCMC():

powers <- 1:chainNum
temp.exponent <- powers[options$temps]
temperatures <- c()
for(CHAIN in 1:chainNum)
{
  temperatures[CHAIN] <- (CHAIN/chainNum)^temp.exponent
}

I think this is fine (although why is the lowest chain temperature not 0, i.e. draws from the prior?), but there is an issue with options$temps, which is set to 5 in this line in agm():

temperatureExponentChains <- 5 # Friel and Pettit

The result is that we are trying to access powers[5], which you can see from the code above only has 4 values when chainNum=4. I suspect that maybe you wanted to set temp.exponent=5; is this it? I don't mind going with that as the default, although I've not read the Friel and Pettit paper to know why it should be used.

[BennMacdonald]

HI Frank,

I will double check the code and get back to you, but yes, it should be that the temp.exponent = 5. This is actually from Calderhead’s MSc thesis (they tried a range of values and found this to be empirically the best), so I will update the comment in the code.

To answer the other question (why we don’t have 0 included I.e. the prior), this is because if 0 is included it can result in a singularity when calculating the log marginal likelihood using thermodynamic integration. To avoid this issue, the value goes close to 0, but never exactly 0.

I will examine the code to see exactly what I was thinking for temp.exponent and get back to you soon.

Thanks for bringing this to my attention.

[BennMacdonald]

Hi Frank,

I checked the code and it would appear that this code

powers <- 1:chainNum

temp.exponent <- powers[options$temps]

temperatures <- c()

for(CHAIN in 1:chainNum)

{

  temperatures[CHAIN] <- (CHAIN/chainNum)^temp.exponent

}

is a remnant of something I was testing. At some point, I decided to test different temperature exponents and the code

temperatureExponentChains <- 5 # Friel and Pettit

referred to the position of the vector “powers” (the different power values I was testing). After I finished, I kept the code in the above fashion, as for my purposes, I never went below 10 chains and so if I wanted to implement Calderhead’s empirical finding that power 5 was superior, I could just set temperatureExponentChains to be 5.

The code should be changed to the following:

temp.exponent <- options$temps

temperatures <- c()

for(CHAIN in 1:chainNum)

{

  temperatures[CHAIN] <- (CHAIN/chainNum)^temp.exponent

}

temperatureExponentChains <- 5 # Calderhead MSc thesis 2008

Do you want me to make the changes or do you want to? If you have a more recent update of the code than what is on GitHub, it might be better for you to update it. I don’t mind either, just let me know what you prefer. If you prefer I do it, please confirm the files on GitHub are the most up to date.

Thanks again for finding this.

[FrankD]

Great, that makes sense! I've already put something like that in as a temporary fix without committing it, so I'll just commit it tomorrow.

[FrankD]

Just pushed the fix, should be fine now.