• A Change-Point Analysis of a Randomized Trial on the Effects of Calcium Supplementation on Blood Pressure by Lawrence Joseph, David B. Wolfson, Roxane du Berger and Roseann M. Lyle (1995) in the book Bayesian Biostatistics, edited by Don A. Berry and Dalene K. Stangl, published by Marcel Dekker.
• Medical applications of the multi-path change-point model by Lawrence Joseph, David B. Wolfson, Roxane du Berger and Roseann M. Lyle (1994)

Abramowitz, M., and Stegun, I.A. (ed.) (1972)
Handbook of Mathematical Functions
New York: Dover Publications Inc.


Cheng, R.C.H (1978)
Generating Beta Variates with Nonintegral Shape Parameters
Communications of the ACM, 21:317-322
(Algorithms BB and BC).


Press, W.H., Flannery, B.P., Teukolsky, S.A., and
Vetterling, W.T. (1990)
Numerical Recipes, The Art of Scientific Computing (Fortran Version)
New York: Cambridge University Press.


Press, W.H., Flannery, B.P., Teukolsky, S.A., and Vetterling, W.T. (1992)
Numerical Recipes in Fortran, The Art of Scientific Computing, Second Edition
New York: Cambridge University Press.

Ex.data       README.MPCP   itp.f         postbet.f     util.f
Ex.input      gen.f         mpcp.f        setup.f
Makefile      header.f      post.f        upd.f

f77 -o mpcpn *.f

f77 -o mpcpp *.f

%cat -n Ex.input
1  37 10 530           !m rows, n columns, seed
2  10000 7000          !# cycles, last cycles retained
3  1000 1 500          !screenview,filedumps,maxfsize(kbytes)
4  f  -99999           !missing (logical) + (-BIG) indicator
5  f   99999           !sample size for each row + (BIG) indicator
6  13 3 2 3 32 23 10 20 22 25       !row indices = min(10,m)
7  10 4 3 1  5  6  7  8  9  9       !column indices = min(10,n)
8  mp.g1                !Gamma  parameter updates - distribution 1
9  mp.g2                !Gamma  parameter updates - distribution 2
10  mp.n1                !Normal parameter updates - distribution 1
11  mp.n2                !Normal parameter updates - distribution 2
12  mp.al                !Dirichlet parameter updates
13  mp.ad                !avg (iterations) dirichlet updates
14  mp.ga                !avg (iterations) E(X),X~Gamma
15  mp.am                !avg E(X2)-E(X1),E(X)=updated Normal mean
16  mp.ap                !Pr(no change)
17  16 16                !m prior gamma shapes
18  16 16
:
:
52  16 16
53  16 16
54  .0025 .0025          !m prior gamma scales
55  .0025 .0025
:
:
89  .0025 .0025
90  .0025 .0025
91  115 115 115 115 115 115 115 115 115 115
92  115 115 115 115 115 115 115 115 115 115
93  115 115 115 115 115 115 115 115 115 115
94  115 115 115 115 115 115 115           !m exp mean dist1
95  115 115 115 115 115 115 115 115 115 115
96  115 115 115 115 115 115 115 115 115 115
97  115 115 115 115 115 115 115 115 115 115
98  115 115 115 115 115 115 115           !m exp mean dist2
99  .25 .25                               !omega 1 and 2
100  .0 .0 .0 .1 .1 .1 .1 .1 .1 .1         !n dirichlet parameters

• Line 1: The first 2 items are number of rows <37> and number of columns <10> of the data matrix; item 3 is the seed number <530> for the random number generator used in the Gibbs sampler.
• Line 2: Number of Gibbs cycles to run <10000> with the last 7000 cycles retained for inference. Note that either or both items set to zero results in an error.
• Line 3: Item 1 sets the cycle numbers (every <1000>) at which the user views intermediate results and is useful for inspection. Item 2 (k = 1 filedumps) requires that each of the last 7000 cycles requested in item 2 of line 2, be used for inference and output to file. If k > 1, then every kth information cycle is used for inference and output to file. if k < 0 or > 7000, k is reset to 1. Finally, because the Bayesian multiple-path change-point method can generate huge amounts of data, item 3 imposes a size limit <500> (approximately 500 kbytes) to files. If the size limit is very small with respect to the number of cycles and the number of filedumps, then only summary files are created.
• Lines 4 and 5: Line 4 is a logical missing value indicator. If set to true (`t'), the program requires that the missing value be a negative number whose value is less than the minimum value in the data matrix. For example, the minimum value in this data matrix is 93, so that if there were missing data we can safely set the missing value indicator to -1 and insert in the matrix element the value of -1. Line 5 is the sample size at each data matrix row. If the sample size varies from row to row -- the logical indicator is set to `t', we choose a value of the variable sample size indicator which exceeds the maximum value in the data matrix. Suppose the maximum value in this data set is 140, and that row 23 has a sample size of 8 instead of 10, and row 31 has a sample size of 5. We choose a value which exceeds 140 (i.e. 141 or 99999) and position this value at columns 9 and 6, respectively for rows 23 and 31. Note that each row must be of fixed length, so that columns 9 and 10 of row 23 and columns 6 to 10 of row 31 must be filled; the simplest is to set the rest of the columns to 141 or less.
• Lines 6 and 7: Line 6 is the set of row indices required for file output of the Gamma and Normal parameter updates. For example, the first element of line 6 is 13, which means the selection of the 13th row. Choosing a row index greater than 37 (see Line 1) results in an error. Line 7 is the set of column indices required for file output of the Dirichlet parameter updates. For example the fourth element of line 7 is 1, which means the selection of the 1st column. Choosing a column index greater than 10 (see Line 1) results in an error. If the data matrix contains more than ten rows/columns, the user selects a maximum of ten row/column indices in any order. On the other hand, if the data matrix contains less than ten rows/columns - say eigth rows/columns - the user selects eight row/column indices in any order. These criteria are necessary to limit the space requirements of the parameter update files.
• Lines 8 to 16: the program outputs a total of nine files associated to Normal data. The first five contain information at each selected cycle for:
  1) Gamma and Normal parameter updates, before and after change-point, at selected rows; and
  2) Dirichlet parameter updates at selected columns.
  The remaining four files contain parameters averaged over selected cycles:
  1) Dirichlet updates at each data column;
  2) Expected values of Gamma parameter updates before and after change-point, at each data row;
  3) Difference between two expected values before and after change-point of Normal parameter updates, at each data matrix row; and
  
  4) The probability of no change is the last column value of the multinomial generation from Dirichlet parameter updates. This information is useful to compute the probability that subjects in one group are more likely to change than subjects in another group using a Mann-Whitney-type statistic.
• Lines 17 to 53: To each row of Ex.data and to each distribution <1 and 2>, set Gamma shape prior parameters.
• Lines 54 to 90: To each row of Ex.data and to each distribution <1 and 2>, set Gamma scale prior parameters.
• Lines 91 to 94: To each row of Ex.data , assign prior data mean values for distribution 1.
• Lines 95 to 98: To each row of Ex.data , assign prior data mean values for distribution 2.
• Line 99: To each distribution, assign a weighting factor, used for Gamma and Normal parameter updates.
• Line 100: Assign to each column of Ex.data , a prior Dirichlet parameter where, by definition, each is greater than zero. If a prior Dirichlet parameter is set to zero at column k, than parameter updating is not performed from columns 1 to k. Note that an error results if a prior Dirichlet parameter is set to zero at the last data matrix column. SPECIAL NOTE: an error results if the variable sample size at a given row is less than k (see description of Line 5).

%cat -n Ex.input
1  285 8 582           !m rows, n columns, seed
2  10000 3000          !# cycles, last cycles retained
3  5000 1 500          !screenview,filedumps,maxfsize(kbytes)
4  t  -99              !missing(logical) + (-BIG) indicator
5  f  110              !sample size for each row + (BIG) indicator
6  246 45 46 76 98 32 54 18 222 271 !row indices = min(10,m)
7    8  4  7  6  5  3  1  1  10   8 !column indices = min(10,n)
8  mp.g1                !Gamma parameter updates - distribution 1
9  mp.g2                !Gamma parameter updates - distribution 2
10  mp.al                !Dirichlet parameter updates
11  mp.ad                !avg (iterations) dirichlet updates
12  mp.ga                !avg (iterations) E(X),X~Gamma
13  mp.am                !avg E(X2)-E(X1), X~Gamma
14  mp.ap                !Pr(no change)
15      1.00000    1.00000           !m prior Gamma shapes
16      1.00000    1.00000
:
:
300      15.0000    15.0000           !m prior Gamma scales
301      15.0000    15.0000
:
:
583      15.0000    15.0000
584      15.0000    15.0000
585  0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0  !n dirichlet parameters

• Lines 1 to 7 are defined as in the Normal format.
• Lines 8 to 14: the program outputs a total of seven files associated to Poisson data. The first three contain information at each selected cycle and row or column for:
  1) Gamma parameter updates, before and after change-point, at selected rows; and
  2) Dirichlet parameter updates at selected columns.
  The remaining four files are parameters averaged over selected cycles:
  1) Dirichlet updates as in Normal;
  2) Expected values of Gamma parameter updates described in the Normal format;
  3) Difference between two expected values before and after change-point of Gamma parameter updates, at each data matrix row; and
  4) The probability of no change described in the Normal format.
• Lines 15 to 584 are the Gamma prior parameter information described in the Normal format.
• Line 585 is equivalently described in Line 100 of the Normal format.

% cat Ex.input >! mpc.par    #required parameter file name
% cat Ex.data >! mpc.dat     #required data file name

% mpcpn

% mpcpp

# Number of rows and columns in data matrix
nrowd _ as.numeric(unix("head -1 mpc.par | awk '{print $1}'"))
ncold _ as.numeric(unix("head -1 mpc.par | awk '{print $2}'"))
# Dirichlet parameter priors
unix("egrep '!n dirichlet' mpc.par | awk '{print}' > tmp")
diric _ scan(file='tmp',n=ncold) ; unix("rm tmp")
# Number of rows and columns output to iterative files
#  (assuming they exist)
unix("dat1=`egrep 'Gamma parameter updates - distribution 1' mpc.par | \
awk '{print $1}'` ; tail -1 $dat1 > tmp")
nrowi _ length(scan(file='tmp'))/3 ; unix("rm tmp")
unix("dat1=`egrep 'Dirichlet parameter' mpc.par | \
awk '{print $1}'`; tail -1 $dat1 > tmp")
ncoli _ length(scan(file='tmp')) ; unix("rm tmp")

labdir _ c('mp.al')
colsel _ scan(file=labdir,skip=7,n=ncoli)   # Column index
marpost _ matrix(scan(file=labdir,skip=18),byrow=T,
ncol=ncoli)
rsum _  nrowd+sum(diric)         # Each row sum is identical
marpost _ marpost/rsum           # Divide each element by row sum
marpost _ apply(marpost,2,mean)  # Take the mean of each column
marpost _ marpost[order(colsel)] # ``Natural'' column order

labdir _ c('mp.ad')
marpost _ scan(labdir,skip=6)

openlook()                       # Or any valid graphics device
labx _ c('measurement period')
laby _ c('change-point \n marginal posterior mean')
plot(marpost,ylim=c(0,max(marpost)+0.05),xaxt='n',
type = 'n', lty = 1, xlab = labx, ylab = laby)
axis(1,at=seq(1,length(marpost),by=1))
for(j in seq(along = marpost)) {
segments(j-0.2, 0, j-0.2, marpost[j])
segments(j-0.2,marpost[j],j+0.2,marpost[j])
segments(j+0.2,marpost[j],j+0.2,0)
}

labdir _ c('mp.al')
xduf _ matrix(scan(file=labdir,skip=18),
byrow=T,ncol=ncoli)
ba _ xduf[,2]   # Second column index - 4th data column

dec _ 100                     # Number of domain points
xq _  seq(0.00001,0.99999,length=dec)
rsum _ nrowd+sum(diric)
bb _ rsum - ba
fxbet _ post.beta(xq,ba,bb)   # Mean Beta Density
post.beta _ function(x,a,b){
# Compute beta distribution at each
#  domain point x averaged over beta
#  parameter vectors a,b
# Assume that Fortran procedure is in dyn.load
res _ .Fortran('postbet',
m = as.integer(length(a)),
n = as.integer(length(x)),
fk = as.double(x),
as.double(x),
as.double(a),
as.double(b))$fk
return(res)
}

openlook()                   # Or any valid graphics device
labx _ c('')
laby _ paste(c('marginal posterior density \n'),
c('given Dirichlet column updates '))
plot(xq[-c(1,length(xq))],fxbet[-c(1,length(fxbet))],
type = 'l',lty = 1, xlab = labx, ylab = laby)

labnoch1 _ c('mp.ap1')
labnoch2 _ c('mp.ap2')
noch1 _ scan(labnoch1,skip=9)
noch2 _ scan(labnoch2,skip=9)
prnoch _ prob.nochange(noch1,noch2)

prob.nochange _ function(pn1,pn2)
{
# Mann-Whitney statistic wm
# wm = wilc - term, wilc = Wilcoxon,
# term = 0.5*k*(k+1), k=vector length of both pn1,pn2
# For large k, wm = wilc*sqrt(Var(wilc)) + e(wilc) - term,
# wilc=uncorrected Z,
# Pr(Pr(pn1)<Pr(pn2)) = 1 - wm/k^2
k <- length(pn1)
if (k!=length(pn2)) stop('sample sizes different')
term <- 0.5 * k * (k + 1)
wilc <- wilcox.test(pn1, pn2, correct = F)$statistic
if(names(wilc) == "rank-sum statistic W") {
wm <- wilc - term
}
else {
ewilc <- 0.5 * k * (2 * k + 1)
svars <- sqrt((1/12) * k^2 * (2 * k + 1))
wm <- wilc * svars + ewilc - term
}
pwm _ 1 - wm/k^2
return(ifelse(pwm>=0&pwm<=1,pwm,NA))
}

labgam1 _ c('mp.g1')
labgam2 _ c('mp.g2')
rowsel _ scan(file=labgam1,skip=9,n=nrowi)          # Row index
k _ 8                    # Choose row index element that exists
if (k>length(rowsel)|k<1) stop('\n \n ****** bad index element')
ind _ (k-1)*3 + 1                       # Matrix column element
#Gamma updates for Distributions 1 and 2
gam1 _ matrix(scan(file=labgam1,skip=20),byrow=T,
ncol=3*length(rowsel))
gam2 _ matrix(scan(file=labgam2,skip=20),byrow=T,
ncol=3*length(rowsel))
#Gamma shape and scale - Distribution 1
a1 _ gam1[,ind]
b1 _ gam1[,ind+1]
#Gamma shape and scale - Distribution 2 (change-points only)
a2 _ gam2[gam2[,ind+2]>0,ind]
b2 _ gam2[gam2[,ind+2]>0,ind+1]
#Domain determination (*****to be refined to taste*****)
#Care should be taken in choosing domain points. In the
#beginning, use only a limited number of domain points.
#Once the location of f>0 is found, add more domain points to
#obtain a smoother plot. This is necessary, because Splus
#can run into memory allocation problems.
dec _ 20
xq1 _ seq(0,1,length=dec)
xq2 _ seq(1,8,length=dec)
xq _ unique(c(xq1,xq2))
#Compute Gamma densities
fx1 _ dens.gam(xq1,a1,b1)
fx2 _ dens.gam(xq2,a2,b2)

dens.gam _ function(x,a,b) {
# compute Gamma density whose domain is x and
# vector elements parameters a,b
# Output is mean density at each domain point.
mrows _ length(x)
ncols _ length(a)
fmat _ t(matrix(x,mrows,ncols))
fmat1 _ apply(fmat,2,f.gam,a,b)
return(apply(fmat1,2,mean))
}
f.gam _ function(x,a,b) {
#Divide the domain points by scale b to make
#density function scaleless (i.e. dgamma(a,1))
#Divide result by scale b to make density function
#scaled (i.e. density gamma(a,b))
#return(dgamma(x/b,a))   #scaleless result
return(dgamma(x/b,a)/b)  #scaled result
}

openlook()                       # Or any valid graphics device
labx _ paste('observation number',rowsel[k])
laby _ paste('probability density')
plot(c(xq[1],xq[length(xq)]),c(0,max(fx1,fx2)+0.01),
type ='n', lty = 1, xlab = labx, ylab = laby)
lines(xq1,fx1,lty=6)  # Rate before change
lines(xq2,fx2,lty=1)  # Rate after change
leg.names _ c('Rate before change',
'Rate after change')
legend(max(xq)*0.5,max(fx1,fx2)*0.5,leg.names,
bty="n", lty=c(6,1))

1. This note is not to be removed;
2. Publications using mpcpn or mpcpp results should reference the manuscripts mentioned above.