Create a posterior log-pdf function for monte carlo sampling the model parameters.

The function returned takes a vector of parameter values as an argument and returns a the log of the posterior pdf calculated by comparing this data to the supplied observed data. There is an optional second argument that is described further below. The log-post function should be usable in any MC sampler. It will also work in an optimizer, provided you configure it to look for a maximum rather than a minimum.

mc.setup(
  obsdata_filename,
  logprior = NULL,
  logfile = NULL,
  chunksize = 10,
  trace_param = FALSE
)

Arguments

obsdata_filename

File name for the oberved data.
logprior

Optional function that takes a parameter vector and returns a log prior probability density. It need not (and should not) support the log-posterior function's optional argument described in the details section.
logfile

Optional file name for logging.
chunksize

See details for description of chunking and why it is needed.
trace_param

Set this to TRUE to generate the data file (.dat file) during the statistical fitting procedure.

Value

A function that computes the log-posterior probability density for an input vector of parameters.

Details

The observed data should have the following columns: Ps, Pn, Y, Qs, Qn, sigQs, sigQn. There is a special provision to process the dataset produced for GCAM

The demand system employs a nonlinear equation solver, which is used to solve for the budget fractions that appear in the demand equations. Although structured as a system of equations, each time is actually independent of all the others, but the solver has no way of knowing this. Because the complexity of the solver scales nonlinearly with the number of equations, we break the times in the dataset up into chunks of modest size.

The optional second argument to the log-posterior function produced by this function gives the number of parameter sets that are concatenated into a single vector. This capability is supplied to make life easier back when we were calling these functions from a C code that used SSE/AVX instructions to vectorize the sampler. Samplers written in R can and should pretend that this argument doesn't exist.

TODO

  • Provide an example dataset for the observed data.

  • Dispense with this chunk business and just run a 1-d solver in a loop. What was I thinking?