At the end of this messageI have pasted code I developed to perform MVNormal
Posterior
Predictive Discriminant Analysis (PPDA). I validated the code by applying
it to a data set found
in G.A.F. Seber's book, "Multivariate Observations", p. 295. The expected
value(s) for the
posterior distributions of the population means matched Seber's analysis
exactly, while the
values obtained for the covariance matrix entries are slightly different,
but still close.
The two "new" observations I used in forming ratios of posterior predictive
distributions were
made up by me.
It is always possible that this code doesn't work like it should, and I
would be especially grateful
for any comments and criticisms. I would be particularly interested if
anyone were to find
errors in the manner I calculated the exponents for the posterior predictive
distributions.
It is my hope that this program benefits those of you who are interested in
Discriminant
Analysis.
Regards,
David Alan Paul, Ph.D.
Battelle Memorial Institute - SDAS
[log in to unmask]
614-424-3176
614-424-4611 (fax)
http://www.battelle.org/statistics
#The ultimate goal is to create a posterior predictive density for each of
two
#populations and use these predictive distributions to form a ratio, on
which we
#shall base the assignment of new observations to one of the two
populations.
#The theory for this approach is given in G.A.F. Seber's book, "Multivariate
#Observations", on pages 292 and 301. (Wiley)
#In this program, N1 represents the total number of multivariate
observations from
#Group1, N2 represents the total number of multivariate observations from
Group2,
#and N3 represents the total number of new observations to be classified.
#The likelihoods being formed correspond to the multivariate normal
distribution,
#and the example data represents the Flea-Beetles data, page 295, Seber (see
#comments before the actual data list for more information).
#It is assumed that the covariance structures of the two populations are
unequal,
#so that our Bayesian analysis corresponds most closely to Quadratic
Discriminant
#Analysis in the frequentist context. Forcing tau to be the same for both
groups
#would not be difficult to do.
#The assignment of new observations to one or the other of the multivariate
#populations is made based on the variable 'criterion'. The classification
rule is
#to assign an observation to Group 1 if the ratio
#(posterior predictive for group 1)/(posterior predictive for group 2)
#is greater than some specified constant, c. In its most general form, this
#constant is given by (p2 * C[1|2]) / (p1 * C[2|1]). Assuming equal prior
#probability of group membership and equal costs of misclassification, c=1.
model
{
for(i in 1:N1)
{
FleaBeetles[i,1:4] ~ dmnorm(mu[1,] ,
tau[1,,])
}
for(i in (N1+1):(N1+N2))
{
FleaBeetles[i,1:4] ~ dmnorm(mu[2,] ,
tau[2,,])
}
#The most difficult part of implementing multivariate predictive posterior
#methods in WinBUGS is that you have to be able to express the posterior
#predictive distribution(s) at the data point(s) of interest - that is,
observed
#values of as-yet unclassified beetles. Since WinBUGS does not yet support
#much in the way of matrix algebra, we are stuck with writing out the
exponent
#of the multivariate normal distribution in longhand. The trick is to be
able
#to use the inprod(v1,v2) function to express this exponent in a convenient
#way. Note:
#
# (1) The vector [ (x[k,1]-mu1), (x[k,2]-mu2), (x[k,3]-mu3), (x[k,4]-mu4)
], where
# the x[k,j] - values come from the new observation(s), is designated
as
# Vec[k,,]. Note that k indicates the relevant population, either 1
or 2.
#
# (2) The columns of the posterior covariance matrices are designated
using
# tau[ population#, row# , column# ]. Thus, tau[1,,2] designates the
second
# column of the posterior precision matrix for population 1.
#
#With these conventions, in mind, the following code calculates the ratios
#of the posterior predictive multivariate normal distributions evaluated at
#the new observations of interest. Of course, the reader will have to
modify
#this code to accommodate data with higher or lower dimensionality. The
form
#of the exponent in a multivariate normal distribution is given as
#(-.5)(x-mu)*V(x-mu) where * indicates a transpose, and V represents the
#precision matrix. Note that V=Inverse(Covariance Matrix). To express
#this exponent using sums of squares, we need to determine the vector
#(x - mu) and then assign vectors to each of the columns of the posterior
#precision matrix. Once this is done, by using the inprod(v1,v2) function
#we can express this exponent rather simply, using the following code:
for(i in 1:N3) #Loop for the number of observations to be classified
{
for(j in 1:4) #Loop for the number of variables in the
data set
{
x.mu.1.vec[i,j] <- Predict[i,j] - mu[1,j]
x.mu.2.vec[i,j] <- Predict[i,j] - mu[2,j]
}
for(k in 1:4) #Loop through a dummy index, necessary for
the
#"matrix multiplication" we are doing
}
temp.1.vec[i,k] <- inprod( x.mu.1.vec[i,] ,
tau[1,,k] )
temp.2.vec[i,k] <- inprod( x.mu.2.vec[i,] ,
tau[2,,k] )
}
exponent1[i] <- inprod( temp.1.vec[i,] , x.mu.1.vec[i,] )
exponent2[i] <- inprod( temp.2.vec[i,] , x.mu.2.vec[i,] )
top[i] <- pow( (2*Pi), -2 ) * pow( exp( logdet( tau[1,,]) ) , .5) * exp(-.5
* exponent1[i])
bot[i] <- pow( (2*Pi), -2 ) * pow( exp( logdet( tau[2,,]) ) , .5) * exp(-.5
* exponent2[i])
criterion[i] <- top[i]/bot[i]
}
for(i in 1:4)
{
for(j in 1:4)
{
sigma1[i,j] <- inverse( tau[1,,] , i
, j )
sigma2[i,j] <- inverse( tau[2,,] , i
, j )
}
}
mu[1, 1:4] ~ dmnorm(mean[1:4], precision[1:4, 1:4])
mu[2, 1:4] ~ dmnorm(mean[1:4], precision[1:4, 1:4])
tau[1, 1:4, 1:4] ~ dwish(R[1:4, 1:4], 4)
tau[2, 1:4, 1:4] ~ dwish(R[1:4, 1:4], 4)
}
#In the following data list, note that the original data, reproduced on page
295 of
#Seber's book (taken from Lubischew's study of flea beetle measurements,
1962), is
#in the form
#
# x1 x2 x3 x4
# 189 245 137 163
# 192 260 132 217
# . . . .
#
#and so on. We have stacked the 19 data points for the haltica oleracea on
top of
#the 20 data points for the haltica carduorum, for a grand total of 39 data
points.
#WinBUGS requires that data in this sort of structure be listed by rows.
The .DIM
#statements at the end of the data structures specify (#rows, #columns).
The
#FleaBeetles structure corresponds to the actual data, while the Predict
structure
#corresponds to two new flea beetle observations we wish to classify.
list(FleaBeetles= structure(.Data = c(189., 245.,137., 163.,192., 260.,
132., 217.,217.,
276., 141., 192.,221., 299., 142., 213.,171., 239., 128., 158.,192., 262.,
147., 173.,213.,
278., 136., 201.,192., 255., 128., 185.,170., 244., 128., 192.,201., 276.,
146., 186.,195.,
242., 128., 192.,205., 263., 147., 192.,180., 252., 121., 167.,192., 283.,
138., 183.,200.,
294., 138., 188.,192., 277., 150., 177.,200., 287., 136., 173.,181., 255.,
146., 183.,192.,
287., 141., 198.,181., 305., 184., 209.,158., 237., 133., 188.,184., 300.,
166., 231.,171.,
273., 162., 213.,181., 297., 163., 224.,181., 308., 160., 223.,177., 301.,
166., 221.,198.,
308., 141., 197.,180., 286., 146., 214.,177., 299., 171., 192.,176., 317.,
166., 213.,192.,
312., 166., 209.,176., 285., 141., 200.,169., 287., 162., 214.,164., 265.,
147., 192.,181.,
308., 157., 204.,192., 276., 154., 209.,181., 278., 149., 235.,175., 271.,
140., 192.,197.,
303., 170., 205.) , .Dim=c(39,4)),
#The reader will notice that the data values to be used for prediction
purposes
#lead to some very interesting (and sharp) posterior distributions.
Predict = structure(.Data=c(160, 280, 140, 193, 187, 278, 160, 197),
.Dim=c(2,4)),
N1=19, N2=20, N3=2, Pi=3.14159,
mean = c(186.8 , 279.2 , 147.5 , 197.9),
# We are specifying a precision matrix for the means that implies low
# precision for them, and relative independence (the means do not
# covary)
precision = structure(.Data = c(1.0E-6, 0, 0, 0, 0,1.0E-6, 0, 0, 0,
0,1.0E-6, 0, 0, 0,0,1.0E-6),
.Dim = c(4,4)),
#According to the .5 WinBUGS manual, it is best to think of R as a prior
#estimate of the magnitude of the covariance matrix; the numbers in
#this matrix were chosen because they roughly match the pooled covariance
#matrix that Seber obtains. MCMC updating will be required in any case
#since the correlations are left unspecified.
R = structure(.Data = c(140, 0, 0, 0, 0, 365, 0, 0, 0, 0, 110, 0, 0, 0, 0,
205), .Dim = c(4, 4))
)
#I used the "GenInits" option to get the initial values for this model.
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|
