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