argiopetech/base

SampleMass Metropolis-Hastings sampler

Closed this issue · 1 comments

The current brute force numerical integration method used by sampleMass is sufficiently computationally complex that it is not useful in the general case. The following Metropolis Hastings sampler pseudocode should be more effective:

set burnin = 50 (say)
foreach sampledParameter
  foreach star
     for (iter = 0; iter < burnin; iter++)
      proposedMass = currentMass + randomStep1
      proposedMassRatio = currentMassRatio + randomStep2
      MHratio = logpost(proposedMass, proposedMassRatio, sampledParameter) - 
                     logpost(currentMass, currentMassRatio, sampledParameter)
      if (log(randomUniform) < MHratio)
        currentMass = proposedMass
        currentMassRatio = proposedMassRatio
      endif
    endfor
  endforeach
endforeach

Implemented in 821c30e554582997e1beb5ec045ad7dce77e2851.