離散データで定義された表面下の体積を計算する方法は？

Question

離散データポイントで表される一連のサーフェスの下のボリュームを決定する必要があります。私のデータでは、各サンプルは、データフレームのリスト内の別個のデータフレームとして格納されます。同様の質問へ

df1 <- data.frame(x=c(2,2,2,3,3,3,4,4,4,5,5,5,6,6,6), 
        y=c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3), 
        z=c(0,2,0,4,6,7,3,2,1,2,7,8,9,4,2)) 

df2 <- data.frame(x=c(2,2,2,3,3,3,4,4,4,5,5,5,6,6,6), 
        y=c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3), 
        z=c(1,1,2,3,5,6,2,1,3,3,8,9,8,3,1)) 

DF <- list(df1,df2) 

回答がいずれかの他の言語（MATLAB、パイソン）、または解答に問題（as here）に対処するために使用可能なスクリプトが含まれていないされています。ここではいくつかの（小）データ例です。 1）データフレーム（DF）のリスト全体に適用されるRの関数として、シンプソンのルールの離散化されたバージョンを書き出します。 2）x、y、zの間の任意の関係を計算し、多変量数値積分を使って表面下の体積を求める（パッケージpracmaではsimpson2d/quad2d、キューブではadaptIntegrateのような関数を用いる）。

複合シンプソンのルール（私が使用したい）の式はhereですが、その複雑さのために、私は働いている二重合計関数を書くことに失敗しています。この式において、I（λ（em）λ（ex））は各x、y格子点の上記データセットのzに等しく、Delta（em）およびDelta（ex）はxとyの間の間隔を表す。

第2のアプローチは、基本的にアプローチfound hereを多変量スプラインフィットに拡張し、予測されたz値を積分の関数として渡します。ここで私はこれまでこのアプローチで試したことがあります：

require(pracma) 

df1.loess <- loess(z ~ x + y, data=DF[[1]]) 
mod.fun <- function(x,y) predict(df1.loess, newdata=x,y) 

simpson2d(mod.fun, x=c(2,6), y=c(1,3)) 

しかし、これは有用な結果をもたらしません。

実際には、私は個々のサンプルに対してほぼ100データフレームのリストを持っているので、リスト内のすべてのデータフレームにわたってこれらの計算を自動化する一連のラッピ関数としてソリューションを表現することが本当に必要です。例は次のようになります。

require(akima) 
DF.splines <- lapply(DF, function(x,y,z) interp(x = "x", y = "y", z = "z", 
               linear=F, nx=4, ny=2)) 

残念ながら、これは欠損値とInfsの例外を生成します。私は、これらの戦略の1つをどのようにして成功させるか、あるいは別の（より単純な）アプローチを活用するための提案については非常に公開しています。 Kriging関数（DiceKrigingパッケージのkmのような）は、数値積分のために渡すことができるより良い適合を生み出すことができますか？

Answer 1

ボリュームサーフェスメッシュは、点を直線で結んで定義されているものとします。その後、エリアA_i

ボリュームを計算する三角形T_i

のそれぞれに対応するz値Z_iを見つけると三角形T_iに

1. 経由(x,y)グリッドの三角形のテッセレーションをその表面の下にボリュームを見つけることができますV_i （T_iおよびZ_iによって定義される）を、V_i=A_i*sum(Z_i)/3（https://en.wikipedia.org/wiki/Prism_(geometry)およびhttps://math.stackexchange.com/questions/2371139/volume-of-truncated-prismを参照）
2. すべて切り捨てプリズムボリュームを合計V_i

は、ボリュームがあなたのテッセレーションによって異なり、テセレーションが一意でないことはないこと、しかし、覚えておいてください。しかし、あなたの問題はポイント間の補間方法を記述していないという意味では完全には定義されていません。したがって、ボリュームを計算するどのようなアプローチでも、追加の前提が必要になります。

私のソリューションのアプローチに戻って、ポイント1と2はgeometryパッケージで達成できます。 ここではいくつかのコード

library(geometry) 

getVolume=function(df) { 
    #find triangular tesselation of (x,y) grid 
    res=delaunayn(as.matrix(df[,-3]),full=TRUE,options="Qz") 
    #calulates sum of truncated prism volumes 
    sum(mapply(function(triPoints,A) A/3*sum(df[triPoints,"z"]), 
      split.data.frame(res$tri,seq_along(res$areas)), 
      res$areas)) 
} 

sapply(DF,getVolume) 
#[1] 32.50000 30.33333 

それが結果はここで、我々は正しい答えを知っている簡単な例一貫しているかどうかを確認するのは難しいので。これは、x軸に沿ってくさびを切り取った、辺の長さが2の立方体です。カットアウト領域は全体積の1/4です。

cutOutCube=expand.grid(c(0,1,2),c(0,1,2)) 
colnames(cutOutCube)=c("x","y") 
cutOutCube$z=ifelse(cutOutCube$x==1,1,2) 

sapply(list(cutOutCube),getVolume) 
#[1] 6 

2^3*(1-1/4)=6から正しいです。

もう1つの健全性チェックは、ボリュームw.r.tの「補数」を計算することによって実行できます。すべてzの値が最大zの値（どちらの場合でもmax(z)=9）に設定されている単純な直方体に変換します。あなたの両方のケースで、単純な立方体の容積は72です。のは、だから、ボリュームおよび補体量はいずれの場合も、右の単純な直方体の体積を与える補完面を定義し、音量を総括し、ボリューム

df1c=df1 
df1c$z=max(df1c$z)-df1c$z 
df2c=df2 
df2c$z=max(df2c$z)-df2c$z 
DFc=list(df1c,df2c) 

sapply(DFc,getVolume)+sapply(DF,getVolume) 
#[1] 72 72 

を補完することはできませ。

Answer 2

パッケージのbarylag2dファンクションに実装されているような "重心ラグランジュ"アプローチでサーフェスを近似することができます。次に、ベクトル化の問題を避けるために、ガウシアン直交ガイドを明示的に適用します。

library(pracma) 

df1 <- data.frame(x=c(2,2,2,3,3,3,4,4,4,5,5,5,6,6,6), 
        y=c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3), 
        z=c(0,2,0,4,6,7,3,2,1,2,7,8,9,4,2)) 

# Define the nodes in x- and y-direction 
xn <- df1$x[c(1,4,7,10,13)] 
yn <- df1$y[1:3] 

# Define the matrix representing the function 
m1 <- matrix(df1$z, nrow=5, byrow=TRUE) 
f <- function(x, y) 
     c(pracma::barylag2d(m1, xn, yn, x, y)) 

# 32 nodes in integration intervals 
n <- 32 
xa <- 2; xb <- 6; ya <- 1; yb <- 3 

# Apply quadrature rules explicitely 
cx <- gaussLegendre(n, xa, xb) 
x <- cx$x; wx <- cx$w 
cy <- gaussLegendre(n, ya, yb) 
y <- cy$x; wy <- cy$w 

# Sum weights * values over all nodes 
I <- 0 
for (i in 1:n) { 
    for (j in 1:n) { 
    I <- I + wx[i] * wy[j] * f(x[i], y[j]) 
    } 
} 
I # 40.37037 

データが与えられた場合、40の積分値は妥当と思われます。 simpson2dまたはquad2dはこの設定では機能しません。

adaptIntegrateがそのように定義された関数fで動作するかどうか試してみてください。