Winbugs 模拟过程的探讨

暗香疏影N01

1413

收藏 2018-04-30

鄙人最近又开始建模编程了，假装自己是业内高手。哈哈

最近在使用贝叶斯统计推断模型WinBUGS模拟水生态系统中溶解氧动力学的问题，希望通过所监测的溶解氧，光照辐射等数据的变化，模拟出模型中的一些关键参数。模型拟合值比较正常，但从WinBUGS输出结果来看，有几个关键的参数仍然不能收敛，也无法做模型比较。

现将模型模拟和调试过程分享在这里。另外，对于遇到的一些Error及自己的理解，也一并贴在这里。欢迎同行讨论，纠正。

model;

{

for( i in 1 : N ) {

lambda ~ dnorm(u,tau) }

for( i in 1 : N ) {

u <- P - R + Rear

P <- mu * IZ

IZ <- I * exp(e)

e <- ( -1) * Kd

Kd <- a + b * Tur

R <- R20 * pow(1.047,T - 20)

Rear <- K * (1 - (Sat + Sat[i + 1]) / 200)

T ~ dnorm(25,1.0E-6)I(20,35)

I ~ dnorm( 0.0,1.0E-6)

Tur ~ dnorm(10, 0.1)I(0,)

Sat ~ dunif(0,200)

}

tau ~ dgamma(0.001,0.001)
R20 ~ dnorm( 0.0,1.0E-6)
K ~ dnorm( 0.0,1.0E-6)
mu ~ dnorm( 0.0,1.0E-6)
b ~ dnorm( 0.0,1.0E-6)
a ~ dnorm( 0.0,1.0E-6)
sigma <- 1 / sqrt(tau)
}

data 是另起一个对话框
list(lambda=c(0.3, 0, 0.03,       0.12,       0.13,       0.12,       0.03,       0.27,       0.29,       0.02,       0.2,       0.25,       0.26,       0.15,       0.16, 0.74,       0.3,       0.4,       0.4,       0.28,       0.15,       -0.15,       -0.07,       0.02,       -0.13,       -0.3,       -0.26,       -0.36,       -0.26,       -0.28,       -0.26,       -0.32,       -0.18,       -0.29,       -0.27,       -0.09,       -0.32,       -0.21,       -0.18,       -0.16,       -0.23,       -0.18,       -0.16,       -0.13,       -0.18,       -0.07,       -0.15,       -0.11,       -0.03,       0.01,       0.03,       0.12,       -0.07,       0.12,       0.08,       0.18,       0.24,       0.3,       0.08,       0.35,       0.27,       0.29,       0.02,       0.2,       0.25,       0.26,       0.24,       0.2,       0.2,       -0.09,       -0.2,       -0.16,       -0.16,       -0.08,       -0.06,       -0.17,       -0.32,       -0.14,       -0.18,       -0.32,       -0.16,       -0.28,       -0.04,       -0.27,       -0.15,       -0.06,       -0.15,       -0.11,       -0.15,       -0.11,       -0.03,       -0.25,       -0.02,       -0.06,       -0.16,       -0.08,       0.08),
  T=c(27.45,27.56,27.33,27.39,27.61,27.61,28,28.4,28.78,29,29.5,29.94,30,30.22,30.56,31,31.4,31.8,31.9,30.82,30.15,30.25,30.25,29.86,29.94,29.76,29.52,29.26,29.02,28.8,28.59,28.33,28.09,27.85,27.63,27.46,27.31,27.17,27.04,26.91,26.77,26.63,26.52,26.41,26.29,26.13,26.01,25.88,25.79,25.79,25.77,25.75,25.76,25.86,25.8,25.85,26.01,26.15,26.27,26.49,26.74,27,27.56,27.87,27.83,28.02,28.36,28.67,29.11,29.45,29.55,29.35,28.67,28.52,28.05,27.98,27.91,27.55,27.26,26.99,26.72,26.46,26.21,25.98,25.71,25.44,25.23,25.01,24.79,24.6,24.41,24.23,24.05,23.88,23.73,23.59,23.48 ),

I=c(24,52,141,138,179,193,148,306,210,423,306,258,235,385,157,382,598,517,427,603,279,196,235,124,132,87,42,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7,20,59,90,109,112,86,143,206,248,244,260,384,361,460,563,398,204,510,614,608,588,523,430,332,246,161,73,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,22,90),

Sat=c(79.2,78.8,78.9,79.3,81,79.9,81.5,82.8,85.2,88.4,92.4,95.6,98.3,103,107.3,117.5,120.4,124.2,128.3,130,138,143.3,143.9,143.1,142.6,140.4,135.8,131.7,126.5,122.5,118.4,114.5,109.8,107.1,103,99.2,94.9,93.4,90.5,88,85.7,82.6,80.2,78,76.2,73.8,72.8,70.6,69.2,68.8,68.9,69.3,71,69.9,71.5,72.8,75.2,78.4,82.4,85.6,88.3,93,97.3,97.5,100.4,104.2,108.3,112.3,115.5,118.4,116.5,115,110,108,106,105.4,102.5,98,95.7,92.9,88.4,86.1,82.2,81.3,77.5,75.5,74.3,72.2,70.6,68.6,67.1,66.4,63.2,62.8,62,59.9,58.8,58.8),
      N=97)
list(mu=0.5, a=1, b=2, K=0.5, R20=1, tau=1, Tur=10)