Readme file
SERIES C
 Applied Statistics
A continuous latent spatial model for crack initiation in bone cement, by E. A. Heron and C. D. Walsh, pages 25–42;
Description of Program
----------------------------
This program uses Markov chain Monte Carlo techniques to estimate the
parameters of the model specified in the above paper. See the paper for
a detailed description of the algorithm. Simulated data are provided.
*******************************************************
* Copyright notice
* ------------------
* These programs are available in the public domain but WITHOUT ANY
* WARRANTY.
* If this software is used in work presented for publication, we ask that
* the above paper be referenced.
*
* Dec. 2007
*
* Thank you.
********************************************************
Input Files
----------------------------
****** DATA
The following data files containing two - SIMULATED - data sets are
supplied together with the stress values at these locations.
cracks_x1.txt
cracks_y1.txt
these contain the x and y coordinate values of the crack locations of the
first set of simulated data and
cracks_x2.txt
cracks_y2.txt
these contain the x and y coordinate values of the crack locations of the
second set of simulated data.
stress_1.txt
stress_2.txt
contain the stress values at the crack locations in the first and second
specimens respectively.
****** GRID APPROXIMATION
In order to calculate the integrals we used a discretisation. We divided each
window into a fine grid and calculated the stress value at the centroid of
each grid segment. Thus each segment on the grid has a vector (C_jg, T_jg),
j indicating medial or lateral window and g indicating which grid segment.
One of C_jg or T_jg will be zero as there cannot be both a compression and
a tension value at a single location.
The approximations are as follows:
int_X_j(C(l)w(dx)) is given by Sum(C_jg*A_jg), where int_X_j is the integral
over X_j and A_jg is the area of the grid segment g in window j. Similarily
int_X_j(T(l)w(dx)) is given by Sum(T_jg*A_jg).
The following files are provided for this discretisation:
centroidX.txt
centroidY.txt
containing the x and y coordinates of the grid,
centroidgridarea.txt
containing the areas of the grid segments, and
centroidgridstress.txt
containing the stress values at the centroids of the grid segments.
To Run the Program
----------------------------
The program has only been run on a unix platform. All files should be stored in
the same directory. A makefile is included called "Makefile" for compiling and
linking the program. To run the program, at the prompt type the following:
1. make clean
and then press return. This command removes any old .o files if they are
there, if there are no .o files then running this command does not cause
any problems. Then type
2. make
and press return. Then type
3. main
and press return.
Files generated by the Program
--------------------------------
When the program is running at each iteration output is written to the following
file
out.txt
which stores the following:
1st column: beta_0 values
2nd column: beta_1 values
3rd column: rho values
Grid1.txt
Grid2.txt
These two files store the posterior mean of the Gamma random fields, 1 and
2 indicating the specimens. See further on for R code that uses this
output to plot the Gamma random fields. These files each contain five
columns of data. First four columns are coordinates for a grid and the
fifth column contains the values of the Gamma random field for each of the
grid segments.
Program Files
--------------------------------
The following program files are included:
******* main.c
This file contains the main function that implements the algorithm
described in the paper by calling functions in the other files included.
******* Input.c
The function in this file inputs the x and y coordinates of the crack
locations for the simulated data sets together with the stress values
at these locations. The stress values are then divided into compression
and tension stresses according to whether they are positive or negative.
The coordinates of the grid used in the discretisation is inputted
along with the stress values at the centroids and the areas of the
grid segments. The stress values for the grid are also divided into
compression and tension stresses according to whether they are postive
or negative.
******* Priors.c
The parameters of the prior distributions for beta_0, beta_1, rho and
the Gamma random field are initialised in the function in this file.
******* Bits.c
Using Bernoulli random variables the function in this file attributes each
of the cracks individually to either compression, tension or the latent factor.
If a crack is attributed to a latent factor this crack location is stored
for later use in the algorithm.
******* Latent.c
The function "latent(long, idum)" in this file simulates the Gamma random field
using the Inverse Levy Measure algorithm as described in the paper. This file also
contains two other functions "CumSum(double, double, int)" and "invE1(double)"
that are used in the simulation of the Gamma random field.
******* Latent_Cracks.c
For every crack that is attributable to the Gamma random field, the function in
this file chooses a corresponding point in the S space according to some
probabilities. For each of these locations a random walk Metropolis step is
then carried out, this is to ensure that the same points in the S space are not
repeatedly chosen.
******* Kernel.c
Given the locations (x1, y1) and (x2, y2) and the parameter rr, the function
in this file returns the Gaussian kernel.
******* Rho.c
The function in this file performs a Gaussian random walk Metropolis step in
order to update the parameter rho. The standard deviation of the proposal
distribution can be adjusted in this file.
******* Grid.c
At a given iteration of the MCMC algorithm the function in this file sums up
the Gamma Random Field components that fall into the grid for the ouput of the
Gamma Random Field. This function is used to obtain the posterior mean of Gamma
random fields for plotting later.
******* Bern.c
Given the input parameter p, the function in this file returns a Bernouli random
variable.
******* Exp.c
The function in this file generates an Exponential random number.
******* Gamma.c
The function in this file generates a Gamma random number with parameters:
shape = alpha and scale = 1.
******* Ran2.c
The function in this file generates a Uniform random number between 0 and 1.
******* Normal.c
The function in this file generates a Normal random number with mean = m and
standard deviation = sd.
******* H.h
This is a header file that contains variable and function declarations.
R Code
--------------------------------
The file "GRF_Plotting.r" contains R code to plot the posterior
mean of the Gamma random fields that are produced by the MCMC algorithm.
This code relies on the use of the R package "splancs". Two text files for
plotting the polygons surrounding the windows are also supplied in "polygonMed_coords.txt" and "polygonLat_coords.txt"
References
--------------------------------
The following references are provided for the algorithms for the random number
generators that are used in the code. The R package "splancs" is used in
the R code for plotting the Gamma random field.
- Press W. Flannery B. Teukolsky, S. and Vetterling, W. "Numerical Recipes
in C: The Art of Scientific Computing". Cambridge University Press, 2nd
edition, 1992.
- Ripley, Brian D. "Stochastic Simulation", New York; Chichester: Wiley,
1987.
- The splancs Package for R: Spatial and Space-Time Point Pattern Analysis,
Authors: Barry Rowlingson and Peter Diggle, adapted and packaged for R by
Roger Bivand.
Elizabeth A. Heron
Department of Statistics
University of Warwick
Coventry
CV4 7AL
UK
E-email: e.heron@warwick.ac.uk
|
