[R] pmnorm to find cdf with linear decision bound

David Kikuchi dwkikuchi at gmail.com
Thu Oct 2 03:05:12 CEST 2014


Hi all,

I need help calculating the error rate of an optimal responder to a 
multidimensional discrimination task. I have been trying to use pmnorm 
to do this, but am not sure about its functionality. I'm working with a 
dataset of responses (binary, attack/not attack) of subjects who were 
asked to discriminate between two different, overlapping categories of 
stimuli. The stimuli subjects viewed were bivariate normal in their 
distributions (squares with different mean blue:yellow ratios and 
sizes). It is easy to find a threshold (a line) on the bivariate plane 
that describes the optimal discrimination function. It is a diagonal. To 
find the error rate committed by this optimal responder, I need to find 
the cumulative distribution function for a bivariate normal that lies 
above the line. However, pmnorm only seems to calculate rectangular 
probabilities, i.e. only uses limits of integration perpendicular to the 
axes. Is there another function I could use?

Thanks,
David

The problem is visualized with the code below:

library(mnormt)
library(lattice)

modelms<- 31.2
mimicms<- 24
sds<- 4*4
modelmc<- 0.7
mimicmc<- 0.4
sdc<- 0.15*0.15
xv<-seq(0,1,0.01)
yv<-seq(10,50,0.1)

ys <- matrix(NA,length(xv),length(coeffs[1,]))
for(i in 1:length(xv)) ys[i,] <- 
(coeffs[1,]+coeffs[2,]*xv[i])/(coeffs[3,]*-1)

mu <- c(modelmc,modelms) #model
sigma <- matrix(c(sdc,0,0,sds),2,2)
z1<-NULL
for(x in xv){
   f <- dmnorm(cbind(x,yv), mu, sigma)
   z1<-rbind(z1,f)
}
contour(xv,yv,z1, nlevels = 5,col = "blue", lty = "solid", lwd = 1.8, 
xlab = "proportion yellow", ylab = "size")

mu <- c(mimicmc,mimicms) #mimic
sigma <- matrix(c(sdc,0,0,sds),2,2)
z2<-NULL
for(x in xv){
   f <- dmnorm(cbind(x,yv), mu, sigma)
   z2<-rbind(z2,f)
}
contour(xv,yv,z2, nlevels = 5,add = TRUE,col = "red", lty = "solid", lwd 
= 1.8)
contour(xv,yv,z2-z1, nlevels = 1,add = TRUE, col = "black", lty = 
"solid", lwd = 2.5)



More information about the R-help mailing list