代码之家  ›  专栏  ›  技术社区  ›  zaphoid

超出数据网格的外推

  •  0
  • zaphoid  · 技术社区  · 11 年前

    我有一个(x,y,z)值的网格,我想要一个函数,当给定(x,y)点位于网格之外时,它可以近似z值。

    我已经尝试使用Akima包(代码块3)解决这个问题,但我似乎无法让interp函数与linear=FALSE选项一起工作,这是外推到网格之外所需的。

    数据:

    # Grid data
    x <- seq(0,1,length.out=6)
    y <- seq(0,1,length.out=6)
    z <- outer(x,y,function(x,y){sqrt(x^2+y^3)})
    

    可视化数据(这不是问题的关键):

    ## Visualize the data - Not important for question ##
      jet.colors <- colorRampPalette( c("Royal Blue", "Lime Green") )
      nbcol <- 100
      color <- jet.colors(nbcol)
      nrz <- nrow(z)
      ncz <- ncol(z)
      zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
      facetcol <- cut(zfacet, nbcol)
      pmat <- persp(x, y, z, d = 1,r = 4,
                    ticktype="detailed", 
                    col = color[facetcol], 
                    main = "Title",
                    xlab="x value", 
                    ylab = "y value", 
                    zlab= "z value",
                    scale=FALSE,
                    expand=0.6, 
                    theta=-40, 
                    phi=25)
    ## End visualization ##
    

    我尝试使用Akima包解决问题

    library(akima)
    
    # Vectorize the grid:
    zz <- as.vector(z)
    # create all combinations of x and y
    xy <- expand.grid(x,y)
    
    # What we want:
    sqrt(0.7^2 + 0.7^3) # c(0.7, 0.7) = 0.9126883
    sqrt(0.7^2 + 1.2^3) # c(0.7, 1.2) = 1.489295
    
    # We get a result for the first point inside the grid,
    # but not for the second one outside the grid.
    # This is expected behaviour when linear=TRUE:
    interp(xy[,1], xy[,2], zz, xo = c(0.7), yo= c(0.7, 1.2), linear=TRUE)
    # = (0.929506, NA)
    
    # When LINEAR = FALSE we get z= 0, 0!!
    interp(xy[,1], xy[,2], zz, xo = c(0.7), yo= c(0.7, 1.2), linear=FALSE, extrap = TRUE)
    # = (0, 0)
    
    # Dropping extrap=TRUE we see that both are actually NA in this case
    interp(xy[,1], xy[,2], zz, xo = c(0.7), yo= c(0.7, 1.2), linear=FALSE)
    # = (NA, NA)
    # What is going on?
    

    使用R 3.1.3和akima_0.5-11。

    1 回复  |  直到 11 年前
        1
  •  1
  •   zaphoid    11 年前

    直接使用interp.old()得到一个结果:

    interp.old(xy[,1], xy[,2], zz, xo = c(0.7), yo= c(0.7, 1.2), ncp=2, extrap=TRUE)
    # = (0.9211469, 1.472434)
    

    我不确定这是否是interp()函数中的错误。interp()是一个包装器,如果linear=TRUE,则调用interp.old(),如果linear=FALSE,则调用interp.new()。interp.new()函数似乎失败。

    请注意,不推荐使用ncp参数。根据手册:

    ncp-已弃用,请使用参数linear。现在只被interp.old()使用。意思是:在计算每个数据点的偏导数时要使用的额外点数。 ncp必须为0(不使用偏导数),或至少为2,但小于数据点的数量(且小于25)

    推荐文章