sir.model
Class SIRState

java.lang.Object
  |
  +--sir.model.SIRState
All Implemented Interfaces:
org.omegahat.Simulation.MCMC.MCMCState, Uncmin_methods

public class SIRState
extends java.lang.Object
implements org.omegahat.Simulation.MCMC.MCMCState, Uncmin_methods

Title: LadyBug
Description: Estimation of SIR parameters. Class corresponds to state of the Markov Chain in the simple SIR example. It consists of beta, gamma, I1, \mathbf{I}. The two latter are stored into the same I array as I[0] and I[1]...I[?], respectively.
Copyright: Copyright (c) 2002

Version:
2 April 2002 Added: Implements Uncmin_methods which will contain a wrapper for the optimization class.
Author:
Michael Höhle <hoehle@dina.kvl.dk>

Field Summary
static int accepted
          Static stuff to keep track of progress
static int betaNaccept
          Static stuff to keep track of progress
static int[] betaaccept
          Static link to an array structure holding individual beta acceptance rates.
static Data data
          Static link to the data.
 double deviance
          Deviance part
static int[] gammaaccept
          Gamma acceptance rate
 double loglik
          LogLikelihood
 boolean nogen
          Boolean telling whether genXYZ has been ran on this object (inits events array)
static int runs
          Static stuff to keep track of progress
 
Constructor Summary
SIRState()
          Constructor - necessary to allocate an array for I
SIRState(Event[] I)
          Constructor - necessary to allocate an array for I
 
Method Summary
 void add(SIRState s, double noOfSamples)
          Add values for posterior mean calculations
 double allBetaBetweenIntegrate(int r)
          Calculate int_{S[i-1]},S[i]} sum{u\in U} \sum\sum \beta^N X^u Y^u for a regime i, i.e.
 double allBetaIntegrate(int r)
          Calculate int_{S[i-1]},S[i]} \lambda_{all}(t) for a regime i, i.e.
 double allBetaWithinIntegrate(int r)
          Calculate int_{S[i-1]},S[i]} sum{u\in U} \beta X^u Y^u for a regime i, i.e.
 double f_to_minimize(double[] theta)
          In case the object is given to a numerical minimizer, this is the function to minimize.
 int genUniformI(java.lang.Double I0, Distributions prob)
          Method to generate an instance of the I vector conditioned on tau.
protected  void genXYZFromLambdaI()
          Generate X, Y, XY, and time by sorting appropriate events.
 void gradient(double[] x, double[] g)
          Gradient function, in case one wants to do numeric optimization of the likelihood.
 void hessian(double[] x, double[][] h)
          Hessian function, in case one wants to do numeric optimization of the likelihood.
 boolean isIValid()
          Verify whether I is valid wrt. tau, i.e. that the constraints are fulfilled
protected  int ISize(Unit u, int r)
          Deduce number of 'I' events within unit and regime time contributing to the likelihood.
protected  int ISize(Unit u, int r, Unit uorigin)
          Deduce number of 'I' events within unit u and regime r contributing to the likelihood originating from section u'.
 java.lang.String ItoString()
          Show I as one line
 double logflambdaI()
          Corrected spatial edition, calculating \int_{I_0}^T \sum_{u\in U} \lambda_\all^u(t)
protected  int tauSize(Unit u, int r)
          Deduce number of 'R' events within unit and regime time contributing to the likelihood.
 java.lang.String toString()
          Debug/Information as usual
 void translateTheta2Beta(double[] theta)
          Function to transfer beta vector from one long array to our 3D array This is used in the minimizing stuff.
 void writeRFile()
          To R-compatible file, i.e. write a .txt file containing the following columns t, X_11, Y_11, ..., X_nn, Y_nn, eventype.
protected  double XYintegrate(Unit u, int r)
          Calculate int_{S[i-1]},S[i]} X(t) Y(t) for a regime i.
protected  double Yintegrate(Unit u, int r)
          Calculate int_{S[i-1]},S[i]} Y_u(t) dt for a unit u and regime r.
 
Methods inherited from class java.lang.Object
clone, equals, finalize, getClass, hashCode, notify, notifyAll, wait, wait, wait
 

Field Detail

loglik

public double loglik
LogLikelihood


nogen

public boolean nogen
Boolean telling whether genXYZ has been ran on this object (inits events array)


deviance

public double deviance
Deviance part


runs

public static int runs
Static stuff to keep track of progress


accepted

public static int accepted
Static stuff to keep track of progress


betaNaccept

public static int betaNaccept
Static stuff to keep track of progress


betaaccept

public static int[] betaaccept
Static link to an array structure holding individual beta acceptance rates.


gammaaccept

public static int[] gammaaccept
Gamma acceptance rate


data

public static Data data
Static link to the data.

Constructor Detail

SIRState

public SIRState()
Constructor - necessary to allocate an array for I


SIRState

public SIRState(Event[] I)
Constructor - necessary to allocate an array for I

Method Detail

translateTheta2Beta

public void translateTheta2Beta(double[] theta)
Function to transfer beta vector from one long array to our 3D array This is used in the minimizing stuff.


f_to_minimize

public double f_to_minimize(double[] theta)
In case the object is given to a numerical minimizer, this is the function to minimize. We will only do maximization of the likelihood, i.e. parameters have to be transmoglyffed to account for different indexing of arrays etc. Thus, this will work as a wrapper to logflambdaI.

Specified by:
f_to_minimize in interface Uncmin_methods

gradient

public void gradient(double[] x,
                     double[] g)
Gradient function, in case one wants to do numeric optimization of the likelihood. Contains an EMPTY body.

Specified by:
gradient in interface Uncmin_methods
Parameters:
x - The parameter (remember Fortran indexing of arrayx (starting at 1))
g - The computed gradient.

hessian

public void hessian(double[] x,
                    double[][] h)
Hessian function, in case one wants to do numeric optimization of the likelihood. Contains an EMPTY body.

Specified by:
hessian in interface Uncmin_methods
Parameters:
x - The parameter (remember indexing of arrayx (starting at 1))

writeRFile

public void writeRFile()
To R-compatible file, i.e. write a .txt file containing the following columns t, X_11, Y_11, ..., X_nn, Y_nn, eventype. Output is saved to the file R-epidemic.txt. Note: Regime switches are currently NOT supported in the R-procedure.


ItoString

public java.lang.String ItoString()
Show I as one line


toString

public java.lang.String toString()
Debug/Information as usual

Overrides:
toString in class java.lang.Object

add

public void add(SIRState s,
                double noOfSamples)
Add values for posterior mean calculations


Yintegrate

protected double Yintegrate(Unit u,
                            int r)
Calculate int_{S[i-1]},S[i]} Y_u(t) dt for a unit u and regime r.

Parameters:
u - The unit to get Y(t) from.
r - The regime defining the time limits. I.e. all points with that number.
Returns:
The integral.

allBetaBetweenIntegrate

public double allBetaBetweenIntegrate(int r)
Calculate int_{S[i-1]},S[i]} sum{u\in U} \sum\sum \beta^N X^u Y^u for a regime i, i.e. integrate over all between unit hazards (coz they have a common beta^N) This function could use allBetaWithin, but direct implementation is faster.

Parameters:
r - The regime defining the time limits. I.e. all points with that number.

allBetaWithinIntegrate

public double allBetaWithinIntegrate(int r)
Calculate int_{S[i-1]},S[i]} sum{u\in U} \beta X^u Y^u for a regime i, i.e. integrate over all within unit hazards (coz they have a common beta) This function could use XYIntegrate, but I think direct implementation is faster.

Parameters:
r - The regime defining the time limits. I.e. all points with that number.

allBetaIntegrate

public double allBetaIntegrate(int r)
Calculate int_{S[i-1]},S[i]} \lambda_{all}(t) for a regime i, i.e. integrate over the total hazard.

Parameters:
r - The regime defining the time limits. I.e. all points with that number.

XYintegrate

protected double XYintegrate(Unit u,
                             int r)
Calculate int_{S[i-1]},S[i]} X(t) Y(t) for a regime i.

Parameters:
u - Unit to take X,Y values in. This only works as long as there is no interaction.
r - The regime defining the time limits. I.e. all points with that number.

ISize

protected int ISize(Unit u,
                    int r)
Deduce number of 'I' events within unit and regime time contributing to the likelihood. This means that the I_0 event does not count. Same as calculating |I^u_r|

Parameters:
u - The unit
r - The regime
Returns:
|I^u_r|

ISize

protected int ISize(Unit u,
                    int r,
                    Unit uorigin)
Deduce number of 'I' events within unit u and regime r contributing to the likelihood originating from section u'.

Parameters:
u - The unit
r - The regime
uorigin - The origin of the infection (a reference)
Returns:
|I^u_r|

tauSize

protected int tauSize(Unit u,
                      int r)
Deduce number of 'R' events within unit and regime time contributing to the likelihood. This means that the I event does not count. Same as calculating |tau^u_r|. Instead of going through all events it here suffices to go through the local unit events.

Parameters:
u - The unit
r - The regime
Returns:
|I^u_r|

isIValid

public boolean isIValid()
Verify whether I is valid wrt. tau, i.e. that the constraints are fulfilled


genXYZFromLambdaI

protected void genXYZFromLambdaI()
Generate X, Y, XY, and time by sorting appropriate events. Note: New Version using Event classes. Before, we sorted by zipping. Now we use general sort on all events. Slower? Probably - use zipper again, which required I to be sorted, but then we have to sort it twice??


logflambdaI

public double logflambdaI()
Corrected spatial edition, calculating \int_{I_0}^T \sum_{u\in U} \lambda_\all^u(t)


genUniformI

public int genUniformI(java.lang.Double I0,
                       Distributions prob)
Method to generate an instance of the I vector conditioned on tau. I.e. I has to be valid wrt. tau. Currently rejection sampling is used for the sampling procedure - should be made more efficient.

Parameters:
I0 - Initial value, if not null then it is kept fixed and used to generate I.
prob - The random number generator
Returns:
Number of tries, before a valid sample could be generated.