Code for implementation of OD etimation
Claudia Tebaldi
Sat, 24 Jan 2009 13:45:19 -0700
Three files:
- "metronet.for" is a fortran 77 program that produces the MCMC simulation
of OD counts.
- "input" is a template for the input file that the program requires.
- "guess.s" contains several functions written in Splus that allow you to
prepare the necessary input to the problem, given only your OD matrix A and
the vector Y of link counts.
The content of "input" is consistent with the 7by12 problem that was
presented in Tebaldi&West (1998).
In more details, the fortran program needs as input
K the number of links
Y the vector of link counts
N the number of OD-pairs
X2 an initial guess for the part of the vector of OD counts that is "free"
X1 the derived part of the OD vector that, given X2, satisfies the link counts in Y
A1 the invertible (K by K) part of the matrix A, after pivoting, as
explained in the article, for example by a QR decomposition
A2 the remaining part of the matrix, K by N-K
invA1 the inverse ofthe matrix A1
RECORD the number of simulated set of values you want to save in the output files.
When presented with a network and link counts you of course have only A
and Y.
>From them you need to derive A1, A2 and an initial guess for
X1 and X2 (see article if this notation doesn't ring a bell immediately).
That is what the Splus functions in the file "guess.s" help you to do.
They are not strictly necessary: you can apply a QR decomposition to A by
any other numerical software of your choice, and you can try to get a
feasible initial guess for the vector X (decomposed in X1 and X2) by other
means than the functions I provide...but if you are familiar with Splus
(actually even better with R (http://www.r-project.org/) which is its
freely available version) you can source the file "guess.s" and work out
the ingredients to the problem by just
1) entering the matrix A and the vector Y that pertain to your problem
2) applying the function get12(A), thus extracting A1, A2 and the pivot
information
3) recovering an initial guess for X2 by applying either the function
guesslow(A1,A2,Y) or guessup(A1,A2,Y)
4) recovering the vector X1 by simula.x1(A1,Y,A2,X2)
4) writing a new "input" file with all the ingredients required by the
fortran program
5) compiling and running metronet.for
Note that there are a few lines commented out in metronet.for that -- if
uncommented and by commenting out the lines that read in the input from
the file -- allow you to enter the information at the prompt, when running
the program.
Also note that the vector X that gets written out by the MCMC simulation
will have the order of the pivot, not the original order of the columns of
the matrix A.
--
Claudia Tebaldi