sir.estimate
Class LadyBug2

java.lang.Object
  |
  +--org.omegahat.Simulation.MCMC.NotifyingMCMCObject
        |
        +--org.omegahat.Simulation.MCMC.BaseMarkovChain
              |
              +--sir.estimate.LadyBug2
All Implemented Interfaces:
org.omegahat.Simulation.MCMC.MarkovChain, org.omegahat.Simulation.MCMC.NotifyingObject, java.lang.Runnable

public class LadyBug2
extends org.omegahat.Simulation.MCMC.BaseMarkovChain

This class implements estimation of the parameters in a SIR epidemic. Four methods are supported:

  1. MCMC analysis as described in O'Neill & Roberts, JRSSA 162. The results of MCMC using burn-in rate, thinning factor, etc. specified by the code are written to a file log.txt. This can then be used by CODA or BOA for analysis in Splus.
  2. Maximum Likelihood
  3. Maximum Marginalized Likelihood
  4. Simulated ML
The class design unfortunately shows clea signs that this has been an incremental project, i.e. it has grown from nothing to something quite huge, without the structure following the same development. E.g. this means that global variables and switches are used instead of nice design patterns.

Version:
2 April 2002
Author:
Michael Höhle <hoehle@dina.kvl.dk>

Field Summary
protected  Data data
          New version - all data related variables are gathered in a Data object
static boolean DEBUG
          Debugging flag - gives extra printed info
static java.text.NumberFormat form
          For formatting of doubles to an output stream.
protected  boolean gotI
          Do we have acces to I time
protected  boolean metropolisBeta
          Do we want the beta's to be updated using Metropolis sampling.
protected  boolean metropolisGamma
          Do we want the gamma's to be updated using Metropolis, i.e.
protected  Distributions prob
          The Pseudo random generator for distributions
 
Fields inherited from class org.omegahat.Simulation.MCMC.NotifyingMCMCObject
listeners
 
Constructor Summary
LadyBug2(org.omegahat.Simulation.RandomGenerators.PRNG prng, java.lang.String lambdaFile, java.lang.String IFile)
          Constructor which calls all the necessary sub-routines, s.a. reading of data, registering of listeners and iteration of the MC.
 
Method Summary
 org.omegahat.Simulation.MCMC.MCMCState generate(org.omegahat.Simulation.MCMC.MCMCState state)
          Generate the next state of the Markov Chain using Gibbs-within-Metropolis.
 org.omegahat.Simulation.MCMC.MCMCState generate2(org.omegahat.Simulation.MCMC.MCMCState state)
          Generate the next state of the Markov Chain using Gibbs-within-Metropolis.
 org.omegahat.Simulation.MCMC.MCMCState initialState(boolean mcmctheta)
          Method to generate the initial state of the Markov Chain.
static void main(java.lang.String[] argv)
          Main function - start a LadyBug2 analyserer based on prompt arguments
 void MCMCEst(int noOfSamples, int thinFactor, int noOfBurnin, double hmeanscale)
          MCMC Estimation method - generates noOfSamples*thinFactor + noOfBurnin samples and estimates beta and gamma based on recovery times alone.
 SIRState MLEst(SIRState state, boolean printEstimates)
          Calculate parameters beta and gamma using ML estimation.
 void MMLEst(int noOfEstimates, int noOfMCSamples, double sumscale)
          Do maximization of the marginalized density f(\bm{lambda}|\bm{beta}).
protected  void readDataFile(java.lang.String lambdaFile, java.lang.String IFile)
          Read information about the epidemic from the denoted files.
protected  void readLogFile(double[] Ebeta, double[] Egamma, double[] EbetaN, double[] beta, double[] gamma, double[] betaN, double[] pdm)
          Method for calculating mean or variance of samples from the logfile.
 void SimMLEst()
          Sim params and then do ML
 
Methods inherited from class org.omegahat.Simulation.MCMC.BaseMarkovChain
getState, iterate, run, step
 
Methods inherited from class org.omegahat.Simulation.MCMC.NotifyingMCMCObject
notifyAll, registerListener, unregisterListener
 
Methods inherited from class java.lang.Object
clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
 

Field Detail

DEBUG

public static boolean DEBUG
Debugging flag - gives extra printed info


gotI

protected boolean gotI
Do we have acces to I time


data

protected Data data
New version - all data related variables are gathered in a Data object


metropolisBeta

protected boolean metropolisBeta
Do we want the beta's to be updated using Metropolis sampling. Usually the answer is no if data.noOfUnits == 1 and yes otherwise. But it is possible to hardcode it.


metropolisGamma

protected boolean metropolisGamma
Do we want the gamma's to be updated using Metropolis, i.e. RandomWalk


form

public static java.text.NumberFormat form
For formatting of doubles to an output stream.


prob

protected Distributions prob
The Pseudo random generator for distributions

Constructor Detail

LadyBug2

public LadyBug2(org.omegahat.Simulation.RandomGenerators.PRNG prng,
                java.lang.String lambdaFile,
                java.lang.String IFile)
Constructor which calls all the necessary sub-routines, s.a. reading of data, registering of listeners and iteration of the MC. Output is written into the file log.txt.

Parameters:
prng - The pseudo random number generator
lambdaFile - Name of the file containing recovery times
Method Detail

initialState

public org.omegahat.Simulation.MCMC.MCMCState initialState(boolean mcmctheta)
Method to generate the initial state of the Markov Chain. Currently beta, gamma are assigned fixed value and I vector is sampled step by step, s.t. I[i-1] <= I[i] <(=) thau[i-1]. This is ensured by drawing I[i] uniformly on (I[i-1],thau[i-1]) THIS IS THE OLD METHOD

Parameters:
mcmctheta - How to interpretate theta, if true then param in exp distribution.

generate

public org.omegahat.Simulation.MCMC.MCMCState generate(org.omegahat.Simulation.MCMC.MCMCState state)
Generate the next state of the Markov Chain using Gibbs-within-Metropolis. In this new version all the business is done by the state.

Specified by:
generate in class org.omegahat.Simulation.MCMC.BaseMarkovChain
Parameters:
state - Current chain of the MC
Returns:
New state of chain.
See Also:
EpiState

generate2

public org.omegahat.Simulation.MCMC.MCMCState generate2(org.omegahat.Simulation.MCMC.MCMCState state)
Generate the next state of the Markov Chain using Gibbs-within-Metropolis. The components beta, gamma, and I1 are obtained by Gibbs sampling. I2,...,Im is obtained by a Metropolis sampler with proposal distribution with a acceptance probability. The overall construction is. New version using unit information as well.

Parameters:
state - Current chain of the MC
Returns:
New state of chain.
See Also:
SIRState

readLogFile

protected void readLogFile(double[] Ebeta,
                           double[] Egamma,
                           double[] EbetaN,
                           double[] beta,
                           double[] gamma,
                           double[] betaN,
                           double[] pdm)
Method for calculating mean or variance of samples from the logfile. The output is mean if Ebeta, Egamma, and EbetaN are null. Otherwise - if provided we calculate the variance.

Parameters:
pdm - P(D|M) - marginal posterior of data estimated by P2 estimator (harmonic mean of likelihoods)

readDataFile

protected void readDataFile(java.lang.String lambdaFile,
                            java.lang.String IFile)
Read information about the epidemic from the denoted files. User always has to denote a file containing the recovery times, if information about infectious times is available, this can be specified as well. New version used JavaCC grammar to parse the recovery specification file.

Parameters:
lambdaFile - The filename containing the recovery times
IFile - Filename of the infectious time (if it exists).

SimMLEst

public void SimMLEst()
Sim params and then do ML


MMLEst

public void MMLEst(int noOfEstimates,
                   int noOfMCSamples,
                   double sumscale)
Do maximization of the marginalized density f(\bm{lambda}|\bm{beta}). This is done by calculating the integral \int f(\bm{lambda},\bm{I},I_0|\bm{beta}) using Monte Carlo integration.


MLEst

public SIRState MLEst(SIRState state,
                      boolean printEstimates)
Calculate parameters beta and gamma using ML estimation.

Parameters:
state - In case we want to supply infection times directly, we provide them. If null then they are taken from the I-file.

MCMCEst

public void MCMCEst(int noOfSamples,
                    int thinFactor,
                    int noOfBurnin,
                    double hmeanscale)
MCMC Estimation method - generates noOfSamples*thinFactor + noOfBurnin samples and estimates beta and gamma based on recovery times alone.

Parameters:
noOfBurnin - Number of samples to use as burnin. Note: No thinning
thinFactor - Thinning,e.g. take each 100th sample.
hmeanscale - Scaling of \hat{P}_2, i.e. scaling when compuing harmonic mean.

main

public static void main(java.lang.String[] argv)
                 throws java.lang.Throwable
Main function - start a LadyBug2 analyserer based on prompt arguments

java.lang.Throwable
See Also:
readDataFile(java.lang.String, java.lang.String)