比较R中两个矩阵的图

怎么样:

par(mfrow = c(1, 3))
image(m1, zlim = c(-1, 1))
contour(m1 - m2, add = TRUE)
image(m2, zlim = c(-1, 1))
contour(m1 - m2, add = TRUE)
image(m1 - m2, zlim = c(0, 1))
contour(m1 - m2, add = TRUE)

这将在区域周围添加等高线图。有点像是在第三个图的区域周围放了一个环(可能需要摆弄等高线的(n)个级别以得到更少的“圆”)。

另一种方法是:

image(m1>m2)

我写了一些代码来做类似的事情。我想通过在阈值0.95以上的区域周围画一个框来突出显示相邻区域,所以我得到了所有阈值0.95以上的网格正方形,并对它们进行了聚类。然后对聚类输出进行一些修改,以获得区域的矩形坐标:

computeHotspots = function(xyz, thresh, minsize=1, margin=1){
### given a list(x,y,z), return a data frame where each row
### is a (xmin,xmax,ymin,ymax) of bounding box of a contiguous area
### over the given threshhold.
### or approximately. lets use the clustering tools in R...

  overs <- which(xyz$z>thresh,arr.ind=T)

  if(length(overs)==0){
    ## found no hotspots
    return(NULL)
  }

  if(length(overs)==2){
    ## found one hotspot
    xRange <- cbind(xyz$x[overs[,1]],xyz$x[overs[,1]])
    yRange <- cbind(xyz$y[overs[,2]],xyz$y[overs[,2]])
  }else{

    oTree <- hclust(dist(overs),method="single")
    oCut <- cutree(oTree,h=10)

    oXYc <- data.frame(x=xyz$x[overs[,1]],y=xyz$y[overs[,2]],oCut)

    xRange <- do.call("rbind",tapply(oXYc[,1],oCut,range))
    yRange <- do.call("rbind",tapply(oXYc[,2],oCut,range))

  }

### add user-margins
 xRange[,1] <- xRange[,1]-margin
 xRange[,2] <- xRange[,2]+margin
 yRange[,1] <- yRange[,1]-margin
 yRange[,2] <- yRange[,2]+margin

## put it all together
 xr <- apply(xRange,1,diff)
 xm <- apply(xRange,1,mean)
 xRange[xr<minsize,1] <- xm[xr<minsize]-(minsize/2)
 xRange[xr<minsize,2] <- xm[xr<minsize]+(minsize/2)

 yr <- apply(yRange,1,diff)
 ym <- apply(yRange,1,mean)
 yRange[yr<minsize,1] <- ym[yr<minsize]-(minsize/2)
 yRange[yr<minsize,2] <- ym[yr<minsize]+(minsize/2)

  cbind(xRange,yRange)

}

x=1:23
y=7:34
m1=list(x=x,y=y,z=outer(x,y,function(x,y){sin(x/3)*cos(y/3)}))
image(m1)
hs = computeHotspots(m1,0.95)

这会给你一个矩形坐标的矩阵:

> hs
  [,1] [,2] [,3] [,4]
1   13   15    8   11
2    3    6   17   20
3   22   24   18   20
4   13   16   27   30

image(m1)
rect(hs[,1],hs[,3],hs[,2],hs[,4])

image(list(x=m1$x,y=m1$y,z=m1$z>0.95))
rect(hs[,1],hs[,3],hs[,2],hs[,4])

当然,你可以修改它来画圆,但是更复杂的形状会很棘手。当感兴趣的区域相当紧凑时,它工作得最好。

巴里