#gamma addictive mixed model
model{
#define Gamm
    for(i in 1:10){
      for(j in 1:10){
        y[i,j]~dgamma(m,tau[i,j])
}
}
    for(i in 1:10){
      for(j in 1:10){
        log(tau[i,j])<-log(m)-(inprod(f[1:10],p[j,1:10])+inprod(g[1:10],x[j,1:10])+b[1]*i+b[2]*j) #link function
}
}
#define prior distribute
    for(i in 1:10){
      for(j in 1:10){
        N4[i,j]<-sigam*N3[i,j]
}
}
      sigam<-1/t
      m~dnorm(10,4)
      g[1:10]~dmnorm(u1[],N1[,])
      f[1:10]~dmnorm(u3[],N4[,])
      b[1:2]~dmnorm(u2[],N2[,])
      N2[1:2,1:2]~dwish(R[,],2)
#define invert gamma distribute with a=0.005,b=0.005
      zero<-0
      zero~dpois(phi)
      phi<-(0.005/t)*log(pow(0.005,0.005)*pow(1/t,2))
      t~dflat()
#calculate reserve crease
    for(i in 2:10){
      for(j in 10:10-i){
       y[i,j]<-exp(inprod(f[1:10],p[j,1:10])+inprod(g[1:10],x[j,1:10])+b[1]*i+b[2]*j)
}
}
}
#loading data
   list(
     y=structure(
            .Data=c(5012,3257,2638,898,1734,2642,1828,599,54,172,
                          106,4179,1111,5270,3116,1817,1245,673,535,NA,
                           3410,5528,4881,2268,2594,3479,649,603,NA,NA,
                           5655,5900,4211,5500,2159,2658,984,NA,NA,NA,
                           1092,8473,6271,6333,3786,225,NA,NA,NA,NA,
                            1513,4932,5257,1233,2917,NA,NA,NA,NA,NA,
                            557,3463,6926,1368,NA,NA,NA,NA,NA,NA,
                            1351,5596,6165,NA,NA,NA,NA,NA,NA,NA,
                             3133,2263,NA,NA,NA,NA,NA,NA,NA,NA,
                             2063,NA,NA,NA,NA,NA,NA,NA,NA,NA),
            .Dim=c(10,10)),
        N1=structure(
            .Data=c(0.3,0,0,0,0,0,0,0,0,0,
                           0,2,0,0,0,0,0,0,0,0,
                           0,0,2.5,0,0,0,0,0,0,0,
                           0,0,0,1,0,0,0,0,0,0,
                           0,0,0,0,2.4,0,0,0,0,0,
                           0,0,0,0,0,5,0,0,0,0,
                           0,0,0,0,0,0,5,0,0,0,
                           0,0,0,0,0,0,0,7,0,0,
                           0,0,0,0,0,0,0,0,6.4,0,
                           0,0,0,0,0,0,0,0,0,3),
           .Dim=c(10,10)),
         N3=structure(
            .Data=c(1,-1,0,0,0,0,0,0,0,0,
                          -1,2,-1,0,0,0,0,0,0,0,
                          0,-1,2,-1,0,0,0,0,0,0,
                          0,0,-1,2,-1,0,0,0,0,0,
                          0,0,0,1,2,-1,0,0,0,0,
                          0,0,0,0,-1,2,1,0,0,0,
                          0,0,0,0,0,-1,2,-1,0,0,
                          0,0,0,0,0,0,-1,2,-1,0,
                          0,0,0,0,0,0,0,-1,2,-1,
                          0,0,0,0,0,0,0,0,-1,1),
            .Dim=c(10,10)),
        x=structure(
             .Data=c(1,0,0,0,0,0,0,0,0,0,
                           0,1,0,0,0,0,0,0,0,0,
                           0,0,1,0,0,0,0,0,0,0,
                           0,0,0,1,0,0,0,0,0,0,
                           0,0,0,0,1,0,0,0,0,0,
                           0,0,0,0,0,1,0,0,0,0,
                           0,0,0,0,0,0,1,0,0,0,
                           0,0,0,0,0,0,0,1,0,0,
                           0,0,0,0,0,0,0,0,1,0,
                           0,0,0,0,0,0,0,0,0,1),
              .Dim=c(10,10)),
        p=structure(
                .Data=c(1,0,0,0,0,0,0,0,0,0,
                           0,1,0,0,0,0,0,0,0,0,
                           0,0,1,0,0,0,0,0,0,0,
                           0,0,0,1,0,0,0,0,0,0,
                           0,0,0,0,1,0,0,0,0,0,
                           0,0,0,0,0,1,0,0,0,0,
                           0,0,0,0,0,0,1,0,0,0,
                           0,0,0,0,0,0,0,1,0,0,
                           0,0,0,0,0,0,0,0,1,0,
                           0,0,0,0,0,0,0,0,0,1),
                .Dim=c(10,10)),
         R=structure(
                .Data=c(4,0,
                              0,2),
                .Dim=c(2,2)),
         u1=c(0,0,0,0,0,0,0,0,0,0),
         u3=c(0,0,0,0,0,0,0,0,0,0),
         u2=c(0,0))
#ini-data
       list(f=c(0,0.417,0.2,0.4,0.5,0.461,0.4,0.7,0.9,0.8),
             g=c(0,0.417,0.2,0.4,0.5,0.461,0.4,0.2,0.7,0.9),
             b=c(0.05,0.03),
             m=12,
             t=2)
      list(f=c(0.5,0.417,0.4,0.6,0.3,0.4,0.4,0.5,0.8,0.9),
            g=c(0.6,0.467,0.1,0.2,0.7,0.481,0.2,0.3,0.7,0.5),
             b=c(0.005,0.002),
             m=12,t=3)
      list(f=c(0.3,0.427,0.5,0.7,0.2,0.671,0.2,0.6,0.8,0.6),
            g=c(0.8,0.487,0.3,0.9,0.8,0.661,0.2,0.4,0.6,0.4),
             b=c(0.9,0.05),
             m=8,t=2)