私は、既知の直径と種のマッピングされた樹木のプロットで、.5m^2 の「ごみトラップ」のネットワークで収集されたごみの量を記述するモデルの適合を最適化しようとしています。選択したモデルには、ごみ生産の相対成長スケーリングと、ごみ移動距離の指数関数的減衰という 2 つの要因があります。
tree1.litter = alpha*gamma^2 * DBH^Beta/(2*pi) * exp(-gamma*z-delta*DBH)
ただし、トラップ データには複数のツリーからの入力が含まれています (これは、タイトルで言及されている「欠落レベル」です)。
Obs.Litter = tree1.litter + tree2.litter + ... + treej.litter + error
これまでのところ、シミュレートされたデータでも非常に複雑な結果が得られました。直径と距離の組み合わせが十分にあれば、関数はある程度十分に制限されるはずです。この分析は、私が模倣している記事で実行されました。ログ（Obs.Litter）の分析も試しましたが、これが道だと思います。しかし、ログ バージョンをコーディングした方法で、パフォーマンスが向上すると期待される結果が得られたかどうかはわかりません。
この時点で、このタイプの「隠れたプロセス」を使用した非線形回帰のフィッティングまたはモデルのフィッティング問題の経験が豊富な人から、あらゆる種類のアドバイス (コードベースまたは概念的) を探しているだけだと思います。データ シミュレーションのコードとさまざまな可能性を以下に示します。OpenBUGS のベイジアン階層モデルを使用してこれらのパラメーターを推定することで、有益な事前情報のみを使用して、もう少し成功しました。
library(plyr)
########################
##Generate Data#########
########################
alpha = 5
Beta = 2
gamma = .2
delta = .02
n = 600 #Number of trees
N.trap = 45 #Number of litter traps
D = rlnorm(n, 2)+5 #generate diameters
Z = runif(n, 0, 25) #generate distances
trap.id = sort(sample(1:N.trap, size = n, replace = T)) #assign trees to traps
tree.lit = (2*pi)^-1*alpha*gamma^2*D^Beta * exp(-gamma*Z-delta*D) #generate litter
log.tree.lit = -(2*pi) + log(alpha) + 2*log(gamma) + Beta*DBH -gamma*Z - delta*D
litter = data.frame(D=D, Z = Z, trap.id = trap.id, tree.lit = tree.lit)
data = ddply(litter, .(trap.id), summarize, trap.lit = sum(tree.lit), n.trees=length(trap.id) )
trap.lit = data[,2]
#####################################
##### Maximum Likelihood Optimization
#####################################
library(bbmle)
log.Litter.Func<-function(alpha, Beta, gamma, delta, sigma, D, Z, N.trap, trap.id, Obs.Litter){
log.Expected.Litter.tree = -log(2*pi) + log(alpha) + 2*log(gamma) + Beta*log(D) -gamma*Z - delta*D
log.Expected.Litter.trap = rep(0, N.trap)
for(i in 1:N.trap){
log.Expected.Litter.trap[i] <- sum(exp(log.Expected.Litter.tree[trap.id==i]))
}
-sum(dlnorm(log(Obs.Litter), log.Expected.Litter.trap, sd=sigma, log=T))
}
Litter.Func<-function(alpha, Beta, gamma, delta, sigma, D, Z, N.trap, trap.id, Obs.Litter){
Expected.Litter.tree = 1/(2*pi) * alpha * gamma^2 * D^Beta *exp(-gamma*Z - delta*D)
Expected.Litter.trap = rep(0, N.trap)
for(i in 1:N.trap){
Expected.Litter.trap[i] <- sum(Expected.Litter.tree[trap.id==i])
}
-sum(dnorm(Obs.Litter, Expected.Litter.trap, sd=sigma, log=T))
}
log.fit<-mle2(log.Litter.Func,
start = list(alpha = 5,gamma = .2,Beta = 2,delta = .02, sigma = 1),
#upper = list(alpha = 20,gamma = 1,Beta = 4,delta = .2,sigma = 20),
#lower = list(alpha = .002,gamma = .002,Beta = .0002,delta = .000000002,sigma = .020),
#method="L-BFGS-B",
data=list(D=D, Z=Z, N.trap=N.trap, trap.id=litter$trap.id, Obs.Litter=trap.lit)
)
fit<-mle2(Litter.Func,
start = list(alpha = 5,gamma = .2,Beta = 2,delta = .02, sigma = 1),
#upper = list(alpha = 20,gamma = 1,Beta = 4,delta = .2,sigma = 20),
#lower = list(alpha = .002,gamma = .002,Beta = .0002,delta = .000000002,sigma = .020),
#method="L-BFGS-B",
data=list(D = D,Z = Z,N.trap=N.trap, trap.id=litter$trap.id,Obs.Litter = trap.lit)
)