#This program fits a hierarchical model to the wind speed data from northern/central England, as outlined in the paper "A hierarchical model for extreme wind speeds". It allows for seasonal (monthly) and site effects when a GPD is fitted to excesses over high thresholds. The approach is to update each component singly, using a Gibbs sampler where conjugacy facilitates straightforward sampling from the conditionals, and a random walk Metropolis step elsewhere.# #Section 1 - reading in the data# #The following commands read in the threshold excesses for each site# #Read in the Sheffield data# exceedshef<-scan("sheffieldexceedances.txt",skip=12) janshef.exceed<-exceedshef[1:497] febshef.exceed<-exceedshef[498:953] marshef.exceed<-exceedshef[954:1506] aprshef.exceed<-exceedshef[1507:1909] mayshef.exceed<-exceedshef[1910:2343] junshef.exceed<-exceedshef[2344:2808] julshef.exceed<-exceedshef[2809:3362] augshef.exceed<-exceedshef[3363:3839] sepshef.exceed<-exceedshef[3840:5022] octshef.exceed<-exceedshef[5023:5460] novshef.exceed<-exceedshef[5461:5927] decshef.exceed<-exceedshef[5928:6427] #Read in the Bradfield data# exceedbrad<-scan("bradfieldexceedances.txt",skip=12) janbrad.exceed<-exceedbrad[1:326] febbrad.exceed<-exceedbrad[327:659] marbrad.exceed<-exceedbrad[660:1472] aprbrad.exceed<-exceedbrad[1473:1689] maybrad.exceed<-exceedbrad[1690:2074] junbrad.exceed<-exceedbrad[2075:2729] julbrad.exceed<-exceedbrad[2730:2951] augbrad.exceed<-exceedbrad[2952:3405] sepbrad.exceed<-exceedbrad[3406:3967] octbrad.exceed<-exceedbrad[3968:4509] novbrad.exceed<-exceedbrad[4510:4863] decbrad.exceed<-exceedbrad[4864:5437] #Read in the Leeming data# exceed2245<-scan("site2245exceedances.txt",skip=12) jan2245.exceed<-exceed2245[1:844] feb2245.exceed<-exceed2245[845:1571] mar2245.exceed<-exceed2245[1572:2518] apr2245.exceed<-exceed2245[2519:3242] may2245.exceed<-exceed2245[3243:4048] jun2245.exceed<-exceed2245[4049:4922] jul2245.exceed<-exceed2245[4923:5763] aug2245.exceed<-exceed2245[5764:6567] sep2245.exceed<-exceed2245[6568:8037] oct2245.exceed<-exceed2245[8038:8737] nov2245.exceed<-exceed2245[8738:9487] dec2245.exceed<-exceed2245[9488:10317] #Read in the Linton data# exceed2255<-scan("site2255exceedances.txt",skip=12) jan2255.exceed<-exceed2255[1:991] feb2255.exceed<-exceed2255[992:1782] mar2255.exceed<-exceed2255[1783:2708] apr2255.exceed<-exceed2255[2709:3442] may2255.exceed<-exceed2255[3443:4254] jun2255.exceed<-exceed2255[4255:5007] jul2255.exceed<-exceed2255[5008:5773] aug2255.exceed<-exceed2255[5774:6512] sep2255.exceed<-exceed2255[6513:7889] oct2255.exceed<-exceed2255[7890:8770] nov2255.exceed<-exceed2255[8771:9563] dec2255.exceed<-exceed2255[9564:10439] #Read in the Scampton data# exceed2416<-scan("site2416exceedances.txt",skip=12) jan2416.exceed<-exceed2416[1:966] feb2416.exceed<-exceed2416[967:1605] mar2416.exceed<-exceed2416[1606:2378] apr2416.exceed<-exceed2416[2379:3114] may2416.exceed<-exceed2416[3115:3961] jun2416.exceed<-exceed2416[3962:4760] jul2416.exceed<-exceed2416[4761:5715] aug2416.exceed<-exceed2416[5716:6577] sep2416.exceed<-exceed2416[6578:8200] oct2416.exceed<-exceed2416[8201:9078] nov2416.exceed<-exceed2416[9079:9872] dec2416.exceed<-exceed2416[9873:10719] #Read in the Cranwell data# exceed2432<-scan("site2432exceedances.txt",skip=12) jan2432.exceed<-exceed2432[1:808] feb2432.exceed<-exceed2432[809:1556] mar2432.exceed<-exceed2432[1557:2431] apr2432.exceed<-exceed2432[2432:3300] may2432.exceed<-exceed2432[3301:4256] jun2432.exceed<-exceed2432[4257:5168] jul2432.exceed<-exceed2432[5169:5903] aug2432.exceed<-exceed2432[5904:6752] sep2432.exceed<-exceed2432[6753:8357] oct2432.exceed<-exceed2432[8358:9216] nov2432.exceed<-exceed2432[9217:10010] dec2432.exceed<-exceed2432[10011:10852] #Read in the Bingley data# exceed4036<-scan("site4036exceedances.txt",skip=12) jan4036.exceed<-exceed4036[1:861] feb4036.exceed<-exceed4036[862:1506] mar4036.exceed<-exceed4036[1507:2319] apr4036.exceed<-exceed4036[2320:3174] may4036.exceed<-exceed4036[3175:3894] jun4036.exceed<-exceed4036[3895:4695] jul4036.exceed<-exceed4036[4696:5486] aug4036.exceed<-exceed4036[5487:6292] sep4036.exceed<-exceed4036[6293:7953] oct4036.exceed<-exceed4036[7954:8629] nov4036.exceed<-exceed4036[8630:9425] dec4036.exceed<-exceed4036[9426:10271] #Read in the Nottingham data# exceed4206<-scan("site4206exceedances.txt",skip=12) jan4206.exceed<-exceed4206[1:985] feb4206.exceed<-exceed4206[986:1759] mar4206.exceed<-exceed4206[1760:2566] apr4206.exceed<-exceed4206[2567:3364] may4206.exceed<-exceed4206[3365:4231] jun4206.exceed<-exceed4206[4232:5144] jul4206.exceed<-exceed4206[5145:6010] aug4206.exceed<-exceed4206[6011:6922] sep4206.exceed<-exceed4206[6923:8483] oct4206.exceed<-exceed4206[8484:9277] nov4206.exceed<-exceed4206[9278:10124] dec4206.exceed<-exceed4206[10125:10938] #Read in the Finningley data# exceed4238<-scan("site4238exceedances.txt",skip=12) jan4238.exceed<-exceed4238[1:1126] feb4238.exceed<-exceed4238[1127:2128] mar4238.exceed<-exceed4238[2129:3080] apr4238.exceed<-exceed4238[3081:3818] may4238.exceed<-exceed4238[3819:5042] jun4238.exceed<-exceed4238[5043:5871] jul4238.exceed<-exceed4238[5872:6786] aug4238.exceed<-exceed4238[6787:7622] sep4238.exceed<-exceed4238[7623:9044] oct4238.exceed<-exceed4238[9045:9987] nov4238.exceed<-exceed4238[9988:10978] dec4238.exceed<-exceed4238[10979:11976] #Read in the Keele data# exceed4617<-scan("site4617exceedances.txt",skip=12) jan4617.exceed<-exceed4617[1:970] feb4617.exceed<-exceed4617[971:1704] mar4617.exceed<-exceed4617[1705:2505] apr4617.exceed<-exceed4617[2506:3281] may4617.exceed<-exceed4617[3282:4311] jun4617.exceed<-exceed4617[4312:5153] jul4617.exceed<-exceed4617[5154:6071] aug4617.exceed<-exceed4617[6072:6957] sep4617.exceed<-exceed4617[6958:8485] oct4617.exceed<-exceed4617[8486:9324] nov4617.exceed<-exceed4617[9325:10245] dec4617.exceed<-exceed4617[10246:11103] #Read in the Aigburth data# exceed7221<-scan("site7221exceedances.txt",skip=12) jan7221.exceed<-exceed7221[1:869] feb7221.exceed<-exceed7221[870:1721] mar7221.exceed<-exceed7221[1722:2480] apr7221.exceed<-exceed7221[2481:3222] may7221.exceed<-exceed7221[3223:4067] jun7221.exceed<-exceed7221[4068:5007] jul7221.exceed<-exceed7221[5008:5992] aug7221.exceed<-exceed7221[5993:6728] sep7221.exceed<-exceed7221[6729:8026] oct7221.exceed<-exceed7221[8027:8961] nov7221.exceed<-exceed7221[8962:9713] dec7221.exceed<-exceed7221[9714:10657] #Read in the Ringway data# exceed7377<-scan("site7377exceedances.txt",skip=12) jan7377.exceed<-exceed7377[1:1009] feb7377.exceed<-exceed7377[1010:1969] mar7377.exceed<-exceed7377[1970:2976] apr7377.exceed<-exceed7377[2977:3843] may7377.exceed<-exceed7377[3844:4717] jun7377.exceed<-exceed7377[4718:5470] jul7377.exceed<-exceed7377[5471:6511] aug7377.exceed<-exceed7377[6512:7442] sep7377.exceed<-exceed7377[7443:8943] oct7377.exceed<-exceed7377[8944:9956] nov7377.exceed<-exceed7377[9957:10864] dec7377.exceed<-exceed7377[10865:11785] #The monthly (seasonal) exceedances# dataset1<-c(jan2245.exceed,jan2255.exceed,jan2416.exceed,jan2432.exceed,jan4036.exceed,jan4206.exceed,jan4238.exceed, jan4617.exceed, jan7221.exceed,jan7377.exceed,janshef.exceed,janbrad.exceed) dataset2<-c(feb2245.exceed,feb2255.exceed,feb2416.exceed,feb2432.exceed,feb4036.exceed,feb4206.exceed,feb4238.exceed, feb4617.exceed, feb7221.exceed,feb7377.exceed,febshef.exceed,febbrad.exceed) dataset3<-c(mar2245.exceed,mar2255.exceed,mar2416.exceed,mar2432.exceed,mar4036.exceed,mar4206.exceed,mar4238.exceed, mar4617.exceed, mar7221.exceed,mar7377.exceed,marshef.exceed,marbrad.exceed) dataset4<-c(apr2245.exceed,apr2255.exceed,apr2416.exceed,apr2432.exceed,apr4036.exceed,apr4206.exceed,apr4238.exceed, apr4617.exceed, apr7221.exceed,apr7377.exceed,aprshef.exceed,aprbrad.exceed) dataset5<-c(may2245.exceed,may2255.exceed,may2416.exceed,may2432.exceed,may4036.exceed,may4206.exceed,may4238.exceed, may4617.exceed, may7221.exceed,may7377.exceed,mayshef.exceed,maybrad.exceed) dataset6<-c(jun2245.exceed,jun2255.exceed,jun2416.exceed,jun2432.exceed,jun4036.exceed,jun4206.exceed,jun4238.exceed, jun4617.exceed, jun7221.exceed,jun7377.exceed,junshef.exceed,junbrad.exceed) dataset7<-c(jul2245.exceed,jul2255.exceed,jul2416.exceed,jul2432.exceed,jul4036.exceed,jul4206.exceed,jul4238.exceed, jul4617.exceed, jul7221.exceed,jul7377.exceed,julshef.exceed,julbrad.exceed) dataset8<-c(aug2245.exceed,aug2255.exceed,aug2416.exceed,aug2432.exceed,aug4036.exceed,aug4206.exceed,aug4238.exceed, aug4617.exceed, aug7221.exceed,aug7377.exceed,augshef.exceed,augbrad.exceed) dataset9<-c(sep2245.exceed,sep2255.exceed,sep2416.exceed,sep2432.exceed,sep4036.exceed,sep4206.exceed,sep4238.exceed, sep4617.exceed, sep7221.exceed,sep7377.exceed,sepshef.exceed,sepbrad.exceed) dataset10<-c(oct2245.exceed,oct2255.exceed,oct2416.exceed,oct2432.exceed,oct4036.exceed,oct4206.exceed,oct4238.exceed, oct4617.exceed, oct7221.exceed,oct7377.exceed,octshef.exceed,octbrad.exceed) dataset11<-c(nov2245.exceed,nov2255.exceed,nov2416.exceed,nov2432.exceed,nov4036.exceed,nov4206.exceed,nov4238.exceed, nov4617.exceed, nov7221.exceed,nov7377.exceed,novshef.exceed,novbrad.exceed) dataset12<-c(dec2245.exceed,dec2255.exceed,dec2416.exceed,dec2432.exceed,dec4036.exceed,dec4206.exceed,dec4238.exceed, dec4617.exceed, dec7221.exceed,dec7377.exceed,decshef.exceed,decbrad.exceed) all.months<-c(dataset1,dataset2,dataset3,dataset4,dataset5,dataset6,dataset7,dataset8,dataset9,dataset10,dataset11,dataset12) all.sites<-c(exceed2245,exceed2255,exceed2416,exceed2432,exceed4036,exceed4206,exceed4238,exceed4617,exceed7221,exceed7377,exceedshef,exceedbrad) #Section 2 - the hierarchical model# #Now the program which fits the hierarchical model# random.effects<-function(n,seasons,sites,tau.etastart,tau.xistart,zeta.etastart,zeta.xistart,bstart,cstart,dstart,estart,fstart,gstart) { #"seasons" is the concatenated vector containing season1, season2, ..., season12, where these sub-vectors contain the excesses for each season# #"sites" is the concatenated vector containing site1, site2, ..., site12, where these sub-vectors contain the excesses for each site# #Section 2.1# #Starting values# gamma.eta1start<-1.0 gamma.eta2start<-1.0 gamma.eta3start<-1.0 gamma.eta4start<-1.0 gamma.eta5start<-1.0 gamma.eta6start<-1.0 gamma.eta7start<-1.0 gamma.eta8start<-1.0 gamma.eta9start<-1.0 gamma.eta10start<-1.0 gamma.eta11start<-1.0 gamma.eta12start<-1.0 gamma.xi1start<-0.05 gamma.xi2start<-0.05 gamma.xi3start<-0.05 gamma.xi4start<-0.05 gamma.xi5start<-0.05 gamma.xi6start<-0.05 gamma.xi7start<-0.05 gamma.xi8start<-0.05 gamma.xi9start<-0.05 gamma.xi10start<-0.05 gamma.xi11start<-0.05 gamma.xi12start<-0.05 epsilon.eta1start<-0.8 epsilon.eta2start<-0.8 epsilon.eta3start<-0.8 epsilon.eta4start<-0.8 epsilon.eta5start<-0.8 epsilon.eta6start<-0.8 epsilon.eta7start<-0.8 epsilon.eta8start<-0.8 epsilon.eta9start<-0.8 epsilon.eta10start<-0.8 epsilon.eta11start<-0.8 epsilon.eta12start<-0.8 epsilon.xi1start<-0.025 epsilon.xi2start<-0.025 epsilon.xi3start<-0.025 epsilon.xi4start<-0.025 epsilon.xi5start<-0.025 epsilon.xi6start<-0.025 epsilon.xi7start<-0.025 epsilon.xi8start<-0.025 epsilon.xi9start<-0.025 epsilon.xi10start<-0.025 epsilon.xi11start<-0.025 epsilon.xi12start<-0.025 #Section 2.2# #Variances# errgamma.eta1<-0.025 errgamma.eta2<-0.027 errgamma.eta3<-0.027 errgamma.eta4<-0.026 errgamma.eta5<-0.025 errgamma.eta6<-0.025 errgamma.eta7<-0.026 errgamma.eta8<-0.026 errgamma.eta9<-0.02 errgamma.eta10<-0.024 errgamma.eta11<-0.025 errgamma.eta12<-0.025 errgamma.xi1<-0.004 errgamma.xi2<-0.0045 errgamma.xi3<-0.0045 errgamma.xi4<-0.0045 errgamma.xi5<-0.0045 errgamma.xi6<-0.0045 errgamma.xi7<-0.0045 errgamma.xi8<-0.0045 errgamma.xi9<-0.0045 errgamma.xi10<-0.0045 errgamma.xi11<-0.0045 errgamma.xi12<-0.0045 errepsilon.eta1<-0.027 errepsilon.eta2<-0.027 errepsilon.eta3<-0.026 errepsilon.eta4<-0.026 errepsilon.eta5<-0.026 errepsilon.eta6<-0.026 errepsilon.eta7<-0.026 errepsilon.eta8<-0.027 errepsilon.eta9<-0.027 errepsilon.eta10<-0.027 errepsilon.eta11<-0.027 errepsilon.eta12<-0.027 errepsilon.xi1<-0.008 errepsilon.xi2<-0.008 errepsilon.xi3<-0.008 errepsilon.xi4<-0.008 errepsilon.xi5<-0.008 errepsilon.xi6<-0.008 errepsilon.xi7<-0.008 errepsilon.xi8<-0.008 errepsilon.xi9<-0.008 errepsilon.xi10<-0.008 errepsilon.xi11<-0.008 errepsilon.xi12<-0.008 k1<-length(dataset1) #length of excessses in January# k2<-length(dataset2) k3<-length(dataset3) k4<-length(dataset4) k5<-length(dataset5) k6<-length(dataset6) k7<-length(dataset7) k8<-length(dataset8) k9<-length(dataset9) k10<-length(dataset10) k11<-length(dataset11) k12<-length(dataset12) gamma.eta1<-gamma.eta1start gamma.eta2<-gamma.eta2start gamma.eta3<-gamma.eta3start gamma.eta4<-gamma.eta4start gamma.eta5<-gamma.eta5start gamma.eta6<-gamma.eta6start gamma.eta7<-gamma.eta7start gamma.eta8<-gamma.eta8start gamma.eta9<-gamma.eta9start gamma.eta10<-gamma.eta10start gamma.eta11<-gamma.eta11start gamma.eta12<-gamma.eta12start gamma.xi1<-gamma.xi1start gamma.xi2<-gamma.xi2start gamma.xi3<-gamma.xi3start gamma.xi4<-gamma.xi4start gamma.xi5<-gamma.xi5start gamma.xi6<-gamma.xi6start gamma.xi7<-gamma.xi7start gamma.xi8<-gamma.xi8start gamma.xi9<-gamma.xi9start gamma.xi10<-gamma.xi10start gamma.xi11<-gamma.xi11start gamma.xi12<-gamma.xi12start epsilon.eta1<-epsilon.eta1start epsilon.eta2<-epsilon.eta2start epsilon.eta3<-epsilon.eta3start epsilon.eta4<-epsilon.eta4start epsilon.eta5<-epsilon.eta5start epsilon.eta6<-epsilon.eta6start epsilon.eta7<-epsilon.eta7start epsilon.eta8<-epsilon.eta8start epsilon.eta9<-epsilon.eta9start epsilon.eta10<-epsilon.eta10start epsilon.eta11<-epsilon.eta11start epsilon.eta12<-epsilon.eta12start epsilon.xi1<-epsilon.xi1start epsilon.xi2<-epsilon.xi2start epsilon.xi3<-epsilon.xi3start epsilon.xi4<-epsilon.xi4start epsilon.xi5<-epsilon.xi5start epsilon.xi6<-epsilon.xi6start epsilon.xi7<-epsilon.xi7start epsilon.xi8<-epsilon.xi8start epsilon.xi9<-epsilon.xi9start epsilon.xi10<-epsilon.xi10start epsilon.xi11<-epsilon.xi11start epsilon.xi12<-epsilon.xi12start taueta<-tau.etastart #initialies tau<-eta^(1) at tau# tauxi<-tau.xistart #initialises tau<-xi^(1) at tau# zetaeta<-zeta.etastart #initialises zeta<-eta^(1) at zeta# zetaxi<-zeta.xistart #initialises zeta<-xi^(1) at zeta# b<-bstart c<-cstart d<-dstart e<-estart f<-fstart g<-gstart #Section 2.3# cangamma.eta1<-vector("numeric") cangamma.eta2<-vector("numeric") cangamma.eta3<-vector("numeric") cangamma.eta4<-vector("numeric") cangamma.eta5<-vector("numeric") cangamma.eta6<-vector("numeric") cangamma.eta7<-vector("numeric") cangamma.eta8<-vector("numeric") cangamma.eta9<-vector("numeric") cangamma.eta10<-vector("numeric") cangamma.eta11<-vector("numeric") cangamma.eta12<-vector("numeric") cangamma.xi1<-vector("numeric") cangamma.xi2<-vector("numeric") cangamma.xi3<-vector("numeric") cangamma.xi4<-vector("numeric") cangamma.xi5<-vector("numeric") cangamma.xi6<-vector("numeric") cangamma.xi7<-vector("numeric") cangamma.xi8<-vector("numeric") cangamma.xi9<-vector("numeric") cangamma.xi10<-vector("numeric") cangamma.xi11<-vector("numeric") cangamma.xi12<-vector("numeric") canepsilon.eta1<-vector("numeric") canepsilon.eta2<-vector("numeric") canepsilon.eta3<-vector("numeric") canepsilon.eta4<-vector("numeric") canepsilon.eta5<-vector("numeric") canepsilon.eta6<-vector("numeric") canepsilon.eta7<-vector("numeric") canepsilon.eta8<-vector("numeric") canepsilon.eta9<-vector("numeric") canepsilon.eta10<-vector("numeric") canepsilon.eta11<-vector("numeric") canepsilon.eta12<-vector("numeric") canepsilon.xi1<-vector("numeric") canepsilon.xi2<-vector("numeric") canepsilon.xi3<-vector("numeric") canepsilon.xi4<-vector("numeric") canepsilon.xi5<-vector("numeric") canepsilon.xi6<-vector("numeric") canepsilon.xi7<-vector("numeric") canepsilon.xi8<-vector("numeric") canepsilon.xi9<-vector("numeric") canepsilon.xi10<-vector("numeric") canepsilon.xi11<-vector("numeric") canepsilon.xi12<-vector("numeric") x1<-vector("numeric") #These x vectors will take the simulated values for the gamma<-eta's# x2<-vector("numeric") #i.e. the seasonal effects for eta# x3<-vector("numeric") x4<-vector("numeric") x5<-vector("numeric") x6<-vector("numeric") x7<-vector("numeric") x8<-vector("numeric") x9<-vector("numeric") x10<-vector("numeric") x11<-vector("numeric") x12<-vector("numeric") y1<-vector("numeric") #These y vectors will take the simulated values for the gamma<-xi's# y2<-vector("numeric") #i.e. the seasonal effects for xi# y3<-vector("numeric") y4<-vector("numeric") y5<-vector("numeric") y6<-vector("numeric") y7<-vector("numeric") y8<-vector("numeric") y9<-vector("numeric") y10<-vector("numeric") y11<-vector("numeric") y12<-vector("numeric") v1<-vector("numeric") #These v vectors will take the simulated values for the epsilon<-eta's# v2<-vector("numeric") #i.e. the site effects for eta# v3<-vector("numeric") v4<-vector("numeric") v5<-vector("numeric") v6<-vector("numeric") v7<-vector("numeric") v8<-vector("numeric") v9<-vector("numeric") v10<-vector("numeric") v11<-vector("numeric") v12<-vector("numeric") w1<-vector("numeric") #These w vectors will take the simulated values for the epsilon<-xi's# w2<-vector("numeric") #i.e. the site effects for xi# w3<-vector("numeric") w4<-vector("numeric") w5<-vector("numeric") w6<-vector("numeric") w7<-vector("numeric") w8<-vector("numeric") w9<-vector("numeric") w10<-vector("numeric") w11<-vector("numeric") w12<-vector("numeric") eta.one.one<-vector("numeric") eta.two.one<-vector("numeric") eta.three.one<-vector("numeric") eta.four.one<-vector("numeric") eta.five.one<-vector("numeric") eta.six.one<-vector("numeric") eta.seven.one<-vector("numeric") eta.eight.one<-vector("numeric") eta.nine.one<-vector("numeric") eta.ten.one<-vector("numeric") eta.eleven.one<-vector("numeric") eta.twelve.one<-vector("numeric") xi.one.one<-vector("numeric") xi.two.one<-vector("numeric") xi.three.one<-vector("numeric") xi.four.one<-vector("numeric") xi.five.one<-vector("numeric") xi.six.one<-vector("numeric") xi.seven.one<-vector("numeric") xi.eight.one<-vector("numeric") xi.nine.one<-vector("numeric") xi.ten.one<-vector("numeric") xi.eleven.one<-vector("numeric") xi.twelve.one<-vector("numeric") eta.one.two<-vector("numeric") eta.two.two<-vector("numeric") eta.three.two<-vector("numeric") eta.four.two<-vector("numeric") eta.five.two<-vector("numeric") eta.six.two<-vector("numeric") eta.seven.two<-vector("numeric") eta.eight.two<-vector("numeric") eta.nine.two<-vector("numeric") eta.ten.two<-vector("numeric") eta.eleven.two<-vector("numeric") eta.twelve.two<-vector("numeric") xi.one.two<-vector("numeric") xi.two.two<-vector("numeric") xi.three.two<-vector("numeric") xi.four.two<-vector("numeric") xi.five.two<-vector("numeric") xi.six.two<-vector("numeric") xi.seven.two<-vector("numeric") xi.eight.two<-vector("numeric") xi.nine.two<-vector("numeric") xi.ten.two<-vector("numeric") xi.eleven.two<-vector("numeric") xi.twelve.two<-vector("numeric") eta.one.three<-vector("numeric") eta.two.three<-vector("numeric") eta.three.three<-vector("numeric") eta.four.three<-vector("numeric") eta.five.three<-vector("numeric") eta.six.three<-vector("numeric") eta.seven.three<-vector("numeric") eta.eight.three<-vector("numeric") eta.nine.three<-vector("numeric") eta.ten.three<-vector("numeric") eta.eleven.three<-vector("numeric") eta.twelve.three<-vector("numeric") xi.one.three<-vector("numeric") xi.two.three<-vector("numeric") xi.three.three<-vector("numeric") xi.four.three<-vector("numeric") xi.five.three<-vector("numeric") xi.six.three<-vector("numeric") xi.seven.three<-vector("numeric") xi.eight.three<-vector("numeric") xi.nine.three<-vector("numeric") xi.ten.three<-vector("numeric") xi.eleven.three<-vector("numeric") xi.twelve.three<-vector("numeric") eta.one.four<-vector("numeric") eta.two.four<-vector("numeric") eta.three.four<-vector("numeric") eta.four.four<-vector("numeric") eta.five.four<-vector("numeric") eta.six.four<-vector("numeric") eta.seven.four<-vector("numeric") eta.eight.four<-vector("numeric") eta.nine.four<-vector("numeric") eta.ten.four<-vector("numeric") eta.eleven.four<-vector("numeric") eta.twelve.four<-vector("numeric") xi.one.four<-vector("numeric") xi.two.four<-vector("numeric") xi.three.four<-vector("numeric") xi.four.four<-vector("numeric") xi.five.four<-vector("numeric") xi.six.four<-vector("numeric") xi.seven.four<-vector("numeric") xi.eight.four<-vector("numeric") xi.nine.four<-vector("numeric") xi.ten.four<-vector("numeric") xi.eleven.four<-vector("numeric") xi.twelve.four<-vector("numeric") eta.one.five<-vector("numeric") eta.two.five<-vector("numeric") eta.three.five<-vector("numeric") eta.four.five<-vector("numeric") eta.five.five<-vector("numeric") eta.six.five<-vector("numeric") eta.seven.five<-vector("numeric") eta.eight.five<-vector("numeric") eta.nine.five<-vector("numeric") eta.ten.five<-vector("numeric") eta.eleven.five<-vector("numeric") eta.twelve.five<-vector("numeric") xi.one.five<-vector("numeric") xi.two.five<-vector("numeric") xi.three.five<-vector("numeric") xi.four.five<-vector("numeric") xi.five.five<-vector("numeric") xi.six.five<-vector("numeric") xi.seven.five<-vector("numeric") xi.eight.five<-vector("numeric") xi.nine.five<-vector("numeric") xi.ten.five<-vector("numeric") xi.eleven.five<-vector("numeric") xi.twelve.five<-vector("numeric") eta.one.six<-vector("numeric") eta.two.six<-vector("numeric") eta.three.six<-vector("numeric") eta.four.six<-vector("numeric") eta.five.six<-vector("numeric") eta.six.six<-vector("numeric") eta.seven.six<-vector("numeric") eta.eight.six<-vector("numeric") eta.nine.six<-vector("numeric") eta.ten.six<-vector("numeric") eta.eleven.six<-vector("numeric") eta.twelve.six<-vector("numeric") xi.one.six<-vector("numeric") xi.two.six<-vector("numeric") xi.three.six<-vector("numeric") xi.four.six<-vector("numeric") xi.five.six<-vector("numeric") xi.six.six<-vector("numeric") xi.seven.six<-vector("numeric") xi.eight.six<-vector("numeric") xi.nine.six<-vector("numeric") xi.ten.six<-vector("numeric") xi.eleven.six<-vector("numeric") xi.twelve.six<-vector("numeric") eta.one.seven<-vector("numeric") eta.two.seven<-vector("numeric") eta.three.seven<-vector("numeric") eta.four.seven<-vector("numeric") eta.five.seven<-vector("numeric") eta.six.seven<-vector("numeric") eta.seven.seven<-vector("numeric") eta.eight.seven<-vector("numeric") eta.nine.seven<-vector("numeric") eta.ten.seven<-vector("numeric") eta.eleven.seven<-vector("numeric") eta.twelve.seven<-vector("numeric") xi.one.seven<-vector("numeric") xi.two.seven<-vector("numeric") xi.three.seven<-vector("numeric") xi.four.seven<-vector("numeric") xi.five.seven<-vector("numeric") xi.six.seven<-vector("numeric") xi.seven.seven<-vector("numeric") xi.eight.seven<-vector("numeric") xi.nine.seven<-vector("numeric") xi.ten.seven<-vector("numeric") xi.eleven.seven<-vector("numeric") xi.twelve.seven<-vector("numeric") eta.one.eight<-vector("numeric") eta.two.eight<-vector("numeric") eta.three.eight<-vector("numeric") eta.four.eight<-vector("numeric") eta.five.eight<-vector("numeric") eta.six.eight<-vector("numeric") eta.seven.eight<-vector("numeric") eta.eight.eight<-vector("numeric") eta.nine.eight<-vector("numeric") eta.ten.eight<-vector("numeric") eta.eleven.eight<-vector("numeric") eta.twelve.eight<-vector("numeric") xi.one.eight<-vector("numeric") xi.two.eight<-vector("numeric") xi.three.eight<-vector("numeric") xi.four.eight<-vector("numeric") xi.five.eight<-vector("numeric") xi.six.eight<-vector("numeric") xi.seven.eight<-vector("numeric") xi.eight.eight<-vector("numeric") xi.nine.eight<-vector("numeric") xi.ten.eight<-vector("numeric") xi.eleven.eight<-vector("numeric") xi.twelve.eight<-vector("numeric") eta.one.nine<-vector("numeric") eta.two.nine<-vector("numeric") eta.three.nine<-vector("numeric") eta.four.nine<-vector("numeric") eta.five.nine<-vector("numeric") eta.six.nine<-vector("numeric") eta.seven.nine<-vector("numeric") eta.eight.nine<-vector("numeric") eta.nine.nine<-vector("numeric") eta.ten.nine<-vector("numeric") eta.eleven.nine<-vector("numeric") eta.twelve.nine<-vector("numeric") xi.one.nine<-vector("numeric") xi.two.nine<-vector("numeric") xi.three.nine<-vector("numeric") xi.four.nine<-vector("numeric") xi.five.nine<-vector("numeric") xi.six.nine<-vector("numeric") xi.seven.nine<-vector("numeric") xi.eight.nine<-vector("numeric") xi.nine.nine<-vector("numeric") xi.ten.nine<-vector("numeric") xi.eleven.nine<-vector("numeric") xi.twelve.nine<-vector("numeric") eta.one.ten<-vector("numeric") eta.two.ten<-vector("numeric") eta.three.ten<-vector("numeric") eta.four.ten<-vector("numeric") eta.five.ten<-vector("numeric") eta.six.ten<-vector("numeric") eta.seven.ten<-vector("numeric") eta.eight.ten<-vector("numeric") eta.nine.ten<-vector("numeric") eta.ten.ten<-vector("numeric") eta.eleven.ten<-vector("numeric") eta.twelve.ten<-vector("numeric") xi.one.ten<-vector("numeric") xi.two.ten<-vector("numeric") xi.three.ten<-vector("numeric") xi.four.ten<-vector("numeric") xi.five.ten<-vector("numeric") xi.six.ten<-vector("numeric") xi.seven.ten<-vector("numeric") xi.eight.ten<-vector("numeric") xi.nine.ten<-vector("numeric") xi.ten.ten<-vector("numeric") xi.eleven.ten<-vector("numeric") xi.twelve.ten<-vector("numeric") eta.one.eleven<-vector("numeric") eta.two.eleven<-vector("numeric") eta.three.eleven<-vector("numeric") eta.four.eleven<-vector("numeric") eta.five.eleven<-vector("numeric") eta.six.eleven<-vector("numeric") eta.seven.eleven<-vector("numeric") eta.eight.eleven<-vector("numeric") eta.nine.eleven<-vector("numeric") eta.ten.eleven<-vector("numeric") eta.eleven.eleven<-vector("numeric") eta.twelve.eleven<-vector("numeric") xi.one.eleven<-vector("numeric") xi.two.eleven<-vector("numeric") xi.three.eleven<-vector("numeric") xi.four.eleven<-vector("numeric") xi.five.eleven<-vector("numeric") xi.six.eleven<-vector("numeric") xi.seven.eleven<-vector("numeric") xi.eight.eleven<-vector("numeric") xi.nine.eleven<-vector("numeric") xi.ten.eleven<-vector("numeric") xi.eleven.eleven<-vector("numeric") xi.twelve.eleven<-vector("numeric") eta.one.twelve<-vector("numeric") eta.two.twelve<-vector("numeric") eta.three.twelve<-vector("numeric") eta.four.twelve<-vector("numeric") eta.five.twelve<-vector("numeric") eta.six.twelve<-vector("numeric") eta.seven.twelve<-vector("numeric") eta.eight.twelve<-vector("numeric") eta.nine.twelve<-vector("numeric") eta.ten.twelve<-vector("numeric") eta.eleven.twelve<-vector("numeric") eta.twelve.twelve<-vector("numeric") xi.one.twelve<-vector("numeric") xi.two.twelve<-vector("numeric") xi.three.twelve<-vector("numeric") xi.four.twelve<-vector("numeric") xi.five.twelve<-vector("numeric") xi.six.twelve<-vector("numeric") xi.seven.twelve<-vector("numeric") xi.eight.twelve<-vector("numeric") xi.nine.twelve<-vector("numeric") xi.ten.twelve<-vector("numeric") xi.eleven.twelve<-vector("numeric") xi.twelve.twelve<-vector("numeric") cangamma.eta1[1]<-gamma.eta1 #initialises the candiates at their starting values# cangamma.eta2[1]<-gamma.eta2 cangamma.eta3[1]<-gamma.eta3 cangamma.eta4[1]<-gamma.eta4 cangamma.eta5[1]<-gamma.eta5 cangamma.eta6[1]<-gamma.eta6 cangamma.eta7[1]<-gamma.eta7 cangamma.eta8[1]<-gamma.eta8 cangamma.eta9[1]<-gamma.eta9 cangamma.eta10[1]<-gamma.eta10 cangamma.eta11[1]<-gamma.eta11 cangamma.eta12[1]<-gamma.eta12 cangamma.xi1[1]<-gamma.xi1 cangamma.xi2[1]<-gamma.xi2 cangamma.xi3[1]<-gamma.xi3 cangamma.xi4[1]<-gamma.xi4 cangamma.xi5[1]<-gamma.xi5 cangamma.xi6[1]<-gamma.xi6 cangamma.xi7[1]<-gamma.xi7 cangamma.xi8[1]<-gamma.xi8 cangamma.xi9[1]<-gamma.xi9 cangamma.xi10[1]<-gamma.xi10 cangamma.xi11[1]<-gamma.xi11 cangamma.xi12[1]<-gamma.xi12 canepsilon.eta1[1]<-epsilon.eta1 canepsilon.eta2[1]<-epsilon.eta2 canepsilon.eta3[1]<-epsilon.eta3 canepsilon.eta4[1]<-epsilon.eta4 canepsilon.eta5[1]<-epsilon.eta5 canepsilon.eta6[1]<-epsilon.eta6 canepsilon.eta7[1]<-epsilon.eta7 canepsilon.eta8[1]<-epsilon.eta8 canepsilon.eta9[1]<-epsilon.eta9 canepsilon.eta10[1]<-epsilon.eta10 canepsilon.eta11[1]<-epsilon.eta11 canepsilon.eta12[1]<-epsilon.eta12 canepsilon.xi1[1]<-epsilon.xi1 canepsilon.xi2[1]<-epsilon.xi2 canepsilon.xi3[1]<-epsilon.xi3 canepsilon.xi4[1]<-epsilon.xi4 canepsilon.xi5[1]<-epsilon.xi5 canepsilon.xi6[1]<-epsilon.xi6 canepsilon.xi7[1]<-epsilon.xi7 canepsilon.xi8[1]<-epsilon.xi8 canepsilon.xi9[1]<-epsilon.xi9 canepsilon.xi10[1]<-epsilon.xi10 canepsilon.xi11[1]<-epsilon.xi11 canepsilon.xi12[1]<-epsilon.xi12 x1[1]<-cangamma.eta1[1] x2[1]<-cangamma.eta2[1] x3[1]<-cangamma.eta3[1] x4[1]<-cangamma.eta4[1] x5[1]<-cangamma.eta5[1] x6[1]<-cangamma.eta6[1] x7[1]<-cangamma.eta7[1] x8[1]<-cangamma.eta8[1] x9[1]<-cangamma.eta9[1] x10[1]<-cangamma.eta10[1] x11[1]<-cangamma.eta11[1] x12[1]<-cangamma.eta12[1] y1[1]<-cangamma.xi1[1] y2[1]<-cangamma.xi2[1] y3[1]<-cangamma.xi3[1] y4[1]<-cangamma.xi4[1] y5[1]<-cangamma.xi5[1] y6[1]<-cangamma.xi6[1] y7[1]<-cangamma.xi7[1] y8[1]<-cangamma.xi8[1] y9[1]<-cangamma.xi9[1] y10[1]<-cangamma.xi10[1] y11[1]<-cangamma.xi11[1] y12[1]<-cangamma.xi12[1] v1[1]<-canepsilon.eta1[1] v2[1]<-canepsilon.eta2[1] v3[1]<-canepsilon.eta3[1] v4[1]<-canepsilon.eta4[1] v5[1]<-canepsilon.eta5[1] v6[1]<-canepsilon.eta6[1] v7[1]<-canepsilon.eta7[1] v8[1]<-canepsilon.eta8[1] v9[1]<-canepsilon.eta9[1] v10[1]<-canepsilon.eta10[1] v11[1]<-canepsilon.eta11[1] v12[1]<-canepsilon.eta12[1] w1[1]<-canepsilon.xi1[1] w2[1]<-canepsilon.xi2[1] w3[1]<-canepsilon.xi3[1] w4[1]<-canepsilon.xi4[1] w5[1]<-canepsilon.xi5[1] w6[1]<-canepsilon.xi6[1] w7[1]<-canepsilon.xi7[1] w8[1]<-canepsilon.xi8[1] w9[1]<-canepsilon.xi9[1] w10[1]<-canepsilon.xi10[1] w11[1]<-canepsilon.xi11[1] w12[1]<-canepsilon.xi12[1] eta.one.one[1]<-x1[1]+v1[1] #Combines seasonal and site effects to recover full parameter estimates eta.two.one[1]<-x2[1]+v1[1] eta.three.one[1]<-x3[1]+v1[1] eta.four.one[1]<-x4[1]+v1[1] eta.five.one[1]<-x5[1]+v1[1] eta.six.one[1]<-x6[1]+v1[1] eta.seven.one[1]<-x7[1]+v1[1] eta.eight.one[1]<-x8[1]+v1[1] eta.nine.one[1]<-x9[1]+v1[1] eta.ten.one[1]<-x10[1]+v1[1] eta.eleven.one[1]<-x11[1]+v1[1] eta.twelve.one[1]<-x12[1]+v1[1] eta.one.two[1]<-x1[1]+v2[1] eta.two.two[1]<-x2[1]+v2[1] eta.three.two[1]<-x3[1]+v2[1] eta.four.two[1]<-x4[1]+v2[1] eta.five.two[1]<-x5[1]+v2[1] eta.six.two[1]<-x6[1]+v2[1] eta.seven.two[1]<-x7[1]+v2[1] eta.eight.two[1]<-x8[1]+v2[1] eta.nine.two[1]<-x9[1]+v2[1] eta.ten.two[1]<-x10[1]+v2[1] eta.eleven.two[1]<-x11[1]+v2[1] eta.twelve.two[1]<-x12[1]+v2[1] eta.one.three[1]<-x1[1]+v3[1] eta.two.three[1]<-x2[1]+v3[1] eta.three.three[1]<-x3[1]+v3[1] eta.four.three[1]<-x4[1]+v3[1] eta.five.three[1]<-x5[1]+v3[1] eta.six.three[1]<-x6[1]+v3[1] eta.seven.three[1]<-x7[1]+v3[1] eta.eight.three[1]<-x8[1]+v3[1] eta.nine.three[1]<-x9[1]+v3[1] eta.ten.three[1]<-x10[1]+v3[1] eta.eleven.three[1]<-x11[1]+v3[1] eta.twelve.three[1]<-x12[1]+v3[1] eta.one.four[1]<-x1[1]+v4[1] eta.two.four[1]<-x2[1]+v4[1] eta.three.four[1]<-x3[1]+v4[1] eta.four.four[1]<-x4[1]+v4[1] eta.five.four[1]<-x5[1]+v4[1] eta.six.four[1]<-x6[1]+v4[1] eta.seven.four[1]<-x7[1]+v4[1] eta.eight.four[1]<-x8[1]+v4[1] eta.nine.four[1]<-x9[1]+v4[1] eta.ten.four[1]<-x10[1]+v4[1] eta.eleven.four[1]<-x11[1]+v4[1] eta.twelve.four[1]<-x12[1]+v4[1] eta.one.five[1]<-x1[1]+v5[1] eta.two.five[1]<-x2[1]+v5[1] eta.three.five[1]<-x3[1]+v5[1] eta.four.five[1]<-x4[1]+v5[1] eta.five.five[1]<-x5[1]+v5[1] eta.six.five[1]<-x6[1]+v5[1] eta.seven.five[1]<-x7[1]+v5[1] eta.eight.five[1]<-x8[1]+v5[1] eta.nine.five[1]<-x9[1]+v5[1] eta.ten.five[1]<-x10[1]+v5[1] eta.eleven.five[1]<-x11[1]+v5[1] eta.twelve.five[1]<-x12[1]+v5[1] eta.one.six[1]<-x1[1]+v6[1] eta.two.six[1]<-x2[1]+v6[1] eta.three.six[1]<-x3[1]+v6[1] eta.four.six[1]<-x4[1]+v6[1] eta.five.six[1]<-x5[1]+v6[1] eta.six.six[1]<-x6[1]+v6[1] eta.seven.six[1]<-x7[1]+v6[1] eta.eight.six[1]<-x8[1]+v6[1] eta.nine.six[1]<-x9[1]+v6[1] eta.ten.six[1]<-x10[1]+v6[1] eta.eleven.six[1]<-x11[1]+v6[1] eta.twelve.six[1]<-x12[1]+v6[1] eta.one.seven[1]<-x1[1]+v7[1] eta.two.seven[1]<-x2[1]+v7[1] eta.three.seven[1]<-x3[1]+v7[1] eta.four.seven[1]<-x4[1]+v7[1] eta.five.seven[1]<-x5[1]+v7[1] eta.six.seven[1]<-x6[1]+v7[1] eta.seven.seven[1]<-x7[1]+v7[1] eta.eight.seven[1]<-x8[1]+v7[1] eta.nine.seven[1]<-x9[1]+v7[1] eta.ten.seven[1]<-x10[1]+v7[1] eta.eleven.seven[1]<-x11[1]+v7[1] eta.twelve.seven[1]<-x12[1]+v7[1] eta.one.eight[1]<-x1[1]+v8[1] eta.two.eight[1]<-x2[1]+v8[1] eta.three.eight[1]<-x3[1]+v8[1] eta.four.eight[1]<-x4[1]+v8[1] eta.five.eight[1]<-x5[1]+v8[1] eta.six.eight[1]<-x6[1]+v8[1] eta.seven.eight[1]<-x7[1]+v8[1] eta.eight.eight[1]<-x8[1]+v8[1] eta.nine.eight[1]<-x9[1]+v8[1] eta.ten.eight[1]<-x10[1]+v8[1] eta.eleven.eight[1]<-x11[1]+v8[1] eta.twelve.eight[1]<-x12[1]+v8[1] eta.one.nine[1]<-x1[1]+v9[1] eta.two.nine[1]<-x2[1]+v9[1] eta.three.nine[1]<-x3[1]+v9[1] eta.four.nine[1]<-x4[1]+v9[1] eta.five.nine[1]<-x5[1]+v9[1] eta.six.nine[1]<-x6[1]+v9[1] eta.seven.nine[1]<-x7[1]+v9[1] eta.eight.nine[1]<-x8[1]+v9[1] eta.nine.nine[1]<-x9[1]+v9[1] eta.ten.nine[1]<-x10[1]+v9[1] eta.eleven.nine[1]<-x11[1]+v9[1] eta.twelve.nine[1]<-x12[1]+v9[1] eta.one.ten[1]<-x1[1]+v10[1] eta.two.ten[1]<-x2[1]+v10[1] eta.three.ten[1]<-x3[1]+v10[1] eta.four.ten[1]<-x4[1]+v10[1] eta.five.ten[1]<-x5[1]+v10[1] eta.six.ten[1]<-x6[1]+v10[1] eta.seven.ten[1]<-x7[1]+v10[1] eta.eight.ten[1]<-x8[1]+v10[1] eta.nine.ten[1]<-x9[1]+v10[1] eta.ten.ten[1]<-x10[1]+v10[1] eta.eleven.ten[1]<-x11[1]+v10[1] eta.twelve.ten[1]<-x12[1]+v10[1] eta.one.eleven[1]<-x1[1]+v11[1] eta.two.eleven[1]<-x2[1]+v11[1] eta.three.eleven[1]<-x3[1]+v11[1] eta.four.eleven[1]<-x4[1]+v11[1] eta.five.eleven[1]<-x5[1]+v11[1] eta.six.eleven[1]<-x6[1]+v11[1] eta.seven.eleven[1]<-x7[1]+v11[1] eta.eight.eleven[1]<-x8[1]+v11[1] eta.nine.eleven[1]<-x9[1]+v11[1] eta.ten.eleven[1]<-x10[1]+v11[1] eta.eleven.eleven[1]<-x11[1]+v11[1] eta.twelve.eleven[1]<-x12[1]+v11[1] eta.one.twelve[1]<-x1[1]+v12[1] eta.two.twelve[1]<-x2[1]+v12[1] eta.three.twelve[1]<-x3[1]+v12[1] eta.four.twelve[1]<-x4[1]+v12[1] eta.five.twelve[1]<-x5[1]+v12[1] eta.six.twelve[1]<-x6[1]+v12[1] eta.seven.twelve[1]<-x7[1]+v12[1] eta.eight.twelve[1]<-x8[1]+v12[1] eta.nine.twelve[1]<-x9[1]+v12[1] eta.ten.twelve[1]<-x10[1]+v12[1] eta.eleven.twelve[1]<-x11[1]+v12[1] eta.twelve.twelve[1]<-x12[1]+v12[1] #----------------------------------------------# xi.one.one[1]<-y1[1]+w1[1] xi.two.one[1]<-y2[1]+w1[1] xi.three.one[1]<-y3[1]+w1[1] xi.four.one[1]<-y4[1]+w1[1] xi.five.one[1]<-y5[1]+w1[1] xi.six.one[1]<-y6[1]+w1[1] xi.seven.one[1]<-y7[1]+w1[1] xi.eight.one[1]<-y8[1]+w1[1] xi.nine.one[1]<-y9[1]+w1[1] xi.ten.one[1]<-y10[1]+w1[1] xi.eleven.one[1]<-y11[1]+w1[1] xi.twelve.one[1]<-y12[1]+w1[1] xi.one.two[1]<-y1[1]+w2[1] xi.two.two[1]<-y2[1]+w2[1] xi.three.two[1]<-y3[1]+w2[1] xi.four.two[1]<-y4[1]+w2[1] xi.five.two[1]<-y5[1]+w2[1] xi.six.two[1]<-y6[1]+w2[1] xi.seven.two[1]<-y7[1]+w2[1] xi.eight.two[1]<-y8[1]+w2[1] xi.nine.two[1]<-y9[1]+w2[1] xi.ten.two[1]<-y10[1]+w2[1] xi.eleven.two[1]<-y11[1]+w2[1] xi.twelve.two[1]<-y12[1]+w2[1] xi.one.three[1]<-y1[1]+w3[1] xi.two.three[1]<-y2[1]+w3[1] xi.three.three[1]<-y3[1]+w3[1] xi.four.three[1]<-y4[1]+w3[1] xi.five.three[1]<-y5[1]+w3[1] xi.six.three[1]<-y6[1]+w3[1] xi.seven.three[1]<-y7[1]+w3[1] xi.eight.three[1]<-y8[1]+w3[1] xi.nine.three[1]<-y9[1]+w3[1] xi.ten.three[1]<-y10[1]+w3[1] xi.eleven.three[1]<-y11[1]+w3[1] xi.twelve.three[1]<-y12[1]+w3[1] xi.one.four[1]<-y1[1]+w4[1] xi.two.four[1]<-y2[1]+w4[1] xi.three.four[1]<-y3[1]+w4[1] xi.four.four[1]<-y4[1]+w4[1] xi.five.four[1]<-y5[1]+w4[1] xi.six.four[1]<-y6[1]+w4[1] xi.seven.four[1]<-y7[1]+w4[1] xi.eight.four[1]<-y8[1]+w4[1] xi.nine.four[1]<-y9[1]+w4[1] xi.ten.four[1]<-y10[1]+w4[1] xi.eleven.four[1]<-y11[1]+w4[1] xi.twelve.four[1]<-y12[1]+w4[1] xi.one.five[1]<-y1[1]+w5[1] xi.two.five[1]<-y2[1]+w5[1] xi.three.five[1]<-y3[1]+w5[1] xi.four.five[1]<-y4[1]+w5[1] xi.five.five[1]<-y5[1]+w5[1] xi.six.five[1]<-y6[1]+w5[1] xi.seven.five[1]<-y7[1]+w5[1] xi.eight.five[1]<-y8[1]+w5[1] xi.nine.five[1]<-y9[1]+w5[1] xi.ten.five[1]<-y10[1]+w5[1] xi.eleven.five[1]<-y11[1]+w5[1] xi.twelve.five[1]<-y12[1]+w5[1] xi.one.six[1]<-y1[1]+w6[1] xi.two.six[1]<-y2[1]+w6[1] xi.three.six[1]<-y3[1]+w6[1] xi.four.six[1]<-y4[1]+w6[1] xi.five.six[1]<-y5[1]+w6[1] xi.six.six[1]<-y6[1]+w6[1] xi.seven.six[1]<-y7[1]+w6[1] xi.eight.six[1]<-y8[1]+w6[1] xi.nine.six[1]<-y9[1]+w6[1] xi.ten.six[1]<-y10[1]+w6[1] xi.eleven.six[1]<-y11[1]+w6[1] xi.twelve.six[1]<-y12[1]+w6[1] xi.one.seven[1]<-y1[1]+w7[1] xi.two.seven[1]<-y2[1]+w7[1] xi.three.seven[1]<-y3[1]+w7[1] xi.four.seven[1]<-y4[1]+w7[1] xi.five.seven[1]<-y5[1]+w7[1] xi.six.seven[1]<-y6[1]+w7[1] xi.seven.seven[1]<-y7[1]+w7[1] xi.eight.seven[1]<-y8[1]+w7[1] xi.nine.seven[1]<-y9[1]+w7[1] xi.ten.seven[1]<-y10[1]+w7[1] xi.eleven.seven[1]<-y11[1]+w7[1] xi.twelve.seven[1]<-y12[1]+w7[1] xi.one.eight[1]<-y1[1]+w8[1] xi.two.eight[1]<-y2[1]+w8[1] xi.three.eight[1]<-y3[1]+w8[1] xi.four.eight[1]<-y4[1]+w8[1] xi.five.eight[1]<-y5[1]+w8[1] xi.six.eight[1]<-y6[1]+w8[1] xi.seven.eight[1]<-y7[1]+w8[1] xi.eight.eight[1]<-y8[1]+w8[1] xi.nine.eight[1]<-y9[1]+w8[1] xi.ten.eight[1]<-y10[1]+w8[1] xi.eleven.eight[1]<-y11[1]+w8[1] xi.twelve.eight[1]<-y12[1]+w8[1] xi.one.nine[1]<-y1[1]+w9[1] xi.two.nine[1]<-y2[1]+w9[1] xi.three.nine[1]<-y3[1]+w9[1] xi.four.nine[1]<-y4[1]+w9[1] xi.five.nine[1]<-y5[1]+w9[1] xi.six.nine[1]<-y6[1]+w9[1] xi.seven.nine[1]<-y7[1]+w9[1] xi.eight.nine[1]<-y8[1]+w9[1] xi.nine.nine[1]<-y9[1]+w9[1] xi.ten.nine[1]<-y10[1]+w9[1] xi.eleven.nine[1]<-y11[1]+w9[1] xi.twelve.nine[1]<-y12[1]+w9[1] xi.one.ten[1]<-y1[1]+w10[1] xi.two.ten[1]<-y2[1]+w10[1] xi.three.ten[1]<-y3[1]+w10[1] xi.four.ten[1]<-y4[1]+w10[1] xi.five.ten[1]<-y5[1]+w10[1] xi.six.ten[1]<-y6[1]+w10[1] xi.seven.ten[1]<-y7[1]+w10[1] xi.eight.ten[1]<-y8[1]+w10[1] xi.nine.ten[1]<-y9[1]+w10[1] xi.ten.ten[1]<-y10[1]+w10[1] xi.eleven.ten[1]<-y11[1]+w10[1] xi.twelve.ten[1]<-y12[1]+w10[1] xi.one.eleven[1]<-y1[1]+w11[1] xi.two.eleven[1]<-y2[1]+w11[1] xi.three.eleven[1]<-y3[1]+w11[1] xi.four.eleven[1]<-y4[1]+w11[1] xi.five.eleven[1]<-y5[1]+w11[1] xi.six.eleven[1]<-y6[1]+w11[1] xi.seven.eleven[1]<-y7[1]+w11[1] xi.eight.eleven[1]<-y8[1]+w11[1] xi.nine.eleven[1]<-y9[1]+w11[1] xi.ten.eleven[1]<-y10[1]+w11[1] xi.eleven.eleven[1]<-y11[1]+w11[1] xi.twelve.eleven[1]<-y12[1]+w11[1] xi.one.twelve[1]<-y1[1]+w12[1] xi.two.twelve[1]<-y2[1]+w12[1] xi.three.twelve[1]<-y3[1]+w12[1] xi.four.twelve[1]<-y4[1]+w12[1] xi.five.twelve[1]<-y5[1]+w12[1] xi.six.twelve[1]<-y6[1]+w12[1] xi.seven.twelve[1]<-y7[1]+w12[1] xi.eight.twelve[1]<-y8[1]+w12[1] xi.nine.twelve[1]<-y9[1]+w12[1] xi.ten.twelve[1]<-y10[1]+w12[1] xi.eleven.twelve[1]<-y11[1]+w12[1] eta.twelve.twelve[1]<-x12[1]+w12[1] #--------------------------------------------------------# aprobgamma.eta1<-vector("numeric") #vectors of acceptance probabilites for parameters in the MCMC scheme# aprobgamma.eta2<-vector("numeric") aprobgamma.eta3<-vector("numeric") aprobgamma.eta4<-vector("numeric") aprobgamma.eta5<-vector("numeric") aprobgamma.eta6<-vector("numeric") aprobgamma.eta7<-vector("numeric") aprobgamma.eta8<-vector("numeric") aprobgamma.eta9<-vector("numeric") aprobgamma.eta10<-vector("numeric") aprobgamma.eta11<-vector("numeric") aprobgamma.eta12<-vector("numeric") aprobgamma.xi1<-vector("numeric") aprobgamma.xi2<-vector("numeric") aprobgamma.xi3<-vector("numeric") aprobgamma.xi4<-vector("numeric") aprobgamma.xi5<-vector("numeric") aprobgamma.xi6<-vector("numeric") aprobgamma.xi7<-vector("numeric") aprobgamma.xi8<-vector("numeric") aprobgamma.xi9<-vector("numeric") aprobgamma.xi10<-vector("numeric") aprobgamma.xi11<-vector("numeric") aprobgamma.xi12<-vector("numeric") aprobepsilon.eta1<-vector("numeric") aprobepsilon.eta2<-vector("numeric") aprobepsilon.eta3<-vector("numeric") aprobepsilon.eta4<-vector("numeric") aprobepsilon.eta5<-vector("numeric") aprobepsilon.eta6<-vector("numeric") aprobepsilon.eta7<-vector("numeric") aprobepsilon.eta8<-vector("numeric") aprobepsilon.eta9<-vector("numeric") aprobepsilon.eta10<-vector("numeric") aprobepsilon.eta11<-vector("numeric") aprobepsilon.eta12<-vector("numeric") aprobepsilon.xi1<-vector("numeric") aprobepsilon.xi2<-vector("numeric") aprobepsilon.xi3<-vector("numeric") aprobepsilon.xi4<-vector("numeric") aprobepsilon.xi5<-vector("numeric") aprobepsilon.xi6<-vector("numeric") aprobepsilon.xi7<-vector("numeric") aprobepsilon.xi8<-vector("numeric") aprobepsilon.xi9<-vector("numeric") aprobepsilon.xi10<-vector("numeric") aprobepsilon.xi11<-vector("numeric") aprobepsilon.xi12<-vector("numeric") tau.eta1<-taueta #initialies tau<-eta^(1) at tau# tau.xi1<-tauxi #initialises tau<-xi^(1) at tau# zeta.eta1<-zetaeta #initialises zeta<-eta^(1) at zeta# zeta.xi1<-zetaxi #initialises zeta<-xi^(1) at zeta# #Section 2.4# loglik<-function(GAMMAETA,GAMMAXI,EPSILONETA,EPSILONXI,DATA,U) { nn<-min(1+(((GAMMAXI+EPSILONXI)*DATA)/((exp(GAMMAETA+EPSILONETA))+((GAMMAXI+EPSILONXI)*U)))) if(((exp(GAMMAETA+EPSILONETA))+(GAMMAXI+EPSILONXI)*U)<1e-05)return(as.double(-1e+06)) if(nn<1e-05)return(as.double(-1e+06)) llh<- (-length(DATA)*log(exp(GAMMAETA+EPSILONETA)+(GAMMAXI+EPSILONXI)*U))-(1+(1/(GAMMAXI+EPSILONXI)))*sum(log(1+(((GAMMAXI+EPSILONXI)*DATA)/(exp(GAMMAETA+EPSILONETA)+((GAMMAXI+EPSILONXI)*U))))) } thresh<-c(35.36897,32.97613,32.97613,25.58329,25.3809,24.58329,24.98806,25.58329,25.58329, 31.3809,32.17852,33.57135,59.777,45.65,51.269,43.396,35.395,37.888,35.142,35.147,35.128,48.396,52.396,53.781,32.826,29.355,30.619,29.413,24.148,24.413,23.677,25.413,25.413,28.619,29.884,30.090,31.914,31.123,31.331,27.748,24.748,23.748,22.957,23.748,23.748,28.123,29.748,30.540,35.458,32.458,33.629,28.629,25.801,24.972,25.143,25.972,25.972,30.629,32.801,33.629,38.541,35.033,35.033,28.033,25.525,25.279,25.525,26.525,26.525,33.033,35.033,37.541,32.914,31.123,31.748,27.748,25.540,23.957,23.957,23.748,23.748,28.331,29.540,31.331,31.612,28.612,30.228,27.152,22.844,23.152,23.152,24.152,24.152,27.228,28.536,30.228,31.612,28.612,30.228,27.152,22.844,23.152,23.152,24.152,24.152,27.228,28.536,30.228,33.603,33.083,32.083,28.322,24.322,24.802,23.802,23.802,23.802,28.083,29.322,32.083,36.100,30.910,33.290,28.480,25.480,24.670,24.860,25.860,25.860,28.910,32.290,32.910,31.619,29.619,29.884,26.148,24.413,23.677,21.413,22.413,22.413,26.884,27.619,30.884) #These thresholds have been obtained for each season, within each site (starting with January, Sheffield), using a mean residual life plot# #Section 2.5# #Now the iterations begin!# for(i in 2:n){ a.eta1<-rnorm(1,(c*b+(zeta.eta1*(v1[i-1]+v2[i-1]+v3[i-1]+v4[i-1]+v5[i-1]+v6[i-1]+v7[i-1]+v8[i-1]+v9[i-1]+v10[i-1]+v11[i-1]+v12[i-1])))/(c+12*zeta.eta1),1/(sqrt(c+12*zeta.eta1))) a.xi1<-rnorm(1,(c*b+(zeta.xi1*(w1[i-1]+w2[i-1]+w3[i-1]+w4[i-1]+w5[i-1]+w6[i-1]+w7[i-1]+w8[i-1]+w9[i-1]+w10[i-1]+w11[i-1]+w12[i-1])))/(c+12*zeta.xi1),1/(sqrt(c+12*zeta.xi1))) zeta.eta2<-rgamma(1,f+6,1/(g+0.5*(((v1[i-1]-a.eta1)^2)+((v2[i-1]-a.eta1)^2)+((v3[i-1]-a.eta1)^2)+((v4[i-1]-a.eta1)^2)+((v5[i-1]-a.eta1)^2)+((v6[i-1]-a.eta1)^2)+((v7[i-1]-a.eta1)^2)+((v8[i-1]-a.eta1)^2)+((v9[i-1]-a.eta1)^2)+((v10[i-1]-a.eta1)^2)+((v11[i-1]-a.eta1)^2)+((v12[i-1]-a.eta1)^2)))) zeta.xi2<-rgamma(1,f+6,1/(g+0.5*(((w1[i-1]-a.xi1)^2)+((w2[i-1]-a.xi1)^2)+((w3[i-1]-a.xi1)^2)+((w4[i-1]-a.xi1)^2)+((w5[i-1]-a.xi1)^2)+((w6[i-1]-a.xi1)^2)+((w7[i-1]-a.xi1)^2)+((w8[i-1]-a.xi1)^2)+((w9[i-1]-a.xi1)^2)+((w10[i-1]-a.xi1)^2)+((w11[i-1]-a.xi1)^2)+((w12[i-1]-a.xi1)^2)))) tau.eta2<-rgamma(1,d+6,1/(e+0.5*(((x1[i-1])^2)+((x2[i-1])^2)+((x3[i-1])^2)+((x4[i-1])^2)+((x5[i-1])^2)+((x6[i-1])^2)+((x7[i-1])^2)+((x8[i-1])^2)+((x9[i-1])^2)+((x10[i-1])^2)+((x11[i-1])^2)+((x12[i-1])^2)))) tau.xi2<-rgamma(1,d+6,1/(e+0.5*(((y1[i-1])^2)+((y2[i-1])^2)+((y3[i-1])^2)+((y4[i-1])^2)+((y5[i-1])^2)+((y6[i-1])^2)+((y7[i-1])^2)+((y8[i-1])^2)+((y9[i-1])^2)+((y10[i-1])^2)+((y11[i-1])^2)+((y12[i-1])^2)))) cangamma.eta1[i]<-x1[i-1]+rnorm(1,0,errgamma.eta1) likely1<-exp(loglik(cangamma.eta1[i],y1[i-1],v1[i-1],w1[i-1],janshef.exceed,thresh[1]) -loglik(x1[i-1],y1[i-1],v1[i-1],w1[i-1],janshef.exceed,thresh[1]) +loglik(cangamma.eta1[i],y1[i-1],v2[i-1],w2[i-1],janbrad.exceed,thresh[13]) -loglik(x1[i-1],y1[i-1],v2[i-1],w2[i-1],janbrad.exceed,thresh[13]) +loglik(cangamma.eta1[i],y1[i-1],v3[i-1],w3[i-1],jan2245.exceed,thresh[25]) -loglik(x1[i-1],y1[i-1],v3[i-1],w3[i-1],jan2245.exceed,thresh[25]) +loglik(cangamma.eta1[i],y1[i-1],v4[i-1],w4[i-1],jan2255.exceed,thresh[37]) -loglik(x1[i-1],y1[i-1],v4[i-1],w4[i-1],jan2255.exceed,thresh[37]) +loglik(cangamma.eta1[i],y1[i-1],v5[i-1],w5[i-1],jan2416.exceed,thresh[49]) -loglik(x1[i-1],y1[i-1],v5[i-1],w5[i-1],jan2416.exceed,thresh[49]) +loglik(cangamma.eta1[i],y1[i-1],v6[i-1],w6[i-1],jan2432.exceed,thresh[61]) -loglik(x1[i-1],y1[i-1],v6[i-1],w6[i-1],jan2432.exceed,thresh[61]) +loglik(cangamma.eta1[i],y1[i-1],v7[i-1],w7[i-1],jan4036.exceed,thresh[73]) -loglik(x1[i-1],y1[i-1],v7[i-1],w7[i-1],jan4036.exceed,thresh[73]) +loglik(cangamma.eta1[i],y1[i-1],v8[i-1],w8[i-1],jan4206.exceed,thresh[85]) -loglik(x1[i-1],y1[i-1],v8[i-1],w8[i-1],jan4206.exceed,thresh[85]) +loglik(cangamma.eta1[i],y1[i-1],v9[i-1],w9[i-1],jan4238.exceed,thresh[97]) -loglik(x1[i-1],y1[i-1],v9[i-1],w9[i-1],jan4238.exceed,thresh[97]) +loglik(cangamma.eta1[i],y1[i-1],v10[i-1],w10[i-1],jan4617.exceed,thresh[109]) -loglik(x1[i-1],y1[i-1],v10[i-1],w10[i-1],jan4617.exceed,thresh[109]) +loglik(cangamma.eta1[i],y1[i-1],v11[i-1],w11[i-1],jan7221.exceed,thresh[121]) -loglik(x1[i-1],y1[i-1],v11[i-1],w11[i-1],jan7221.exceed,thresh[121]) +loglik(cangamma.eta1[i],y1[i-1],v12[i-1],w12[i-1],jan7377.exceed,thresh[133]) -loglik(x1[i-1],y1[i-1],v12[i-1],w12[i-1],jan7377.exceed,thresh[133])) aprobgamma.eta1[i]<-min(1,(likely1*((dnorm(cangamma.eta1[i],0,1/(sqrt(tau.eta2)))/(dnorm(x1[i-1],0,1/(sqrt(tau.eta2)))))))) u<-runif(1) if(u=aprobgamma.eta1[i]){ x1[i]<-x1[i-1] } cangamma.xi1[i]<-y1[i-1]+rnorm(1,0,errgamma.xi1) alikely1<-exp(loglik(x1[i],cangamma.xi1[i],v1[i-1],w1[i-1],janshef.exceed,thresh[1]) -loglik(x1[i],y1[i-1],v1[i-1],w1[i-1],janshef.exceed,thresh[1]) +loglik(x1[i],cangamma.xi1[i],v2[i-1],w2[i-1],janbrad.exceed,thresh[13]) -loglik(x1[i],y1[i-1],v2[i-1],w2[i-1],janbrad.exceed,thresh[13]) +loglik(x1[i],cangamma.xi1[i],v3[i-1],w3[i-1],jan2245.exceed,thresh[25]) -loglik(x1[i],y1[i-1],v3[i-1],w3[i-1],jan2245.exceed,thresh[25]) +loglik(x1[i],cangamma.xi1[i],v4[i-1],w4[i-1],jan2255.exceed,thresh[37]) -loglik(x1[i],y1[i-1],v4[i-1],w4[i-1],jan2255.exceed,thresh[37]) +loglik(x1[i],cangamma.xi1[i],v5[i-1],w5[i-1],jan2416.exceed,thresh[49]) -loglik(x1[i],y1[i-1],v5[i-1],w5[i-1],jan2416.exceed,thresh[49]) +loglik(x1[i],cangamma.xi1[i],v6[i-1],w6[i-1],jan2432.exceed,thresh[61]) -loglik(x1[i],y1[i-1],v6[i-1],w6[i-1],jan2432.exceed,thresh[61]) +loglik(x1[i],cangamma.xi1[i],v7[i-1],w7[i-1],jan4036.exceed,thresh[73]) -loglik(x1[i],y1[i-1],v7[i-1],w7[i-1],jan4036.exceed,thresh[73]) +loglik(x1[i],cangamma.xi1[i],v8[i-1],w8[i-1],jan4206.exceed,thresh[85]) -loglik(x1[i],y1[i-1],v8[i-1],w8[i-1],jan4206.exceed,thresh[85]) +loglik(x1[i],cangamma.xi1[i],v9[i-1],w9[i-1],jan4238.exceed,thresh[97]) -loglik(x1[i],y1[i-1],v9[i-1],w9[i-1],jan4238.exceed,thresh[97]) +loglik(x1[i],cangamma.xi1[i],v10[i-1],w10[i-1],jan4617.exceed,thresh[109]) -loglik(x1[i],y1[i-1],v10[i-1],w10[i-1],jan4617.exceed,thresh[109]) +loglik(x1[i],cangamma.xi1[i],v11[i-1],w11[i-1],jan7221.exceed,thresh[121]) -loglik(x1[i],y1[i-1],v11[i-1],w11[i-1],jan7221.exceed,thresh[121]) +loglik(x1[i],cangamma.xi1[i],v12[i-1],w12[i-1],jan7377.exceed,thresh[133]) -loglik(x1[i],y1[i-1],v12[i-1],w12[i-1],jan7377.exceed,thresh[133])) aprobgamma.xi1[i]<-min(1,(alikely1*((dnorm(cangamma.xi1[i],0,1/(sqrt(tau.xi2)))/(dnorm(y1[i-1],0,1/(sqrt(tau.xi2)))))))) if(u=aprobgamma.xi1[i]){ y1[i]<-y1[i-1] } canepsilon.eta1[i]<-v1[i-1]+rnorm(1,0,errepsilon.eta1) blikely1<-exp(loglik(x1[i],y1[i],canepsilon.eta1[i],w1[i-1],janshef.exceed,thresh[1]) -loglik(x1[i],y1[i],v1[i-1],w1[i-1],janshef.exceed,thresh[1]) +loglik(x2[i-1],y2[i-1],canepsilon.eta1[i],w1[i-1],febshef.exceed,thresh[2]) -loglik(x2[i-1],y2[i-1],v1[i-1],w1[i-1],febshef.exceed,thresh[2]) +loglik(x3[i-1],y3[i-1],canepsilon.eta1[i],w1[i-1],marshef.exceed,thresh[3]) -loglik(x3[i-1],y3[i-1],v1[i-1],w1[i-1],marshef.exceed,thresh[3]) +loglik(x4[i-1],y4[i-1],canepsilon.eta1[i],w1[i-1],aprshef.exceed,thresh[4]) -loglik(x4[i-1],y4[i-1],v1[i-1],w1[i-1],aprshef.exceed,thresh[4]) +loglik(x5[i-1],y5[i-1],canepsilon.eta1[i],w1[i-1],mayshef.exceed,thresh[5]) -loglik(x5[i-1],y5[i-1],v1[i-1],w1[i-1],mayshef.exceed,thresh[5]) +loglik(x6[i-1],y6[i-1],canepsilon.eta1[i],w1[i-1],junshef.exceed,thresh[6]) -loglik(x6[i-1],y6[i-1],v1[i-1],w1[i-1],junshef.exceed,thresh[6]) +loglik(x7[i-1],y7[i-1],canepsilon.eta1[i],w1[i-1],julshef.exceed,thresh[7]) -loglik(x7[i-1],y7[i-1],v1[i-1],w1[i-1],julshef.exceed,thresh[7]) +loglik(x8[i-1],y8[i-1],canepsilon.eta1[i],w1[i-1],augshef.exceed,thresh[8]) -loglik(x8[i-1],y8[i-1],v1[i-1],w1[i-1],augshef.exceed,thresh[8]) +loglik(x9[i-1],y9[i-1],canepsilon.eta1[i],w1[i-1],sepshef.exceed,thresh[9]) -loglik(x9[i-1],y9[i-1],v1[i-1],w1[i-1],sepshef.exceed,thresh[9]) +loglik(x10[i-1],y10[i-1],canepsilon.eta1[i],w1[i-1],octshef.exceed,thresh[10]) -loglik(x10[i-1],y10[i-1],v1[i-1],w1[i-1],octshef.exceed,thresh[10]) +loglik(x11[i-1],y11[i-1],canepsilon.eta1[i],w1[i-1],novshef.exceed,thresh[11]) -loglik(x11[i-1],y11[i-1],v1[i-1],w1[i-1],novshef.exceed,thresh[11]) +loglik(x12[i-1],y12[i-1],canepsilon.eta1[i],w1[i-1],decshef.exceed,thresh[12]) -loglik(x12[i-1],y12[i-1],v1[i-1],w1[i-1],decshef.exceed,thresh[12])) aprobepsilon.eta1[i]<-min(1,(blikely1*((dnorm(canepsilon.eta1[i],a.eta1,1/(sqrt(zeta.eta2)))/(dnorm(v1[i-1],a.eta1,1/(sqrt(zeta.eta2)))))))) if(u=aprobepsilon.eta1[i]){ v1[i]<-v1[i-1] } canepsilon.xi1[i]<-w1[i-1]+rnorm(1,0,errepsilon.xi1) clikely1<-exp(loglik(x1[i],y1[i],v1[i],canepsilon.xi1[i],janshef.exceed,thresh[1]) -loglik(x1[i],y1[i],v1[i],w1[i-1],janshef.exceed,thresh[1]) +loglik(x2[i-1],y2[i-1],v1[i],canepsilon.xi1[i],febshef.exceed,thresh[2]) -loglik(x2[i-1],y2[i-1],v1[i],w1[i-1],febshef.exceed,thresh[2]) +loglik(x3[i-1],y3[i-1],v1[i],canepsilon.xi1[i],marshef.exceed,thresh[3]) -loglik(x3[i-1],y3[i-1],v1[i],w1[i-1],marshef.exceed,thresh[3]) +loglik(x4[i-1],y4[i-1],v1[i],canepsilon.xi1[i],aprshef.exceed,thresh[4]) -loglik(x4[i-1],y4[i-1],v1[i],w1[i-1],aprshef.exceed,thresh[4]) +loglik(x5[i-1],y5[i-1],v1[i],canepsilon.xi1[i],mayshef.exceed,thresh[5]) -loglik(x5[i-1],y5[i-1],v1[i],w1[i-1],mayshef.exceed,thresh[5]) +loglik(x6[i-1],y6[i-1],v1[i],canepsilon.xi1[i],junshef.exceed,thresh[6]) -loglik(x6[i-1],y6[i-1],v1[i],w1[i-1],junshef.exceed,thresh[6]) +loglik(x7[i-1],y7[i-1],v1[i],canepsilon.xi1[i],julshef.exceed,thresh[7]) -loglik(x7[i-1],y7[i-1],v1[i],w1[i-1],julshef.exceed,thresh[7]) +loglik(x8[i-1],y8[i-1],v1[i],canepsilon.xi1[i],augshef.exceed,thresh[8]) -loglik(x8[i-1],y8[i-1],v1[i],w1[i-1],augshef.exceed,thresh[8]) +loglik(x9[i-1],y9[i-1],v1[i],canepsilon.xi1[i],sepshef.exceed,thresh[9]) -loglik(x9[i-1],y9[i-1],v1[i],w1[i-1],sepshef.exceed,thresh[9]) +loglik(x10[i-1],y10[i-1],v1[i],canepsilon.xi1[i],octshef.exceed,thresh[10]) -loglik(x10[i-1],y10[i-1],v1[i],w1[i-1],octshef.exceed,thresh[10]) +loglik(x11[i-1],y11[i-1],v1[i],canepsilon.xi1[i],novshef.exceed,thresh[11]) -loglik(x11[i-1],y11[i-1],v1[i],w1[i-1],novshef.exceed,thresh[11]) +loglik(x12[i-1],y12[i-1],v1[i],canepsilon.xi1[i],decshef.exceed,thresh[12]) -loglik(x12[i-1],y12[i-1],v1[i],w1[i-1],decshef.exceed,thresh[12])) aprobepsilon.xi1[i]<-min(1,(clikely1*((dnorm(canepsilon.xi1[i],a.xi1,1/(sqrt(zeta.xi2)))/(dnorm(w1[i-1],a.xi1/(sqrt(zeta.xi2)))))))) if(u=aprobepsilon.xi1[i]){ w1[i]<-w1[i-1] } #eta.one.one[i]<-x1[i]+v1[i] #xi.one.one[i]<-y1[i]+w1[i] #***FIRST UPDATE DONE!***# a.eta2<-rnorm(1,(c*b+(zeta.eta2*(v1[i]+v2[i-1]+v3[i-1]+v4[i-1]+v5[i-1]+v6[i-1]+v7[i-1]+v8[i-1]+v9[i-1]+v10[i-1]+v11[i-1]+v12[i-1])))/(c+12*zeta.eta2),1/(sqrt(c+12*zeta.eta2))) a.xi2<-rnorm(1,(c*b+(zeta.xi2*(w1[i]+w2[i-1]+w3[i-1]+w4[i-1]+w5[i-1]+w6[i-1]+w7[i-1]+w8[i-1]+w9[i-1]+w10[i-1]+w11[i-1]+w12[i-1])))/(c+12*zeta.xi2),1/(sqrt(c+12*zeta.xi2))) zeta.eta3<-rgamma(1,f+6,1/(g+0.5*(((v1[i]-a.eta2)^2)+((v2[i-1]-a.eta2)^2)+((v3[i-1]-a.eta2)^2)+((v4[i-1]-a.eta2)^2)+((v5[i-1]-a.eta2)^2)+((v6[i-1]-a.eta2)^2)+((v7[i-1]-a.eta2)^2)+((v8[i-1]-a.eta2)^2)+((v9[i-1]-a.eta2)^2)+((v10[i-1]-a.eta2)^2)+((v11[i-1]-a.eta2)^2)+((v12[i-1]-a.eta2)^2)))) zeta.xi3<-rgamma(1,f+6,1/(g+0.5*(((w1[i]-a.xi2)^2)+((w2[i-1]-a.xi2)^2)+((w3[i-1]-a.xi2)^2)+((w4[i-1]-a.xi2)^2)+((w5[i-1]-a.xi2)^2)+((w6[i-1]-a.xi2)^2)+((w7[i-1]-a.xi2)^2)+((w8[i-1]-a.xi2)^2)+((w9[i-1]-a.xi2)^2)+((w10[i-1]-a.xi2)^2)+((w11[i-1]-a.xi2)^2)+((w12[i-1]-a.xi2)^2)))) tau.eta3<-rgamma(1,d+6,1/(e+0.5*(((x1[i])^2)+((x2[i-1])^2)+((x3[i-1])^2)+((x4[i-1])^2)+((x5[i-1])^2)+((x6[i-1])^2)+((x7[i-1])^2)+((x8[i-1])^2)+((x9[i-1])^2)+((x10[i-1])^2)+((x11[i-1])^2)+((x12[i-1])^2)))) tau.xi3<-rgamma(1,d+6,1/(e+0.5*(((y1[i])^2)+((y2[i-1])^2)+((y3[i-1])^2)+((y4[i-1])^2)+((y5[i-1])^2)+((y6[i-1])^2)+((y7[i-1])^2)+((y8[i-1])^2)+((y9[i-1])^2)+((y10[i-1])^2)+((y11[i-1])^2)+((y12[i-1])^2)))) cangamma.eta2[i]<-x2[i-1]+rnorm(1,0,errgamma.eta2) likely2<-exp(loglik(cangamma.eta2[i],y2[i-1],v1[i],w1[i],febshef.exceed,thresh[2]) -loglik(x2[i-1],y2[i-1],v1[i],w1[i],febshef.exceed,thresh[2]) +loglik(cangamma.eta2[i],y2[i-1],v2[i-1],w2[i-1],febbrad.exceed,thresh[14]) -loglik(x2[i-1],y2[i-1],v2[i-1],w2[i-1],febbrad.exceed,thresh[14]) +loglik(cangamma.eta2[i],y2[i-1],v3[i-1],w3[i-1],feb2245.exceed,thresh[26]) -loglik(x2[i-1],y2[i-1],v3[i-1],w3[i-1],feb2245.exceed,thresh[26]) +loglik(cangamma.eta2[i],y2[i-1],v4[i-1],w4[i-1],feb2255.exceed,thresh[38]) -loglik(x2[i-1],y2[i-1],v4[i-1],w4[i-1],feb2255.exceed,thresh[38]) +loglik(cangamma.eta2[i],y2[i-1],v5[i-1],w5[i-1],feb2416.exceed,thresh[50]) -loglik(x2[i-1],y2[i-1],v5[i-1],w5[i-1],feb2416.exceed,thresh[50]) +loglik(cangamma.eta2[i],y2[i-1],v6[i-1],w6[i-1],feb2432.exceed,thresh[62]) -loglik(x2[i-1],y2[i-1],v6[i-1],w6[i-1],feb2432.exceed,thresh[62]) +loglik(cangamma.eta2[i],y2[i-1],v7[i-1],w7[i-1],feb4036.exceed,thresh[74]) -loglik(x2[i-1],y2[i-1],v7[i-1],w7[i-1],feb4036.exceed,thresh[74]) +loglik(cangamma.eta2[i],y2[i-1],v8[i-1],w8[i-1],feb4206.exceed,thresh[86]) -loglik(x2[i-1],y2[i-1],v8[i-1],w8[i-1],feb4206.exceed,thresh[86]) +loglik(cangamma.eta2[i],y2[i-1],v9[i-1],w9[i-1],feb4238.exceed,thresh[98]) -loglik(x2[i-1],y2[i-1],v9[i-1],w9[i-1],feb4238.exceed,thresh[98]) +loglik(cangamma.eta2[i],y2[i-1],v10[i-1],w10[i-1],feb4617.exceed,thresh[110]) -loglik(x2[i-1],y2[i-1],v10[i-1],w10[i-1],feb4617.exceed,thresh[110]) +loglik(cangamma.eta2[i],y2[i-1],v11[i-1],w11[i-1],feb7221.exceed,thresh[122]) -loglik(x2[i-1],y2[i-1],v11[i-1],w11[i-1],feb7221.exceed,thresh[122]) +loglik(cangamma.eta2[i],y2[i-1],v12[i-1],w12[i-1],feb7377.exceed,thresh[134]) -loglik(x2[i-1],y2[i-1],v12[i-1],w12[i-1],feb7377.exceed,thresh[134])) aprobgamma.eta2[i]<-min(1,(likely2*((dnorm(cangamma.eta2[i],0,1/(sqrt(tau.eta3)))/(dnorm(x2[i-1],0,1/(sqrt(tau.eta3)))))))) u<-runif(1) if(u=aprobgamma.eta2[i]){ x2[i]<-x2[i-1] } cangamma.xi2[i]<-y2[i-1]+rnorm(1,0,errgamma.xi2) alikely2<-exp(loglik(x2[i],cangamma.xi2[i],v1[i],w1[i],febshef.exceed,thresh[2]) -loglik(x2[i],y2[i-1],v1[i],w1[i],febshef.exceed,thresh[2]) +loglik(x2[i],cangamma.xi2[i],v2[i-1],w2[i-1],febbrad.exceed,thresh[14]) -loglik(x2[i],y2[i-1],v2[i-1],w2[i-1],febbrad.exceed,thresh[14]) +loglik(x2[i],cangamma.xi2[i],v3[i-1],w3[i-1],feb2245.exceed,thresh[26]) -loglik(x2[i],y2[i-1],v3[i-1],w3[i-1],feb2245.exceed,thresh[26]) +loglik(x2[i],cangamma.xi2[i],v4[i-1],w4[i-1],feb2255.exceed,thresh[38]) -loglik(x2[i],y2[i-1],v4[i-1],w4[i-1],feb2255.exceed,thresh[38]) +loglik(x2[i],cangamma.xi2[i],v5[i-1],w5[i-1],feb2416.exceed,thresh[50]) -loglik(x2[i],y2[i-1],v5[i-1],w5[i-1],feb2416.exceed,thresh[50]) +loglik(x2[i],cangamma.xi2[i],v6[i-1],w6[i-1],feb2432.exceed,thresh[62]) -loglik(x2[i],y2[i-1],v6[i-1],w6[i-1],feb2432.exceed,thresh[62]) +loglik(x2[i],cangamma.xi2[i],v7[i-1],w7[i-1],feb4036.exceed,thresh[74]) -loglik(x2[i],y2[i-1],v7[i-1],w7[i-1],feb4036.exceed,thresh[74]) +loglik(x2[i],cangamma.xi2[i],v8[i-1],w8[i-1],feb4206.exceed,thresh[86]) -loglik(x2[i],y2[i-1],v8[i-1],w8[i-1],feb4206.exceed,thresh[86]) +loglik(x2[i],cangamma.xi2[i],v9[i-1],w9[i-1],feb4238.exceed,thresh[98]) -loglik(x2[i],y2[i-1],v9[i-1],w9[i-1],feb4238.exceed,thresh[98]) +loglik(x2[i],cangamma.xi2[i],v10[i-1],w10[i-1],feb4617.exceed,thresh[110]) -loglik(x2[i],y2[i-1],v10[i-1],w10[i-1],feb4617.exceed,thresh[110]) +loglik(x2[i],cangamma.xi2[i],v11[i-1],w11[i-1],feb7221.exceed,thresh[122]) -loglik(x2[i],y2[i-1],v11[i-1],w11[i-1],feb7221.exceed,thresh[122]) +loglik(x2[i],cangamma.xi2[i],v12[i-1],w12[i-1],feb7377.exceed,thresh[134]) -loglik(x2[i],y2[i-1],v12[i-1],w12[i-1],feb7377.exceed,thresh[134])) aprobgamma.xi2[i]<-min(1,(alikely2*((dnorm(cangamma.xi2[i],0,1/(sqrt(tau.xi3)))/(dnorm(y2[i-1],0,1/(sqrt(tau.xi3)))))))) if(u=aprobgamma.xi2[i]){ y2[i]<-y2[i-1] } canepsilon.eta2[i]<-v2[i-1]+rnorm(1,0,errepsilon.eta2) blikely2<-exp(loglik(x1[i],y1[i],canepsilon.eta2[i],w2[i-1],janbrad.exceed,thresh[13]) -loglik(x1[i],y1[i],v2[i-1],w2[i-1],janbrad.exceed,thresh[13]) +loglik(x2[i],y2[i],canepsilon.eta2[i],w2[i-1],febbrad.exceed,thresh[14]) -loglik(x2[i],y2[i],v2[i-1],w2[i-1],febbrad.exceed,thresh[14]) +loglik(x3[i-1],y3[i-1],canepsilon.eta2[i],w2[i-1],marbrad.exceed,thresh[15]) -loglik(x3[i-1],y3[i-1],v2[i-1],w2[i-1],marbrad.exceed,thresh[15]) +loglik(x4[i-1],y4[i-1],canepsilon.eta2[i],w2[i-1],aprbrad.exceed,thresh[16]) -loglik(x4[i-1],y4[i-1],v2[i-1],w2[i-1],aprbrad.exceed,thresh[16]) +loglik(x5[i-1],y5[i-1],canepsilon.eta2[i],w2[i-1],maybrad.exceed,thresh[17]) -loglik(x5[i-1],y5[i-1],v2[i-1],w2[i-1],maybrad.exceed,thresh[17]) +loglik(x6[i-1],y6[i-1],canepsilon.eta2[i],w2[i-1],junbrad.exceed,thresh[18]) -loglik(x6[i-1],y6[i-1],v2[i-1],w2[i-1],junbrad.exceed,thresh[18]) +loglik(x7[i-1],y7[i-1],canepsilon.eta2[i],w2[i-1],julbrad.exceed,thresh[19]) -loglik(x7[i-1],y7[i-1],v2[i-1],w2[i-1],julbrad.exceed,thresh[19]) +loglik(x8[i-1],y8[i-1],canepsilon.eta2[i],w2[i-1],augbrad.exceed,thresh[20]) -loglik(x8[i-1],y8[i-1],v2[i-1],w2[i-1],augbrad.exceed,thresh[20]) +loglik(x9[i-1],y9[i-1],canepsilon.eta2[i],w2[i-1],sepbrad.exceed,thresh[21]) -loglik(x9[i-1],y9[i-1],v2[i-1],w2[i-1],sepbrad.exceed,thresh[21]) +loglik(x10[i-1],y10[i-1],canepsilon.eta2[i],w2[i-1],octbrad.exceed,thresh[22]) -loglik(x10[i-1],y10[i-1],v2[i-1],w2[i-1],octbrad.exceed,thresh[22]) +loglik(x11[i-1],y11[i-1],canepsilon.eta2[i],w2[i-1],novbrad.exceed,thresh[23]) -loglik(x11[i-1],y11[i-1],v2[i-1],w2[i-1],novbrad.exceed,thresh[23]) +loglik(x12[i-1],y12[i-1],canepsilon.eta2[i],w2[i-1],decbrad.exceed,thresh[24]) -loglik(x12[i-1],y12[i-1],v2[i-1],w2[i-1],decbrad.exceed,thresh[24])) aprobepsilon.eta2[i]<-min(1,(blikely2*((dnorm(canepsilon.eta2[i],a.eta2,1/(sqrt(zeta.eta3)))/(dnorm(v2[i-1],a.eta2,1/(sqrt(zeta.eta3)))))))) if(u=aprobepsilon.eta2[i]){ v2[i]<-v2[i-1] } canepsilon.xi2[i]<-w2[i-1]+rnorm(1,0,errepsilon.xi2) clikely2<-exp(loglik(x1[i],y1[i],v2[i],canepsilon.xi2[i],janbrad.exceed,thresh[13]) -loglik(x1[i],y1[i],v2[i],w2[i-1],janbrad.exceed,thresh[13]) +loglik(x2[i],y2[i],v2[i],canepsilon.xi2[i],febbrad.exceed,thresh[14]) -loglik(x2[i],y2[i],v2[i],w2[i-1],febbrad.exceed,thresh[14]) +loglik(x3[i-1],y3[i-1],v2[i],canepsilon.xi2[i],marbrad.exceed,thresh[15]) -loglik(x3[i-1],y3[i-1],v2[i],w2[i-1],marbrad.exceed,thresh[15]) +loglik(x4[i-1],y4[i-1],v2[i],canepsilon.xi2[i],aprbrad.exceed,thresh[16]) -loglik(x4[i-1],y4[i-1],v2[i],w2[i-1],aprbrad.exceed,thresh[16]) +loglik(x5[i-1],y5[i-1],v2[i],canepsilon.xi2[i],maybrad.exceed,thresh[17]) -loglik(x5[i-1],y5[i-1],v2[i],w2[i-1],maybrad.exceed,thresh[17]) +loglik(x6[i-1],y6[i-1],v2[i],canepsilon.xi2[i],junbrad.exceed,thresh[18]) -loglik(x6[i-1],y6[i-1],v2[i],w2[i-1],junbrad.exceed,thresh[18]) +loglik(x7[i-1],y7[i-1],v2[i],canepsilon.xi2[i],julbrad.exceed,thresh[19]) -loglik(x7[i-1],y7[i-1],v2[i],w2[i-1],julbrad.exceed,thresh[19]) +loglik(x8[i-1],y8[i-1],v2[i],canepsilon.xi2[i],augbrad.exceed,thresh[20]) -loglik(x8[i-1],y8[i-1],v2[i],w2[i-1],augbrad.exceed,thresh[20]) +loglik(x9[i-1],y9[i-1],v2[i],canepsilon.xi2[i],sepbrad.exceed,thresh[21]) -loglik(x9[i-1],y9[i-1],v2[i],w2[i-1],sepbrad.exceed,thresh[21]) +loglik(x10[i-1],y10[i-1],v2[i],canepsilon.xi2[i],octbrad.exceed,thresh[22]) -loglik(x10[i-1],y10[i-1],v2[i],w2[i-1],octbrad.exceed,thresh[22]) +loglik(x11[i-1],y11[i-1],v2[i],canepsilon.xi2[i],novbrad.exceed,thresh[23]) -loglik(x11[i-1],y11[i-1],v2[i],w2[i-1],novbrad.exceed,thresh[23]) +loglik(x12[i-1],y12[i-1],v2[i],canepsilon.xi2[i],decbrad.exceed,thresh[24]) -loglik(x12[i-1],y12[i-1],v2[i],w2[i-1],decbrad.exceed,thresh[24])) aprobepsilon.xi2[i]<-min(1,(clikely2*((dnorm(canepsilon.xi2[i],a.xi2,1/(sqrt(zeta.xi3)))/(dnorm(w2[i-1],a.xi2/(sqrt(zeta.xi3)))))))) if(u=aprobepsilon.xi2[i]){ w2[i]<-w2[i-1] } #eta.one.one[i]<-x1[i]+v1[i] #xi.one.one[i]<-y1[i]+w1[i] #***SECOND UPDATE DONE!***# a.eta3<-rnorm(1,(c*b+(zeta.eta3*(v1[i]+v2[i]+v3[i-1]+v4[i-1]+v5[i-1]+v6[i-1]+v7[i-1]+v8[i-1]+v9[i-1]+v10[i-1]+v11[i-1]+v12[i-1])))/(c+12*zeta.eta3),1/(sqrt(c+12*zeta.eta3))) a.xi3<-rnorm(1,(c*b+(zeta.xi3*(w1[i]+w2[i]+w3[i-1]+w4[i-1]+w5[i-1]+w6[i-1]+w7[i-1]+w8[i-1]+w9[i-1]+w10[i-1]+w11[i-1]+w12[i-1])))/(c+12*zeta.xi3),1/(sqrt(c+12*zeta.xi3))) zeta.eta4<-rgamma(1,f+6,1/(g+0.5*(((v1[i]-a.eta3)^2)+((v2[i]-a.eta3)^2)+((v3[i-1]-a.eta3)^2)+((v4[i-1]-a.eta3)^2)+((v5[i-1]-a.eta3)^2)+((v6[i-1]-a.eta3)^2)+((v7[i-1]-a.eta3)^2)+((v8[i-1]-a.eta3)^2)+((v9[i-1]-a.eta3)^2)+((v10[i-1]-a.eta3)^2)+((v11[i-1]-a.eta3)^2)+((v12[i-1]-a.eta3)^2)))) zeta.xi4<-rgamma(1,f+6,1/(g+0.5*(((w1[i]-a.xi3)^2)+((w2[i]-a.xi3)^2)+((w3[i-1]-a.xi3)^2)+((w4[i-1]-a.xi3)^2)+((w5[i-1]-a.xi3)^2)+((w6[i-1]-a.xi3)^2)+((w7[i-1]-a.xi3)^2)+((w8[i-1]-a.xi3)^2)+((w9[i-1]-a.xi3)^2)+((w10[i-1]-a.xi3)^2)+((w11[i-1]-a.xi3)^2)+((w12[i-1]-a.xi3)^2)))) tau.eta4<-rgamma(1,d+6,1/(e+0.5*(((x1[i])^2)+((x2[i])^2)+((x3[i-1])^2)+((x4[i-1])^2)+((x5[i-1])^2)+((x6[i-1])^2)+((x7[i-1])^2)+((x8[i-1])^2)+((x9[i-1])^2)+((x10[i-1])^2)+((x11[i-1])^2)+((x12[i-1])^2)))) tau.xi4<-rgamma(1,d+6,1/(e+0.5*(((y1[i])^2)+((y2[i])^2)+((y3[i-1])^2)+((y4[i-1])^2)+((y5[i-1])^2)+((y6[i-1])^2)+((y7[i-1])^2)+((y8[i-1])^2)+((y9[i-1])^2)+((y10[i-1])^2)+((y11[i-1])^2)+((y12[i-1])^2)))) cangamma.eta3[i]<-x3[i-1]+rnorm(1,0,errgamma.eta3) likely3<-exp(loglik(cangamma.eta3[i],y3[i-1],v1[i],w1[i],marshef.exceed,thresh[3]) -loglik(x3[i-1],y3[i-1],v1[i],w1[i],marshef.exceed,thresh[3]) +loglik(cangamma.eta3[i],y3[i-1],v2[i],w2[i],marbrad.exceed,thresh[15]) -loglik(x3[i-1],y3[i-1],v2[i],w2[i],marbrad.exceed,thresh[15]) +loglik(cangamma.eta3[i],y3[i-1],v3[i-1],w3[i-1],mar2245.exceed,thresh[27]) -loglik(x3[i-1],y3[i-1],v3[i-1],w3[i-1],mar2245.exceed,thresh[27]) +loglik(cangamma.eta3[i],y3[i-1],v4[i-1],w4[i-1],mar2255.exceed,thresh[39]) -loglik(x3[i-1],y3[i-1],v4[i-1],w4[i-1],mar2255.exceed,thresh[39]) +loglik(cangamma.eta3[i],y3[i-1],v5[i-1],w5[i-1],mar2416.exceed,thresh[51]) -loglik(x3[i-1],y3[i-1],v5[i-1],w5[i-1],mar2416.exceed,thresh[51]) +loglik(cangamma.eta3[i],y3[i-1],v6[i-1],w6[i-1],mar2432.exceed,thresh[63]) -loglik(x3[i-1],y3[i-1],v6[i-1],w6[i-1],mar2432.exceed,thresh[63]) +loglik(cangamma.eta3[i],y3[i-1],v7[i-1],w7[i-1],mar4036.exceed,thresh[75]) -loglik(x3[i-1],y3[i-1],v7[i-1],w7[i-1],mar4036.exceed,thresh[75]) +loglik(cangamma.eta3[i],y3[i-1],v8[i-1],w8[i-1],mar4206.exceed,thresh[87]) -loglik(x3[i-1],y3[i-1],v8[i-1],w8[i-1],mar4206.exceed,thresh[87]) +loglik(cangamma.eta3[i],y3[i-1],v9[i-1],w9[i-1],mar4238.exceed,thresh[99]) -loglik(x3[i-1],y3[i-1],v9[i-1],w9[i-1],mar4238.exceed,thresh[99]) +loglik(cangamma.eta3[i],y3[i-1],v10[i-1],w10[i-1],mar4617.exceed,thresh[111]) -loglik(x3[i-1],y3[i-1],v10[i-1],w10[i-1],mar4617.exceed,thresh[111]) +loglik(cangamma.eta3[i],y3[i-1],v11[i-1],w11[i-1],mar7221.exceed,thresh[123]) -loglik(x3[i-1],y3[i-1],v11[i-1],w11[i-1],mar7221.exceed,thresh[123]) +loglik(cangamma.eta3[i],y3[i-1],v12[i-1],w12[i-1],mar7377.exceed,thresh[135]) -loglik(x3[i-1],y3[i-1],v12[i-1],w12[i-1],mar7377.exceed,thresh[135])) aprobgamma.eta3[i]<-min(1,(likely3*((dnorm(cangamma.eta3[i],0,1/(sqrt(tau.eta4)))/(dnorm(x3[i-1],0,1/(sqrt(tau.eta4)))))))) u<-runif(1) if(u=aprobgamma.eta3[i]){ x3[i]<-x3[i-1] } cangamma.xi3[i]<-y3[i-1]+rnorm(1,0,errgamma.xi3) alikely3<-exp(loglik(x3[i],cangamma.xi3[i],v1[i],w1[i],marshef.exceed,thresh[3]) -loglik(x3[i],y3[i-1],v1[i],w1[i],marshef.exceed,thresh[3]) +loglik(x3[i],cangamma.xi3[i],v2[i],w2[i],marbrad.exceed,thresh[15]) -loglik(x3[i],y3[i-1],v2[i],w2[i],marbrad.exceed,thresh[15]) +loglik(x3[i],cangamma.xi3[i],v3[i-1],w3[i-1],mar2245.exceed,thresh[27]) -loglik(x3[i],y3[i-1],v3[i-1],w3[i-1],mar2245.exceed,thresh[27]) +loglik(x3[i],cangamma.xi3[i],v4[i-1],w4[i-1],mar2255.exceed,thresh[39]) -loglik(x3[i],y3[i-1],v4[i-1],w4[i-1],mar2255.exceed,thresh[39]) +loglik(x3[i],cangamma.xi3[i],v5[i-1],w5[i-1],mar2416.exceed,thresh[51]) -loglik(x3[i],y3[i-1],v5[i-1],w5[i-1],mar2416.exceed,thresh[51]) +loglik(x3[i],cangamma.xi3[i],v6[i-1],w6[i-1],mar2432.exceed,thresh[63]) -loglik(x3[i],y3[i-1],v6[i-1],w6[i-1],mar2432.exceed,thresh[63]) +loglik(x3[i],cangamma.xi3[i],v7[i-1],w7[i-1],mar4036.exceed,thresh[75]) -loglik(x3[i],y3[i-1],v7[i-1],w7[i-1],mar4036.exceed,thresh[75]) +loglik(x3[i],cangamma.xi3[i],v8[i-1],w8[i-1],mar4206.exceed,thresh[87]) -loglik(x3[i],y3[i-1],v8[i-1],w8[i-1],mar4206.exceed,thresh[87]) +loglik(x3[i],cangamma.xi3[i],v9[i-1],w9[i-1],mar4238.exceed,thresh[99]) -loglik(x3[i],y3[i-1],v9[i-1],w9[i-1],mar4238.exceed,thresh[99]) +loglik(x3[i],cangamma.xi3[i],v10[i-1],w10[i-1],mar4617.exceed,thresh[111]) -loglik(x3[i],y3[i-1],v10[i-1],w10[i-1],mar4617.exceed,thresh[111]) +loglik(x3[i],cangamma.xi3[i],v11[i-1],w11[i-1],mar7221.exceed,thresh[123]) -loglik(x3[i],y3[i-1],v11[i-1],w11[i-1],mar7221.exceed,thresh[123]) +loglik(x3[i],cangamma.xi3[i],v12[i-1],w12[i-1],mar7377.exceed,thresh[135]) -loglik(x3[i],y3[i-1],v12[i-1],w12[i-1],mar7377.exceed,thresh[135])) aprobgamma.xi3[i]<-min(1,(alikely3*((dnorm(cangamma.xi3[i],0,1/(sqrt(tau.xi4)))/(dnorm(y3[i-1],0,1/(sqrt(tau.xi4)))))))) if(u=aprobgamma.xi3[i]){ y3[i]<-y3[i-1] } canepsilon.eta3[i]<-v3[i-1]+rnorm(1,0,errepsilon.eta3) blikely3<-exp(loglik(x1[i],y1[i],canepsilon.eta3[i],w3[i-1],jan2245.exceed,thresh[25]) -loglik(x1[i],y1[i],v3[i-1],w3[i-1],jan2245.exceed,thresh[25]) +loglik(x2[i],y2[i],canepsilon.eta3[i],w3[i-1],feb2245.exceed,thresh[26]) -loglik(x2[i],y2[i],v3[i-1],w3[i-1],feb2245.exceed,thresh[26]) +loglik(x3[i],y3[i],canepsilon.eta3[i],w3[i-1],mar2245.exceed,thresh[27]) -loglik(x3[i],y3[i],v3[i-1],w3[i-1],mar2245.exceed,thresh[27]) +loglik(x4[i-1],y4[i-1],canepsilon.eta3[i],w3[i-1],apr2245.exceed,thresh[28]) -loglik(x4[i-1],y4[i-1],v3[i-1],w3[i-1],apr2245.exceed,thresh[28]) +loglik(x5[i-1],y5[i-1],canepsilon.eta3[i],w3[i-1],may2245.exceed,thresh[29]) -loglik(x5[i-1],y5[i-1],v3[i-1],w3[i-1],may2245.exceed,thresh[29]) +loglik(x6[i-1],y6[i-1],canepsilon.eta3[i],w3[i-1],jun2245.exceed,thresh[30]) -loglik(x6[i-1],y6[i-1],v3[i-1],w3[i-1],jun2245.exceed,thresh[30]) +loglik(x7[i-1],y7[i-1],canepsilon.eta3[i],w3[i-1],jul2245.exceed,thresh[31]) -loglik(x7[i-1],y7[i-1],v3[i-1],w3[i-1],jul2245.exceed,thresh[31]) +loglik(x8[i-1],y8[i-1],canepsilon.eta3[i],w3[i-1],aug2245.exceed,thresh[32]) -loglik(x8[i-1],y8[i-1],v3[i-1],w3[i-1],aug2245.exceed,thresh[32]) +loglik(x9[i-1],y9[i-1],canepsilon.eta3[i],w3[i-1],sep2245.exceed,thresh[33]) -loglik(x9[i-1],y9[i-1],v3[i-1],w3[i-1],sep2245.exceed,thresh[33]) +loglik(x10[i-1],y10[i-1],canepsilon.eta3[i],w3[i-1],oct2245.exceed,thresh[34]) -loglik(x10[i-1],y10[i-1],v3[i-1],w3[i-1],oct2245.exceed,thresh[34]) +loglik(x11[i-1],y11[i-1],canepsilon.eta3[i],w3[i-1],nov2245.exceed,thresh[35]) -loglik(x11[i-1],y11[i-1],v3[i-1],w3[i-1],nov2245.exceed,thresh[35]) +loglik(x12[i-1],y12[i-1],canepsilon.eta3[i],w3[i-1],dec2245.exceed,thresh[36]) -loglik(x12[i-1],y12[i-1],v3[i-1],w3[i-1],dec2245.exceed,thresh[36])) aprobepsilon.eta3[i]<-min(1,(blikely3*((dnorm(canepsilon.eta3[i],a.eta3,1/(sqrt(zeta.eta4)))/(dnorm(v3[i-1],a.eta3,1/(sqrt(zeta.eta4)))))))) if(u=aprobepsilon.eta3[i]){ v3[i]<-v3[i-1] } canepsilon.xi3[i]<-w3[i-1]+rnorm(1,0,errepsilon.xi3) clikely3<-exp(loglik(x1[i],y1[i],v3[i],canepsilon.xi3[i],jan2245.exceed,thresh[25]) -loglik(x1[i],y1[i],v3[i],w3[i-1],jan2245.exceed,thresh[25]) +loglik(x2[i],y2[i],v3[i],canepsilon.xi3[i],feb2245.exceed,thresh[26]) -loglik(x2[i],y2[i],v3[i],w3[i-1],feb2245.exceed,thresh[26]) +loglik(x3[i],y3[i],v3[i],canepsilon.xi3[i],mar2245.exceed,thresh[27]) -loglik(x3[i],y3[i],v3[i],w3[i-1],mar2245.exceed,thresh[27]) +loglik(x4[i-1],y4[i-1],v3[i],canepsilon.xi3[i],apr2245.exceed,thresh[28]) -loglik(x4[i-1],y4[i-1],v3[i],w3[i-1],apr2245.exceed,thresh[28]) +loglik(x5[i-1],y5[i-1],v3[i],canepsilon.xi3[i],may2245.exceed,thresh[29]) -loglik(x5[i-1],y5[i-1],v3[i],w3[i-1],may2245.exceed,thresh[29]) +loglik(x6[i-1],y6[i-1],v3[i],canepsilon.xi3[i],jun2245.exceed,thresh[30]) -loglik(x6[i-1],y6[i-1],v3[i],w3[i-1],jun2245.exceed,thresh[30]) +loglik(x7[i-1],y7[i-1],v3[i],canepsilon.xi3[i],jul2245.exceed,thresh[31]) -loglik(x7[i-1],y7[i-1],v3[i],w3[i-1],jul2245.exceed,thresh[31]) +loglik(x8[i-1],y8[i-1],v3[i],canepsilon.xi3[i],aug2245.exceed,thresh[32]) -loglik(x8[i-1],y8[i-1],v3[i],w3[i-1],aug2245.exceed,thresh[32]) +loglik(x9[i-1],y9[i-1],v3[i],canepsilon.xi3[i],sep2245.exceed,thresh[33]) -loglik(x9[i-1],y9[i-1],v3[i],w3[i-1],sep2245.exceed,thresh[33]) +loglik(x10[i-1],y10[i-1],v3[i],canepsilon.xi3[i],oct2245.exceed,thresh[34]) -loglik(x10[i-1],y10[i-1],v3[i],w3[i-1],oct2245.exceed,thresh[34]) +loglik(x11[i-1],y11[i-1],v3[i],canepsilon.xi3[i],nov2245.exceed,thresh[35]) -loglik(x11[i-1],y11[i-1],v3[i],w3[i-1],nov2245.exceed,thresh[35]) +loglik(x12[i-1],y12[i-1],v3[i],canepsilon.xi3[i],dec2245.exceed,thresh[36]) -loglik(x12[i-1],y12[i-1],v3[i],w3[i-1],dec2245.exceed,thresh[36])) aprobepsilon.xi3[i]<-min(1,(clikely3*((dnorm(canepsilon.xi3[i],a.xi3,1/(sqrt(zeta.xi4)))/(dnorm(w3[i-1],a.xi3/(sqrt(zeta.xi4)))))))) if(u=aprobepsilon.xi3[i]){ w3[i]<-w3[i-1] } #eta.one.one[i]<-x1[i]+v1[i] #xi.one.one[i]<-y1[i]+w1[i] #***THIRD UPDATE DONE!***# a.eta4<-rnorm(1,(c*b+(zeta.eta4*(v1[i]+v2[i]+v3[i]+v4[i-1]+v5[i-1]+v6[i-1]+v7[i-1]+v8[i-1]+v9[i-1]+v10[i-1]+v11[i-1]+v12[i-1])))/(c+12*zeta.eta4),1/(sqrt(c+12*zeta.eta4))) a.xi4<-rnorm(1,(c*b+(zeta.xi4*(w1[i]+w2[i]+w3[i]+w4[i-1]+w5[i-1]+w6[i-1]+w7[i-1]+w8[i-1]+w9[i-1]+w10[i-1]+w11[i-1]+w12[i-1])))/(c+12*zeta.xi4),1/(sqrt(c+12*zeta.xi4))) zeta.eta5<-rgamma(1,f+6,1/(g+0.5*(((v1[i]-a.eta4)^2)+((v2[i]-a.eta4)^2)+((v3[i]-a.eta4)^2)+((v4[i-1]-a.eta4)^2)+((v5[i-1]-a.eta4)^2)+((v6[i-1]-a.eta4)^2)+((v7[i-1]-a.eta4)^2)+((v8[i-1]-a.eta4)^2)+((v9[i-1]-a.eta4)^2)+((v10[i-1]-a.eta4)^2)+((v11[i-1]-a.eta4)^2)+((v12[i-1]-a.eta4)^2)))) zeta.xi5<-rgamma(1,f+6,1/(g+0.5*(((w1[i]-a.xi4)^2)+((w2[i]-a.xi4)^2)+((w3[i]-a.xi4)^2)+((w4[i-1]-a.xi4)^2)+((w5[i-1]-a.xi4)^2)+((w6[i-1]-a.xi4)^2)+((w7[i-1]-a.xi4)^2)+((w8[i-1]-a.xi4)^2)+((w9[i-1]-a.xi4)^2)+((w10[i-1]-a.xi4)^2)+((w11[i-1]-a.xi4)^2)+((w12[i-1]-a.xi4)^2)))) tau.eta5<-rgamma(1,d+6,1/(e+0.5*(((x1[i])^2)+((x2[i])^2)+((x3[i])^2)+((x4[i-1])^2)+((x5[i-1])^2)+((x6[i-1])^2)+((x7[i-1])^2)+((x8[i-1])^2)+((x9[i-1])^2)+((x10[i-1])^2)+((x11[i-1])^2)+((x12[i-1])^2)))) tau.xi5<-rgamma(1,d+6,1/(e+0.5*(((y1[i])^2)+((y2[i])^2)+((y3[i])^2)+((y4[i-1])^2)+((y5[i-1])^2)+((y6[i-1])^2)+((y7[i-1])^2)+((y8[i-1])^2)+((y9[i-1])^2)+((y10[i-1])^2)+((y11[i-1])^2)+((y12[i-1])^2)))) cangamma.eta4[i]<-x4[i-1]+rnorm(1,0,errgamma.eta4) likely4<-exp(loglik(cangamma.eta4[i],y4[i-1],v1[i],w1[i],aprshef.exceed,thresh[4]) -loglik(x4[i-1],y4[i-1],v1[i],w1[i],aprshef.exceed,thresh[4]) +loglik(cangamma.eta4[i],y4[i-1],v2[i],w2[i],aprbrad.exceed,thresh[16]) -loglik(x4[i-1],y4[i-1],v2[i],w2[i],aprbrad.exceed,thresh[16]) +loglik(cangamma.eta4[i],y4[i-1],v3[i],w3[i],apr2245.exceed,thresh[28]) -loglik(x4[i-1],y4[i-1],v3[i],w3[i],apr2245.exceed,thresh[28]) +loglik(cangamma.eta4[i],y4[i-1],v4[i-1],w4[i-1],apr2255.exceed,thresh[40]) -loglik(x4[i-1],y4[i-1],v4[i-1],w4[i-1],apr2255.exceed,thresh[40]) +loglik(cangamma.eta4[i],y4[i-1],v5[i-1],w5[i-1],apr2416.exceed,thresh[52]) -loglik(x4[i-1],y4[i-1],v5[i-1],w5[i-1],apr2416.exceed,thresh[52]) +loglik(cangamma.eta4[i],y4[i-1],v6[i-1],w6[i-1],apr2432.exceed,thresh[64]) -loglik(x4[i-1],y4[i-1],v6[i-1],w6[i-1],apr2432.exceed,thresh[64]) +loglik(cangamma.eta4[i],y4[i-1],v7[i-1],w7[i-1],apr4036.exceed,thresh[76]) -loglik(x4[i-1],y4[i-1],v7[i-1],w7[i-1],apr4036.exceed,thresh[76]) +loglik(cangamma.eta4[i],y4[i-1],v8[i-1],w8[i-1],apr4206.exceed,thresh[88]) -loglik(x4[i-1],y4[i-1],v8[i-1],w8[i-1],apr4206.exceed,thresh[88]) +loglik(cangamma.eta4[i],y4[i-1],v9[i-1],w9[i-1],apr4238.exceed,thresh[100]) -loglik(x4[i-1],y4[i-1],v9[i-1],w9[i-1],apr4238.exceed,thresh[100]) +loglik(cangamma.eta4[i],y4[i-1],v10[i-1],w10[i-1],apr4617.exceed,thresh[112]) -loglik(x4[i-1],y4[i-1],v10[i-1],w10[i-1],apr4617.exceed,thresh[112]) +loglik(cangamma.eta4[i],y4[i-1],v11[i-1],w11[i-1],apr7221.exceed,thresh[124]) -loglik(x4[i-1],y4[i-1],v11[i-1],w11[i-1],apr7221.exceed,thresh[124]) +loglik(cangamma.eta4[i],y4[i-1],v12[i-1],w12[i-1],apr7377.exceed,thresh[136]) -loglik(x4[i-1],y4[i-1],v12[i-1],w12[i-1],apr7377.exceed,thresh[136])) aprobgamma.eta4[i]<-min(1,(likely4*((dnorm(cangamma.eta4[i],0,1/(sqrt(tau.eta5)))/(dnorm(x4[i-1],0,1/(sqrt(tau.eta5)))))))) u<-runif(1) if(u=aprobgamma.eta4[i]){ x4[i]<-x4[i-1] } cangamma.xi4[i]<-y4[i-1]+rnorm(1,0,errgamma.xi4) alikely4<-exp(loglik(x4[i],cangamma.xi4[i],v1[i],w1[i],aprshef.exceed,thresh[4]) -loglik(x4[i],y4[i-1],v1[i],w1[i],aprshef.exceed,thresh[4]) +loglik(x4[i],cangamma.xi4[i],v2[i],w2[i],aprbrad.exceed,thresh[16]) -loglik(x4[i],y4[i-1],v2[i],w2[i],aprbrad.exceed,thresh[16]) +loglik(x4[i],cangamma.xi4[i],v3[i],w3[i],apr2245.exceed,thresh[28]) -loglik(x4[i],y4[i-1],v3[i],w3[i],apr2245.exceed,thresh[28]) +loglik(x4[i],cangamma.xi4[i],v4[i-1],w4[i-1],apr2255.exceed,thresh[40]) -loglik(x4[i],y4[i-1],v4[i-1],w4[i-1],apr2255.exceed,thresh[40]) +loglik(x4[i],cangamma.xi4[i],v5[i-1],w5[i-1],apr2416.exceed,thresh[52]) -loglik(x4[i],y4[i-1],v5[i-1],w5[i-1],apr2416.exceed,thresh[52]) +loglik(x4[i],cangamma.xi4[i],v6[i-1],w6[i-1],apr2432.exceed,thresh[64]) -loglik(x4[i],y4[i-1],v6[i-1],w6[i-1],apr2432.exceed,thresh[64]) +loglik(x4[i],cangamma.xi4[i],v7[i-1],w7[i-1],apr4036.exceed,thresh[76]) -loglik(x4[i],y4[i-1],v7[i-1],w7[i-1],apr4036.exceed,thresh[76]) +loglik(x4[i],cangamma.xi4[i],v8[i-1],w8[i-1],apr4206.exceed,thresh[88]) -loglik(x4[i],y4[i-1],v8[i-1],w8[i-1],apr4206.exceed,thresh[88]) +loglik(x4[i],cangamma.xi4[i],v9[i-1],w9[i-1],apr4238.exceed,thresh[100]) -loglik(x4[i],y4[i-1],v9[i-1],w9[i-1],apr4238.exceed,thresh[100]) +loglik(x4[i],cangamma.xi4[i],v10[i-1],w10[i-1],apr4617.exceed,thresh[112]) -loglik(x4[i],y4[i-1],v10[i-1],w10[i-1],apr4617.exceed,thresh[112]) +loglik(x4[i],cangamma.xi4[i],v11[i-1],w11[i-1],apr7221.exceed,thresh[124]) -loglik(x4[i],y4[i-1],v11[i-1],w11[i-1],apr7221.exceed,thresh[124]) +loglik(x4[i],cangamma.xi4[i],v12[i-1],w12[i-1],apr7377.exceed,thresh[136]) -loglik(x4[i],y4[i-1],v12[i-1],w12[i-1],apr7377.exceed,thresh[136])) aprobgamma.xi4[i]<-min(1,(alikely4*((dnorm(cangamma.xi4[i],0,1/(sqrt(tau.xi5)))/(dnorm(y4[i-1],0,1/(sqrt(tau.xi5)))))))) if(u=aprobgamma.xi4[i]){ y4[i]<-y4[i-1] } canepsilon.eta4[i]<-v4[i-1]+rnorm(1,0,errepsilon.eta4) blikely4<-exp(loglik(x1[i],y1[i],canepsilon.eta4[i],w4[i-1],jan2255.exceed,thresh[37]) -loglik(x1[i],y1[i],v4[i-1],w4[i-1],jan2255.exceed,thresh[37]) +loglik(x2[i],y2[i],canepsilon.eta4[i],w4[i-1],feb2255.exceed,thresh[38]) -loglik(x2[i],y2[i],v4[i-1],w4[i-1],feb2255.exceed,thresh[38]) +loglik(x3[i],y3[i],canepsilon.eta4[i],w4[i-1],mar2255.exceed,thresh[39]) -loglik(x3[i],y3[i],v4[i-1],w4[i-1],mar2255.exceed,thresh[39]) +loglik(x4[i],y4[i],canepsilon.eta4[i],w4[i-1],apr2255.exceed,thresh[40]) -loglik(x4[i],y4[i],v4[i-1],w4[i-1],apr2255.exceed,thresh[40]) +loglik(x5[i-1],y5[i-1],canepsilon.eta4[i],w4[i-1],may2255.exceed,thresh[41]) -loglik(x5[i-1],y5[i-1],v4[i-1],w4[i-1],may2255.exceed,thresh[41]) +loglik(x6[i-1],y6[i-1],canepsilon.eta4[i],w4[i-1],jun2255.exceed,thresh[42]) -loglik(x6[i-1],y6[i-1],v4[i-1],w4[i-1],jun2255.exceed,thresh[42]) +loglik(x7[i-1],y7[i-1],canepsilon.eta4[i],w4[i-1],jul2255.exceed,thresh[43]) -loglik(x7[i-1],y7[i-1],v4[i-1],w4[i-1],jul2255.exceed,thresh[43]) +loglik(x8[i-1],y8[i-1],canepsilon.eta4[i],w4[i-1],aug2255.exceed,thresh[44]) -loglik(x8[i-1],y8[i-1],v4[i-1],w4[i-1],aug2255.exceed,thresh[44]) +loglik(x9[i-1],y9[i-1],canepsilon.eta4[i],w4[i-1],sep2255.exceed,thresh[45]) -loglik(x9[i-1],y9[i-1],v4[i-1],w4[i-1],sep2255.exceed,thresh[45]) +loglik(x10[i-1],y10[i-1],canepsilon.eta4[i],w4[i-1],oct2255.exceed,thresh[46]) -loglik(x10[i-1],y10[i-1],v4[i-1],w4[i-1],oct2255.exceed,thresh[46]) +loglik(x11[i-1],y11[i-1],canepsilon.eta4[i],w4[i-1],nov2255.exceed,thresh[47]) -loglik(x11[i-1],y11[i-1],v4[i-1],w4[i-1],nov2255.exceed,thresh[47]) +loglik(x12[i-1],y12[i-1],canepsilon.eta4[i],w4[i-1],dec2255.exceed,thresh[48]) -loglik(x12[i-1],y12[i-1],v4[i-1],w4[i-1],dec2255.exceed,thresh[48])) aprobepsilon.eta4[i]<-min(1,(blikely4*((dnorm(canepsilon.eta4[i],a.eta4,1/(sqrt(zeta.eta5)))/(dnorm(v4[i-1],a.eta4,1/(sqrt(zeta.eta5)))))))) if(u=aprobepsilon.eta4[i]){ v4[i]<-v4[i-1] } canepsilon.xi4[i]<-w4[i-1]+rnorm(1,0,errepsilon.xi4) clikely4<-exp(loglik(x1[i],y1[i],v4[i],canepsilon.xi4[i],jan2255.exceed,thresh[37]) -loglik(x1[i],y1[i],v4[i],w4[i-1],jan2255.exceed,thresh[37]) +loglik(x2[i],y2[i],v4[i],canepsilon.xi4[i],feb2255.exceed,thresh[38]) -loglik(x2[i],y2[i],v4[i],w4[i-1],feb2255.exceed,thresh[38]) +loglik(x3[i],y3[i],v4[i],canepsilon.xi4[i],mar2255.exceed,thresh[39]) -loglik(x3[i],y3[i],v4[i],w4[i-1],mar2255.exceed,thresh[39]) +loglik(x4[i],y4[i],v4[i],canepsilon.xi4[i],apr2255.exceed,thresh[40]) -loglik(x4[i],y4[i],v4[i],w4[i-1],apr2255.exceed,thresh[40]) +loglik(x5[i-1],y5[i-1],v4[i],canepsilon.xi4[i],may2255.exceed,thresh[41]) -loglik(x5[i-1],y5[i-1],v4[i],w4[i-1],may2255.exceed,thresh[41]) +loglik(x6[i-1],y6[i-1],v4[i],canepsilon.xi4[i],jun2255.exceed,thresh[42]) -loglik(x6[i-1],y6[i-1],v4[i],w4[i-1],jun2255.exceed,thresh[42]) +loglik(x7[i-1],y7[i-1],v4[i],canepsilon.xi4[i],jul2255.exceed,thresh[43]) -loglik(x7[i-1],y7[i-1],v4[i],w4[i-1],jul2255.exceed,thresh[43]) +loglik(x8[i-1],y8[i-1],v4[i],canepsilon.xi4[i],aug2255.exceed,thresh[44]) -loglik(x8[i-1],y8[i-1],v4[i],w4[i-1],aug2255.exceed,thresh[44]) +loglik(x9[i-1],y9[i-1],v4[i],canepsilon.xi4[i],sep2255.exceed,thresh[45]) -loglik(x9[i-1],y9[i-1],v4[i],w4[i-1],sep2255.exceed,thresh[45]) +loglik(x10[i-1],y10[i-1],v4[i],canepsilon.xi4[i],oct2255.exceed,thresh[46]) -loglik(x10[i-1],y10[i-1],v4[i],w4[i-1],oct2255.exceed,thresh[46]) +loglik(x11[i-1],y11[i-1],v4[i],canepsilon.xi4[i],nov2255.exceed,thresh[47]) -loglik(x11[i-1],y11[i-1],v4[i],w4[i-1],nov2255.exceed,thresh[47]) +loglik(x12[i-1],y12[i-1],v4[i],canepsilon.xi4[i],dec2255.exceed,thresh[48]) -loglik(x12[i-1],y12[i-1],v4[i],w4[i-1],dec2255.exceed,thresh[48])) aprobepsilon.xi4[i]<-min(1,(clikely4*((dnorm(canepsilon.xi4[i],a.xi4,1/(sqrt(zeta.xi5)))/(dnorm(w4[i-1],a.xi4/(sqrt(zeta.xi5)))))))) if(u=aprobepsilon.xi4[i]){ w4[i]<-w4[i-1] } #eta.one.one[i]<-x1[i]+v1[i] #xi.one.one[i]<-y1[i]+w1[i] #***FOURTH UPDATE DONE!***# a.eta5<-rnorm(1,(c*b+(zeta.eta5*(v1[i]+v2[i]+v3[i]+v4[i]+v5[i-1]+v6[i-1]+v7[i-1]+v8[i-1]+v9[i-1]+v10[i-1]+v11[i-1]+v12[i-1])))/(c+12*zeta.eta5),1/(sqrt(c+12*zeta.eta5))) a.xi5<-rnorm(1,(c*b+(zeta.xi5*(w1[i]+w2[i]+w3[i]+w4[i]+w5[i-1]+w6[i-1]+w7[i-1]+w8[i-1]+w9[i-1]+w10[i-1]+w11[i-1]+w12[i-1])))/(c+12*zeta.xi5),1/(sqrt(c+12*zeta.xi5))) zeta.eta6<-rgamma(1,f+6,1/(g+0.5*(((v1[i]-a.eta5)^2)+((v2[i]-a.eta5)^2)+((v3[i]-a.eta5)^2)+((v4[i]-a.eta5)^2)+((v5[i-1]-a.eta5)^2)+((v6[i-1]-a.eta5)^2)+((v7[i-1]-a.eta5)^2)+((v8[i-1]-a.eta5)^2)+((v9[i-1]-a.eta5)^2)+((v10[i-1]-a.eta5)^2)+((v11[i-1]-a.eta5)^2)+((v12[i-1]-a.eta5)^2)))) zeta.xi6<-rgamma(1,f+6,1/(g+0.5*(((w1[i]-a.xi5)^2)+((w2[i]-a.xi5)^2)+((w3[i]-a.xi5)^2)+((w4[i]-a.xi5)^2)+((w5[i-1]-a.xi5)^2)+((w6[i-1]-a.xi5)^2)+((w7[i-1]-a.xi5)^2)+((w8[i-1]-a.xi5)^2)+((w9[i-1]-a.xi5)^2)+((w10[i-1]-a.xi5)^2)+((w11[i-1]-a.xi5)^2)+((w12[i-1]-a.xi5)^2)))) tau.eta6<-rgamma(1,d+6,1/(e+0.5*(((x1[i])^2)+((x2[i])^2)+((x3[i])^2)+((x4[i])^2)+((x5[i-1])^2)+((x6[i-1])^2)+((x7[i-1])^2)+((x8[i-1])^2)+((x9[i-1])^2)+((x10[i-1])^2)+((x11[i-1])^2)+((x12[i-1])^2)))) tau.xi6<-rgamma(1,d+6,1/(e+0.5*(((y1[i])^2)+((y2[i])^2)+((y3[i])^2)+((y4[i])^2)+((y5[i-1])^2)+((y6[i-1])^2)+((y7[i-1])^2)+((y8[i-1])^2)+((y9[i-1])^2)+((y10[i-1])^2)+((y11[i-1])^2)+((y12[i-1])^2)))) cangamma.eta5[i]<-x5[i-1]+rnorm(1,0,errgamma.eta5) likely5<-exp(loglik(cangamma.eta5[i],y5[i-1],v1[i],w1[i],mayshef.exceed,thresh[5]) -loglik(x5[i-1],y5[i-1],v1[i],w1[i],mayshef.exceed,thresh[5]) +loglik(cangamma.eta5[i],y5[i-1],v2[i],w2[i],maybrad.exceed,thresh[17]) -loglik(x5[i-1],y5[i-1],v2[i],w2[i],maybrad.exceed,thresh[17]) +loglik(cangamma.eta5[i],y5[i-1],v3[i],w3[i],may2245.exceed,thresh[29]) -loglik(x5[i-1],y5[i-1],v3[i],w3[i],may2245.exceed,thresh[29]) +loglik(cangamma.eta5[i],y5[i-1],v4[i],w4[i],may2255.exceed,thresh[41]) -loglik(x5[i-1],y5[i-1],v4[i-1],w4[i-1],may2255.exceed,thresh[41]) +loglik(cangamma.eta5[i],y5[i-1],v5[i-1],w5[i-1],may2416.exceed,thresh[53]) -loglik(x5[i-1],y5[i-1],v5[i-1],w5[i-1],may2416.exceed,thresh[53]) +loglik(cangamma.eta5[i],y5[i-1],v6[i-1],w6[i-1],may2432.exceed,thresh[65]) -loglik(x5[i-1],y5[i-1],v6[i-1],w6[i-1],may2432.exceed,thresh[65]) +loglik(cangamma.eta5[i],y5[i-1],v7[i-1],w7[i-1],may4036.exceed,thresh[77]) -loglik(x5[i-1],y5[i-1],v7[i-1],w7[i-1],may4036.exceed,thresh[77]) +loglik(cangamma.eta5[i],y5[i-1],v8[i-1],w8[i-1],may4206.exceed,thresh[89]) -loglik(x5[i-1],y5[i-1],v8[i-1],w8[i-1],may4206.exceed,thresh[89]) +loglik(cangamma.eta5[i],y5[i-1],v9[i-1],w9[i-1],may4238.exceed,thresh[101]) -loglik(x5[i-1],y5[i-1],v9[i-1],w9[i-1],may4238.exceed,thresh[101]) +loglik(cangamma.eta5[i],y5[i-1],v10[i-1],w10[i-1],may4617.exceed,thresh[113]) -loglik(x5[i-1],y5[i-1],v10[i-1],w10[i-1],may4617.exceed,thresh[113]) +loglik(cangamma.eta5[i],y5[i-1],v11[i-1],w11[i-1],may7221.exceed,thresh[125]) -loglik(x5[i-1],y5[i-1],v11[i-1],w11[i-1],may7221.exceed,thresh[125]) +loglik(cangamma.eta5[i],y5[i-1],v12[i-1],w12[i-1],may7377.exceed,thresh[137]) -loglik(x5[i-1],y5[i-1],v12[i-1],w12[i-1],may7377.exceed,thresh[137])) aprobgamma.eta5[i]<-min(1,(likely5*((dnorm(cangamma.eta5[i],0,1/(sqrt(tau.eta6)))/(dnorm(x5[i-1],0,1/(sqrt(tau.eta6)))))))) u<-runif(1) if(u=aprobgamma.eta5[i]){ x5[i]<-x5[i-1] } cangamma.xi5[i]<-y5[i-1]+rnorm(1,0,errgamma.xi5) alikely5<-exp(loglik(x5[i],cangamma.xi5[i],v1[i],w1[i],mayshef.exceed,thresh[5]) -loglik(x5[i],y5[i-1],v1[i],w1[i],mayshef.exceed,thresh[5]) +loglik(x5[i],cangamma.xi5[i],v2[i],w2[i],maybrad.exceed,thresh[17]) -loglik(x5[i],y5[i-1],v2[i],w2[i],maybrad.exceed,thresh[17]) +loglik(x5[i],cangamma.xi5[i],v3[i],w3[i],may2245.exceed,thresh[29]) -loglik(x5[i],y5[i-1],v3[i],w3[i],may2245.exceed,thresh[29]) +loglik(x5[i],cangamma.xi5[i],v4[i],w4[i],may2255.exceed,thresh[41]) -loglik(x5[i],y5[i-1],v4[i],w4[i],may2255.exceed,thresh[41]) +loglik(x5[i],cangamma.xi5[i],v5[i-1],w5[i-1],may2416.exceed,thresh[53]) -loglik(x5[i],y5[i-1],v5[i-1],w5[i-1],may2416.exceed,thresh[53]) +loglik(x5[i],cangamma.xi5[i],v6[i-1],w6[i-1],may2432.exceed,thresh[65]) -loglik(x5[i],y5[i-1],v6[i-1],w6[i-1],may2432.exceed,thresh[65]) +loglik(x5[i],cangamma.xi5[i],v7[i-1],w7[i-1],may4036.exceed,thresh[77]) -loglik(x5[i],y5[i-1],v7[i-1],w7[i-1],may4036.exceed,thresh[77]) +loglik(x5[i],cangamma.xi5[i],v8[i-1],w8[i-1],may4206.exceed,thresh[89]) -loglik(x5[i],y5[i-1],v8[i-1],w8[i-1],may4206.exceed,thresh[89]) +loglik(x5[i],cangamma.xi5[i],v9[i-1],w9[i-1],may4238.exceed,thresh[101]) -loglik(x5[i],y5[i-1],v9[i-1],w9[i-1],may4238.exceed,thresh[101]) +loglik(x5[i],cangamma.xi5[i],v10[i-1],w10[i-1],may4617.exceed,thresh[113]) -loglik(x5[i],y5[i-1],v10[i-1],w10[i-1],may4617.exceed,thresh[113]) +loglik(x5[i],cangamma.xi5[i],v11[i-1],w11[i-1],may7221.exceed,thresh[125]) -loglik(x5[i],y5[i-1],v11[i-1],w11[i-1],may7221.exceed,thresh[125]) +loglik(x5[i],cangamma.xi5[i],v12[i-1],w12[i-1],may7377.exceed,thresh[137]) -loglik(x5[i],y5[i-1],v12[i-1],w12[i-1],may7377.exceed,thresh[137])) aprobgamma.xi5[i]<-min(1,(alikely5*((dnorm(cangamma.xi5[i],0,1/(sqrt(tau.xi6)))/(dnorm(y5[i-1],0,1/(sqrt(tau.xi6)))))))) if(u=aprobgamma.xi5[i]){ y5[i]<-y5[i-1] } canepsilon.eta5[i]<-v5[i-1]+rnorm(1,0,errepsilon.eta5) blikely5<-exp(loglik(x1[i],y1[i],canepsilon.eta5[i],w5[i-1],jan2416.exceed,thresh[49]) -loglik(x1[i],y1[i],v5[i-1],w5[i-1],jan2416.exceed,thresh[49]) +loglik(x2[i],y2[i],canepsilon.eta5[i],w5[i-1],feb2416.exceed,thresh[50]) -loglik(x2[i],y2[i],v5[i-1],w5[i-1],feb2416.exceed,thresh[50]) +loglik(x3[i],y3[i],canepsilon.eta5[i],w5[i-1],mar2416.exceed,thresh[51]) -loglik(x3[i],y3[i],v5[i-1],w5[i-1],mar2416.exceed,thresh[51]) +loglik(x4[i],y4[i],canepsilon.eta5[i],w5[i-1],apr2416.exceed,thresh[52]) -loglik(x4[i],y4[i],v5[i-1],w5[i-1],apr2416.exceed,thresh[52]) +loglik(x5[i],y5[i],canepsilon.eta5[i],w5[i-1],may2416.exceed,thresh[53]) -loglik(x5[i],y5[i],v5[i-1],w5[i-1],may2416.exceed,thresh[53]) +loglik(x6[i-1],y6[i-1],canepsilon.eta5[i],w5[i-1],jun2416.exceed,thresh[54]) -loglik(x6[i-1],y6[i-1],v5[i-1],w5[i-1],jun2416.exceed,thresh[54]) +loglik(x7[i-1],y7[i-1],canepsilon.eta5[i],w5[i-1],jul2416.exceed,thresh[55]) -loglik(x7[i-1],y7[i-1],v5[i-1],w5[i-1],jul2416.exceed,thresh[55]) +loglik(x8[i-1],y8[i-1],canepsilon.eta5[i],w5[i-1],aug2416.exceed,thresh[56]) -loglik(x8[i-1],y8[i-1],v5[i-1],w5[i-1],aug2416.exceed,thresh[56]) +loglik(x9[i-1],y9[i-1],canepsilon.eta5[i],w5[i-1],sep2416.exceed,thresh[57]) -loglik(x9[i-1],y9[i-1],v5[i-1],w5[i-1],sep2416.exceed,thresh[57]) +loglik(x10[i-1],y10[i-1],canepsilon.eta5[i],w5[i-1],oct2416.exceed,thresh[58]) -loglik(x10[i-1],y10[i-1],v5[i-1],w5[i-1],oct2416.exceed,thresh[58]) +loglik(x11[i-1],y11[i-1],canepsilon.eta5[i],w5[i-1],nov2416.exceed,thresh[59]) -loglik(x11[i-1],y11[i-1],v5[i-1],w5[i-1],nov2416.exceed,thresh[59]) +loglik(x12[i-1],y12[i-1],canepsilon.eta5[i],w5[i-1],dec2416.exceed,thresh[60]) -loglik(x12[i-1],y12[i-1],v5[i-1],w5[i-1],dec2416.exceed,thresh[60])) aprobepsilon.eta5[i]<-min(1,(blikely5*((dnorm(canepsilon.eta5[i],a.eta5,1/(sqrt(zeta.eta6)))/(dnorm(v5[i-1],a.eta5,1/(sqrt(zeta.eta6)))))))) if(u=aprobepsilon.eta5[i]){ v5[i]<-v5[i-1] } canepsilon.xi5[i]<-w5[i-1]+rnorm(1,0,errepsilon.xi5) clikely5<-exp(loglik(x1[i],y1[i],v5[i],canepsilon.xi5[i],jan2416.exceed,thresh[49]) -loglik(x1[i],y1[i],v5[i],w5[i-1],jan2416.exceed,thresh[49]) +loglik(x2[i],y2[i],v5[i],canepsilon.xi5[i],feb2416.exceed,thresh[50]) -loglik(x2[i],y2[i],v5[i],w5[i-1],feb2416.exceed,thresh[50]) +loglik(x3[i],y3[i],v5[i],canepsilon.xi5[i],mar2416.exceed,thresh[51]) -loglik(x3[i],y3[i],v5[i],w5[i-1],mar2416.exceed,thresh[51]) +loglik(x4[i],y4[i],v5[i],canepsilon.xi5[i],apr2416.exceed,thresh[52]) -loglik(x4[i],y4[i],v5[i],w5[i-1],apr2416.exceed,thresh[52]) +loglik(x5[i],y5[i],v5[i],canepsilon.xi5[i],may2416.exceed,thresh[53]) -loglik(x5[i],y5[i],v5[i],w5[i-1],may2416.exceed,thresh[53]) +loglik(x6[i-1],y6[i-1],v5[i],canepsilon.xi5[i],jun2416.exceed,thresh[54]) -loglik(x6[i-1],y6[i-1],v5[i],w5[i-1],jun2416.exceed,thresh[54]) +loglik(x7[i-1],y7[i-1],v5[i],canepsilon.xi5[i],jul2416.exceed,thresh[55]) -loglik(x7[i-1],y7[i-1],v5[i],w5[i-1],jul2416.exceed,thresh[55]) +loglik(x8[i-1],y8[i-1],v5[i],canepsilon.xi5[i],aug2416.exceed,thresh[56]) -loglik(x8[i-1],y8[i-1],v5[i],w5[i-1],aug2416.exceed,thresh[56]) +loglik(x9[i-1],y9[i-1],v5[i],canepsilon.xi5[i],sep2416.exceed,thresh[57]) -loglik(x9[i-1],y9[i-1],v5[i],w5[i-1],sep2416.exceed,thresh[57]) +loglik(x10[i-1],y10[i-1],v5[i],canepsilon.xi5[i],oct2416.exceed,thresh[58]) -loglik(x10[i-1],y10[i-1],v5[i],w5[i-1],oct2416.exceed,thresh[58]) +loglik(x11[i-1],y11[i-1],v5[i],canepsilon.xi5[i],nov2416.exceed,thresh[59]) -loglik(x11[i-1],y11[i-1],v5[i],w5[i-1],nov2416.exceed,thresh[59]) +loglik(x12[i-1],y12[i-1],v5[i],canepsilon.xi5[i],dec2416.exceed,thresh[60]) -loglik(x12[i-1],y12[i-1],v5[i],w5[i-1],dec2416.exceed,thresh[60])) aprobepsilon.xi5[i]<-min(1,(clikely5*((dnorm(canepsilon.xi5[i],a.xi5,1/(sqrt(zeta.xi6)))/(dnorm(w5[i-1],a.xi5/(sqrt(zeta.xi6)))))))) if(u=aprobepsilon.xi5[i]){ w5[i]<-w5[i-1] } #eta.one.one[i]<-x1[i]+v1[i] #xi.one.one[i]<-y1[i]+w1[i] #***FIFTH UPDATE DONE!***# a.eta6<-rnorm(1,(c*b+(zeta.eta6*(v1[i]+v2[i]+v3[i]+v4[i]+v5[i]+v6[i-1]+v7[i-1]+v8[i-1]+v9[i-1]+v10[i-1]+v11[i-1]+v12[i-1])))/(c+12*zeta.eta6),1/(sqrt(c+12*zeta.eta6))) a.xi6<-rnorm(1,(c*b+(zeta.xi6*(w1[i]+w2[i]+w3[i]+w4[i]+w5[i]+w6[i-1]+w7[i-1]+w8[i-1]+w9[i-1]+w10[i-1]+w11[i-1]+w12[i-1])))/(c+12*zeta.xi6),1/(sqrt(c+12*zeta.xi6))) zeta.eta7<-rgamma(1,f+6,1/(g+0.5*(((v1[i]-a.eta6)^2)+((v2[i]-a.eta6)^2)+((v3[i]-a.eta6)^2)+((v4[i]-a.eta6)^2)+((v5[i]-a.eta6)^2)+((v6[i-1]-a.eta6)^2)+((v7[i-1]-a.eta6)^2)+((v8[i-1]-a.eta6)^2)+((v9[i-1]-a.eta6)^2)+((v10[i-1]-a.eta6)^2)+((v11[i-1]-a.eta6)^2)+((v12[i-1]-a.eta6)^2)))) zeta.xi7<-rgamma(1,f+6,1/(g+0.5*(((w1[i]-a.xi6)^2)+((w2[i]-a.xi6)^2)+((w3[i]-a.xi6)^2)+((w4[i]-a.xi6)^2)+((w5[i]-a.xi6)^2)+((w6[i-1]-a.xi6)^2)+((w7[i-1]-a.xi6)^2)+((w8[i-1]-a.xi6)^2)+((w9[i-1]-a.xi6)^2)+((w10[i-1]-a.xi6)^2)+((w11[i-1]-a.xi6)^2)+((w12[i-1]-a.xi6)^2)))) tau.eta7<-rgamma(1,d+6,1/(e+0.5*(((x1[i])^2)+((x2[i])^2)+((x3[i])^2)+((x4[i])^2)+((x5[i])^2)+((x6[i-1])^2)+((x7[i-1])^2)+((x8[i-1])^2)+((x9[i-1])^2)+((x10[i-1])^2)+((x11[i-1])^2)+((x12[i-1])^2)))) tau.xi7<-rgamma(1,d+6,1/(e+0.5*(((y1[i])^2)+((y2[i])^2)+((y3[i])^2)+((y4[i])^2)+((y5[i])^2)+((y6[i-1])^2)+((y7[i-1])^2)+((y8[i-1])^2)+((y9[i-1])^2)+((y10[i-1])^2)+((y11[i-1])^2)+((y12[i-1])^2)))) cangamma.eta6[i]<-x6[i-1]+rnorm(1,0,errgamma.eta6) likely6<-exp(loglik(cangamma.eta6[i],y6[i-1],v1[i],w1[i],junshef.exceed,thresh[6]) -loglik(x6[i-1],y6[i-1],v1[i],w1[i],junshef.exceed,thresh[6]) +loglik(cangamma.eta6[i],y6[i-1],v2[i],w2[i],junbrad.exceed,thresh[18]) -loglik(x6[i-1],y6[i-1],v2[i],w2[i],junbrad.exceed,thresh[18]) +loglik(cangamma.eta6[i],y6[i-1],v3[i],w3[i],jun2245.exceed,thresh[30]) -loglik(x6[i-1],y6[i-1],v3[i],w3[i],jun2245.exceed,thresh[30]) +loglik(cangamma.eta6[i],y6[i-1],v4[i],w4[i],jun2255.exceed,thresh[42]) -loglik(x6[i-1],y6[i-1],v4[i-1],w4[i-1],jun2255.exceed,thresh[42]) +loglik(cangamma.eta6[i],y6[i-1],v5[i],w5[i],jun2416.exceed,thresh[54]) -loglik(x6[i-1],y6[i-1],v5[i],w5[i],jun2416.exceed,thresh[54]) +loglik(cangamma.eta6[i],y6[i-1],v6[i-1],w6[i-1],jun2432.exceed,thresh[66]) -loglik(x6[i-1],y6[i-1],v6[i-1],w6[i-1],jun2432.exceed,thresh[66]) +loglik(cangamma.eta6[i],y6[i-1],v7[i-1],w7[i-1],jun4036.exceed,thresh[78]) -loglik(x6[i-1],y6[i-1],v7[i-1],w7[i-1],jun4036.exceed,thresh[78]) +loglik(cangamma.eta6[i],y6[i-1],v8[i-1],w8[i-1],jun4206.exceed,thresh[90]) -loglik(x6[i-1],y6[i-1],v8[i-1],w8[i-1],jun4206.exceed,thresh[90]) +loglik(cangamma.eta6[i],y6[i-1],v9[i-1],w9[i-1],jun4238.exceed,thresh[102]) -loglik(x6[i-1],y6[i-1],v9[i-1],w9[i-1],jun4238.exceed,thresh[102]) +loglik(cangamma.eta6[i],y6[i-1],v10[i-1],w10[i-1],jun4617.exceed,thresh[114]) -loglik(x6[i-1],y6[i-1],v10[i-1],w10[i-1],jun4617.exceed,thresh[114]) +loglik(cangamma.eta6[i],y6[i-1],v11[i-1],w11[i-1],jun7221.exceed,thresh[126]) -loglik(x6[i-1],y6[i-1],v11[i-1],w11[i-1],jun7221.exceed,thresh[126]) +loglik(cangamma.eta6[i],y6[i-1],v12[i-1],w12[i-1],jun7377.exceed,thresh[138]) -loglik(x6[i-1],y6[i-1],v12[i-1],w12[i-1],jun7377.exceed,thresh[138])) aprobgamma.eta6[i]<-min(1,(likely6*((dnorm(cangamma.eta6[i],0,1/(sqrt(tau.eta7)))/(dnorm(x6[i-1],0,1/(sqrt(tau.eta7)))))))) u<-runif(1) if(u=aprobgamma.eta6[i]){ x6[i]<-x6[i-1] } cangamma.xi6[i]<-y6[i-1]+rnorm(1,0,errgamma.xi6) alikely6<-exp(loglik(x6[i],cangamma.xi6[i],v1[i],w1[i],junshef.exceed,thresh[6]) -loglik(x6[i],y6[i-1],v1[i],w1[i],junshef.exceed,thresh[6]) +loglik(x6[i],cangamma.xi6[i],v2[i],w2[i],junbrad.exceed,thresh[18]) -loglik(x6[i],y6[i-1],v2[i],w2[i],junbrad.exceed,thresh[18]) +loglik(x6[i],cangamma.xi6[i],v3[i],w3[i],jun2245.exceed,thresh[30]) -loglik(x6[i],y6[i-1],v3[i],w3[i],jun2245.exceed,thresh[30]) +loglik(x6[i],cangamma.xi6[i],v4[i],w4[i],jun2255.exceed,thresh[42]) -loglik(x6[i],y6[i-1],v4[i],w4[i],jun2255.exceed,thresh[42]) +loglik(x6[i],cangamma.xi6[i],v5[i],w5[i],jun2416.exceed,thresh[54]) -loglik(x6[i],y6[i-1],v5[i],w5[i],jun2416.exceed,thresh[54]) +loglik(x6[i],cangamma.xi6[i],v6[i-1],w6[i-1],jun2432.exceed,thresh[66]) -loglik(x6[i],y6[i-1],v6[i-1],w6[i-1],jun2432.exceed,thresh[66]) +loglik(x6[i],cangamma.xi6[i],v7[i-1],w7[i-1],jun4036.exceed,thresh[78]) -loglik(x6[i],y6[i-1],v7[i-1],w7[i-1],jun4036.exceed,thresh[78]) +loglik(x6[i],cangamma.xi6[i],v8[i-1],w8[i-1],jun4206.exceed,thresh[90]) -loglik(x6[i],y6[i-1],v8[i-1],w8[i-1],jun4206.exceed,thresh[90]) +loglik(x6[i],cangamma.xi6[i],v9[i-1],w9[i-1],jun4238.exceed,thresh[102]) -loglik(x6[i],y6[i-1],v9[i-1],w9[i-1],jun4238.exceed,thresh[102]) +loglik(x6[i],cangamma.xi6[i],v10[i-1],w10[i-1],jun4617.exceed,thresh[114]) -loglik(x6[i],y6[i-1],v10[i-1],w10[i-1],jun4617.exceed,thresh[114]) +loglik(x6[i],cangamma.xi6[i],v11[i-1],w11[i-1],jun7221.exceed,thresh[126]) -loglik(x6[i],y6[i-1],v11[i-1],w11[i-1],jun7221.exceed,thresh[126]) +loglik(x6[i],cangamma.xi6[i],v12[i-1],w12[i-1],jun7377.exceed,thresh[138]) -loglik(x6[i],y6[i-1],v12[i-1],w12[i-1],jun7377.exceed,thresh[138])) aprobgamma.xi6[i]<-min(1,(alikely6*((dnorm(cangamma.xi6[i],0,1/(sqrt(tau.xi7)))/(dnorm(y6[i-1],0,1/(sqrt(tau.xi7)))))))) if(u=aprobgamma.xi6[i]){ y6[i]<-y6[i-1] } canepsilon.eta6[i]<-v6[i-1]+rnorm(1,0,errepsilon.eta6) blikely6<-exp(loglik(x1[i],y1[i],canepsilon.eta6[i],w6[i-1],jan2432.exceed,thresh[61]) -loglik(x1[i],y1[i],v6[i-1],w6[i-1],jan2432.exceed,thresh[61]) +loglik(x2[i],y2[i],canepsilon.eta6[i],w6[i-1],feb2432.exceed,thresh[62]) -loglik(x2[i],y2[i],v6[i-1],w6[i-1],feb2432.exceed,thresh[62]) +loglik(x3[i],y3[i],canepsilon.eta6[i],w6[i-1],mar2432.exceed,thresh[63]) -loglik(x3[i],y3[i],v6[i-1],w6[i-1],mar2432.exceed,thresh[63]) +loglik(x4[i],y4[i],canepsilon.eta6[i],w6[i-1],apr2432.exceed,thresh[64]) -loglik(x4[i],y4[i],v6[i-1],w6[i-1],apr2432.exceed,thresh[64]) +loglik(x5[i],y5[i],canepsilon.eta6[i],w6[i-1],may2432.exceed,thresh[65]) -loglik(x5[i],y5[i],v6[i-1],w6[i-1],may2432.exceed,thresh[65]) +loglik(x6[i],y6[i],canepsilon.eta6[i],w6[i-1],jun2432.exceed,thresh[66]) -loglik(x6[i],y6[i],v6[i-1],w6[i-1],jun2432.exceed,thresh[66]) +loglik(x7[i-1],y7[i-1],canepsilon.eta6[i],w6[i-1],jul2432.exceed,thresh[67]) -loglik(x7[i-1],y7[i-1],v6[i-1],w6[i-1],jul2432.exceed,thresh[67]) +loglik(x8[i-1],y8[i-1],canepsilon.eta6[i],w6[i-1],aug2432.exceed,thresh[68]) -loglik(x8[i-1],y8[i-1],v6[i-1],w6[i-1],aug2432.exceed,thresh[68]) +loglik(x9[i-1],y9[i-1],canepsilon.eta6[i],w6[i-1],sep2432.exceed,thresh[69]) -loglik(x9[i-1],y9[i-1],v6[i-1],w6[i-1],sep2432.exceed,thresh[69]) +loglik(x10[i-1],y10[i-1],canepsilon.eta6[i],w6[i-1],oct2432.exceed,thresh[70]) -loglik(x10[i-1],y10[i-1],v6[i-1],w6[i-1],oct2432.exceed,thresh[70]) +loglik(x11[i-1],y11[i-1],canepsilon.eta6[i],w6[i-1],nov2432.exceed,thresh[71]) -loglik(x11[i-1],y11[i-1],v6[i-1],w6[i-1],nov2432.exceed,thresh[71]) +loglik(x12[i-1],y12[i-1],canepsilon.eta6[i],w6[i-1],dec2432.exceed,thresh[72]) -loglik(x12[i-1],y12[i-1],v6[i-1],w6[i-1],dec2432.exceed,thresh[72])) aprobepsilon.eta6[i]<-min(1,(blikely6*((dnorm(canepsilon.eta6[i],a.eta6,1/(sqrt(zeta.eta7)))/(dnorm(v6[i-1],a.eta6,1/(sqrt(zeta.eta7)))))))) if(u=aprobepsilon.eta6[i]){ v6[i]<-v6[i-1] } canepsilon.xi6[i]<-w6[i-1]+rnorm(1,0,errepsilon.xi6) clikely6<-exp(loglik(x1[i],y1[i],v6[i],canepsilon.xi6[i],jan2432.exceed,thresh[61]) -loglik(x1[i],y1[i],v6[i],w6[i-1],jan2432.exceed,thresh[61]) +loglik(x2[i],y2[i],v6[i],canepsilon.xi6[i],feb2432.exceed,thresh[62]) -loglik(x2[i],y2[i],v6[i],w6[i-1],feb2432.exceed,thresh[62]) +loglik(x3[i],y3[i],v6[i],canepsilon.xi6[i],mar2432.exceed,thresh[63]) -loglik(x3[i],y3[i],v6[i],w6[i-1],mar2432.exceed,thresh[63]) +loglik(x4[i],y4[i],v6[i],canepsilon.xi6[i],apr2432.exceed,thresh[64]) -loglik(x4[i],y4[i],v6[i],w6[i-1],apr2432.exceed,thresh[64]) +loglik(x5[i],y5[i],v6[i],canepsilon.xi6[i],may2432.exceed,thresh[65]) -loglik(x5[i],y5[i],v6[i],w6[i-1],may2432.exceed,thresh[65]) +loglik(x6[i],y6[i],v6[i],canepsilon.xi6[i],jun2432.exceed,thresh[66]) -loglik(x6[i],y6[i],v6[i],w6[i-1],jun2432.exceed,thresh[66]) +loglik(x7[i-1],y7[i-1],v6[i],canepsilon.xi6[i],jul2432.exceed,thresh[67]) -loglik(x7[i-1],y7[i-1],v6[i],w6[i-1],jul2432.exceed,thresh[67]) +loglik(x8[i-1],y8[i-1],v6[i],canepsilon.xi6[i],aug2432.exceed,thresh[68]) -loglik(x8[i-1],y8[i-1],v6[i],w6[i-1],aug2432.exceed,thresh[68]) +loglik(x9[i-1],y9[i-1],v6[i],canepsilon.xi6[i],sep2432.exceed,thresh[69]) -loglik(x9[i-1],y9[i-1],v6[i],w6[i-1],sep2432.exceed,thresh[69]) +loglik(x10[i-1],y10[i-1],v6[i],canepsilon.xi6[i],oct2432.exceed,thresh[70]) -loglik(x10[i-1],y10[i-1],v6[i],w6[i-1],oct2432.exceed,thresh[70]) +loglik(x11[i-1],y11[i-1],v6[i],canepsilon.xi6[i],nov2432.exceed,thresh[71]) -loglik(x11[i-1],y11[i-1],v6[i],w6[i-1],nov2432.exceed,thresh[71]) +loglik(x12[i-1],y12[i-1],v6[i],canepsilon.xi6[i],dec2432.exceed,thresh[72]) -loglik(x12[i-1],y12[i-1],v6[i],w6[i-1],dec2432.exceed,thresh[72])) aprobepsilon.xi6[i]<-min(1,(clikely6*((dnorm(canepsilon.xi6[i],a.xi6,1/(sqrt(zeta.xi7)))/(dnorm(w6[i-1],a.xi6/(sqrt(zeta.xi7)))))))) if(u=aprobepsilon.xi6[i]){ w6[i]<-w6[i-1] } #eta.one.one[i]<-x1[i]+v1[i] #xi.one.one[i]<-y1[i]+w1[i] #***SIXTH UPDATE DONE!***# a.eta7<-rnorm(1,(c*b+(zeta.eta7*(v1[i]+v2[i]+v3[i]+v4[i]+v5[i]+v6[i]+v7[i-1]+v8[i-1]+v9[i-1]+v10[i-1]+v11[i-1]+v12[i-1])))/(c+12*zeta.eta7),1/(sqrt(c+12*zeta.eta7))) a.xi7<-rnorm(1,(c*b+(zeta.xi7*(w1[i]+w2[i]+w3[i]+w4[i]+w5[i]+w6[i]+w7[i-1]+w8[i-1]+w9[i-1]+w10[i-1]+w11[i-1]+w12[i-1])))/(c+12*zeta.xi7),1/(sqrt(c+12*zeta.xi7))) zeta.eta8<-rgamma(1,f+6,1/(g+0.5*(((v1[i]-a.eta7)^2)+((v2[i]-a.eta7)^2)+((v3[i]-a.eta7)^2)+((v4[i]-a.eta7)^2)+((v5[i]-a.eta7)^2)+((v6[i]-a.eta7)^2)+((v7[i-1]-a.eta7)^2)+((v8[i-1]-a.eta7)^2)+((v9[i-1]-a.eta7)^2)+((v10[i-1]-a.eta7)^2)+((v11[i-1]-a.eta7)^2)+((v12[i-1]-a.eta7)^2)))) zeta.xi8<-rgamma(1,f+6,1/(g+0.5*(((w1[i]-a.xi7)^2)+((w2[i]-a.xi7)^2)+((w3[i]-a.xi7)^2)+((w4[i]-a.xi7)^2)+((w5[i]-a.xi7)^2)+((w6[i]-a.xi7)^2)+((w7[i-1]-a.xi7)^2)+((w8[i-1]-a.xi7)^2)+((w9[i-1]-a.xi7)^2)+((w10[i-1]-a.xi7)^2)+((w11[i-1]-a.xi7)^2)+((w12[i-1]-a.xi7)^2)))) tau.eta8<-rgamma(1,d+6,1/(e+0.5*(((x1[i])^2)+((x2[i])^2)+((x3[i])^2)+((x4[i])^2)+((x5[i])^2)+((x6[i])^2)+((x7[i-1])^2)+((x8[i-1])^2)+((x9[i-1])^2)+((x10[i-1])^2)+((x11[i-1])^2)+((x12[i-1])^2)))) tau.xi8<-rgamma(1,d+6,1/(e+0.5*(((y1[i])^2)+((y2[i])^2)+((y3[i])^2)+((y4[i])^2)+((y5[i])^2)+((y6[i])^2)+((y7[i-1])^2)+((y8[i-1])^2)+((y9[i-1])^2)+((y10[i-1])^2)+((y11[i-1])^2)+((y12[i-1])^2)))) cangamma.eta7[i]<-x7[i-1]+rnorm(1,0,errgamma.eta7) likely7<-exp(loglik(cangamma.eta7[i],y7[i-1],v1[i],w1[i],julshef.exceed,thresh[7]) -loglik(x7[i-1],y7[i-1],v1[i],w1[i],julshef.exceed,thresh[7]) +loglik(cangamma.eta7[i],y7[i-1],v2[i],w2[i],julbrad.exceed,thresh[19]) -loglik(x7[i-1],y7[i-1],v2[i],w2[i],julbrad.exceed,thresh[19]) +loglik(cangamma.eta7[i],y7[i-1],v3[i],w3[i],jul2245.exceed,thresh[31]) -loglik(x7[i-1],y7[i-1],v3[i],w3[i],jul2245.exceed,thresh[31]) +loglik(cangamma.eta7[i],y7[i-1],v4[i],w4[i],jul2255.exceed,thresh[43]) -loglik(x7[i-1],y7[i-1],v4[i-1],w4[i-1],jul2255.exceed,thresh[43]) +loglik(cangamma.eta7[i],y7[i-1],v5[i],w5[i],jul2416.exceed,thresh[55]) -loglik(x7[i-1],y7[i-1],v5[i],w5[i],jul2416.exceed,thresh[55]) +loglik(cangamma.eta7[i],y7[i-1],v6[i],w6[i],jul2432.exceed,thresh[67]) -loglik(x7[i-1],y7[i-1],v6[i],w6[i],jul2432.exceed,thresh[67]) +loglik(cangamma.eta7[i],y7[i-1],v7[i-1],w7[i-1],jul4036.exceed,thresh[79]) -loglik(x7[i-1],y7[i-1],v7[i-1],w7[i-1],jul4036.exceed,thresh[79]) +loglik(cangamma.eta7[i],y7[i-1],v8[i-1],w8[i-1],jul4206.exceed,thresh[91]) -loglik(x7[i-1],y7[i-1],v8[i-1],w8[i-1],jul4206.exceed,thresh[91]) +loglik(cangamma.eta7[i],y7[i-1],v9[i-1],w9[i-1],jul4238.exceed,thresh[103]) -loglik(x7[i-1],y7[i-1],v9[i-1],w9[i-1],jul4238.exceed,thresh[103]) +loglik(cangamma.eta7[i],y7[i-1],v10[i-1],w10[i-1],jul4617.exceed,thresh[115]) -loglik(x7[i-1],y7[i-1],v10[i-1],w10[i-1],jul4617.exceed,thresh[115]) +loglik(cangamma.eta7[i],y7[i-1],v11[i-1],w11[i-1],jul7221.exceed,thresh[127]) -loglik(x7[i-1],y7[i-1],v11[i-1],w11[i-1],jul7221.exceed,thresh[127]) +loglik(cangamma.eta7[i],y7[i-1],v12[i-1],w12[i-1],jul7377.exceed,thresh[139]) -loglik(x7[i-1],y7[i-1],v12[i-1],w12[i-1],jul7377.exceed,thresh[139])) aprobgamma.eta7[i]<-min(1,(likely7*((dnorm(cangamma.eta7[i],0,1/(sqrt(tau.eta8)))/(dnorm(x7[i-1],0,1/(sqrt(tau.eta8)))))))) u<-runif(1) if(u=aprobgamma.eta7[i]){ x7[i]<-x7[i-1] } cangamma.xi7[i]<-y7[i-1]+rnorm(1,0,errgamma.xi7) alikely7<-exp(loglik(x7[i],cangamma.xi7[i],v1[i],w1[i],julshef.exceed,thresh[7]) -loglik(x7[i],y7[i-1],v1[i],w1[i],julshef.exceed,thresh[7]) +loglik(x7[i],cangamma.xi7[i],v2[i],w2[i],julbrad.exceed,thresh[19]) -loglik(x7[i],y7[i-1],v2[i],w2[i],julbrad.exceed,thresh[19]) +loglik(x7[i],cangamma.xi7[i],v3[i],w3[i],jul2245.exceed,thresh[31]) -loglik(x7[i],y7[i-1],v3[i],w3[i],jul2245.exceed,thresh[31]) +loglik(x7[i],cangamma.xi7[i],v4[i],w4[i],jul2255.exceed,thresh[43]) -loglik(x7[i],y7[i-1],v4[i],w4[i],jul2255.exceed,thresh[43]) +loglik(x7[i],cangamma.xi7[i],v5[i],w5[i],jul2416.exceed,thresh[55]) -loglik(x7[i],y7[i-1],v5[i],w5[i],jul2416.exceed,thresh[55]) +loglik(x7[i],cangamma.xi7[i],v6[i],w6[i],jul2432.exceed,thresh[67]) -loglik(x7[i],y7[i-1],v6[i],w6[i],jul2432.exceed,thresh[67]) +loglik(x7[i],cangamma.xi7[i],v7[i-1],w7[i-1],jul4036.exceed,thresh[79]) -loglik(x7[i],y7[i-1],v7[i-1],w7[i-1],jul4036.exceed,thresh[79]) +loglik(x7[i],cangamma.xi7[i],v8[i-1],w8[i-1],jul4206.exceed,thresh[91]) -loglik(x7[i],y7[i-1],v8[i-1],w8[i-1],jul4206.exceed,thresh[91]) +loglik(x7[i],cangamma.xi7[i],v9[i-1],w9[i-1],jul4238.exceed,thresh[103]) -loglik(x7[i],y7[i-1],v9[i-1],w9[i-1],jul4238.exceed,thresh[103]) +loglik(x7[i],cangamma.xi7[i],v10[i-1],w10[i-1],jul4617.exceed,thresh[115]) -loglik(x7[i],y7[i-1],v10[i-1],w10[i-1],jul4617.exceed,thresh[115]) +loglik(x7[i],cangamma.xi7[i],v11[i-1],w11[i-1],jul7221.exceed,thresh[127]) -loglik(x7[i],y7[i-1],v11[i-1],w11[i-1],jul7221.exceed,thresh[127]) +loglik(x7[i],cangamma.xi7[i],v12[i-1],w12[i-1],jul7377.exceed,thresh[139]) -loglik(x7[i],y7[i-1],v12[i-1],w12[i-1],jul7377.exceed,thresh[139])) aprobgamma.xi7[i]<-min(1,(alikely7*((dnorm(cangamma.xi7[i],0,1/(sqrt(tau.xi8)))/(dnorm(y7[i-1],0,1/(sqrt(tau.xi8)))))))) if(u=aprobgamma.xi7[i]){ y7[i]<-y7[i-1] } canepsilon.eta7[i]<-v7[i-1]+rnorm(1,0,errepsilon.eta7) blikely7<-exp(loglik(x1[i],y1[i],canepsilon.eta7[i],w7[i-1],jan4036.exceed,thresh[73]) -loglik(x1[i],y1[i],v7[i-1],w7[i-1],jan4036.exceed,thresh[73]) +loglik(x2[i],y2[i],canepsilon.eta7[i],w7[i-1],feb4036.exceed,thresh[74]) -loglik(x2[i],y2[i],v7[i-1],w7[i-1],feb4036.exceed,thresh[74]) +loglik(x3[i],y3[i],canepsilon.eta7[i],w7[i-1],mar4036.exceed,thresh[75]) -loglik(x3[i],y3[i],v7[i-1],w7[i-1],mar4036.exceed,thresh[75]) +loglik(x4[i],y4[i],canepsilon.eta7[i],w7[i-1],apr4036.exceed,thresh[76]) -loglik(x4[i],y4[i],v7[i-1],w7[i-1],apr4036.exceed,thresh[76]) +loglik(x5[i],y5[i],canepsilon.eta7[i],w7[i-1],may4036.exceed,thresh[77]) -loglik(x5[i],y5[i],v7[i-1],w7[i-1],may4036.exceed,thresh[77]) +loglik(x6[i],y6[i],canepsilon.eta7[i],w7[i-1],jun4036.exceed,thresh[78]) -loglik(x6[i],y6[i],v7[i-1],w7[i-1],jun4036.exceed,thresh[78]) +loglik(x7[i],y7[i],canepsilon.eta7[i],w7[i-1],jul4036.exceed,thresh[79]) -loglik(x7[i],y7[i],v7[i-1],w7[i-1],jul4036.exceed,thresh[79]) +loglik(x8[i-1],y8[i-1],canepsilon.eta7[i],w7[i-1],aug4036.exceed,thresh[80]) -loglik(x8[i-1],y8[i-1],v7[i-1],w7[i-1],aug4036.exceed,thresh[80]) +loglik(x9[i-1],y9[i-1],canepsilon.eta7[i],w7[i-1],sep4036.exceed,thresh[81]) -loglik(x9[i-1],y9[i-1],v7[i-1],w7[i-1],sep4036.exceed,thresh[81]) +loglik(x10[i-1],y10[i-1],canepsilon.eta7[i],w7[i-1],oct4036.exceed,thresh[82]) -loglik(x10[i-1],y10[i-1],v7[i-1],w7[i-1],oct4036.exceed,thresh[82]) +loglik(x11[i-1],y11[i-1],canepsilon.eta7[i],w7[i-1],nov4036.exceed,thresh[83]) -loglik(x11[i-1],y11[i-1],v7[i-1],w7[i-1],nov4036.exceed,thresh[83]) +loglik(x12[i-1],y12[i-1],canepsilon.eta7[i],w7[i-1],dec4036.exceed,thresh[84]) -loglik(x12[i-1],y12[i-1],v7[i-1],w7[i-1],dec4036.exceed,thresh[84])) aprobepsilon.eta7[i]<-min(1,(blikely7*((dnorm(canepsilon.eta7[i],a.eta7,1/(sqrt(zeta.eta8)))/(dnorm(v7[i-1],a.eta7,1/(sqrt(zeta.eta8)))))))) if(u=aprobepsilon.eta7[i]){ v7[i]<-v7[i-1] } canepsilon.xi7[i]<-w7[i-1]+rnorm(1,0,errepsilon.xi7) clikely7<-exp(loglik(x1[i],y1[i],v7[i],canepsilon.xi7[i],jan4036.exceed,thresh[73]) -loglik(x1[i],y1[i],v7[i],w7[i-1],jan4036.exceed,thresh[73]) +loglik(x2[i],y2[i],v7[i],canepsilon.xi7[i],feb4036.exceed,thresh[74]) -loglik(x2[i],y2[i],v7[i],w7[i-1],feb4036.exceed,thresh[74]) +loglik(x3[i],y3[i],v7[i],canepsilon.xi7[i],mar4036.exceed,thresh[75]) -loglik(x3[i],y3[i],v7[i],w7[i-1],mar4036.exceed,thresh[75]) +loglik(x4[i],y4[i],v7[i],canepsilon.xi7[i],apr4036.exceed,thresh[76]) -loglik(x4[i],y4[i],v7[i],w7[i-1],apr4036.exceed,thresh[76]) +loglik(x5[i],y5[i],v7[i],canepsilon.xi7[i],may4036.exceed,thresh[77]) -loglik(x5[i],y5[i],v7[i],w7[i-1],may4036.exceed,thresh[77]) +loglik(x6[i],y6[i],v7[i],canepsilon.xi7[i],jun4036.exceed,thresh[78]) -loglik(x6[i],y6[i],v7[i],w7[i-1],jun4036.exceed,thresh[78]) +loglik(x7[i],y7[i],v7[i],canepsilon.xi7[i],jul4036.exceed,thresh[79]) -loglik(x7[i],y7[i],v7[i],w7[i-1],jul4036.exceed,thresh[79]) +loglik(x8[i-1],y8[i-1],v7[i],canepsilon.xi7[i],aug4036.exceed,thresh[80]) -loglik(x8[i-1],y8[i-1],v7[i],w7[i-1],aug4036.exceed,thresh[80]) +loglik(x9[i-1],y9[i-1],v7[i],canepsilon.xi7[i],sep4036.exceed,thresh[81]) -loglik(x9[i-1],y9[i-1],v7[i],w7[i-1],sep4036.exceed,thresh[81]) +loglik(x10[i-1],y10[i-1],v7[i],canepsilon.xi7[i],oct4036.exceed,thresh[82]) -loglik(x10[i-1],y10[i-1],v7[i],w7[i-1],oct4036.exceed,thresh[82]) +loglik(x11[i-1],y11[i-1],v7[i],canepsilon.xi7[i],nov4036.exceed,thresh[83]) -loglik(x11[i-1],y11[i-1],v7[i],w7[i-1],nov4036.exceed,thresh[83]) +loglik(x12[i-1],y12[i-1],v7[i],canepsilon.xi7[i],dec4036.exceed,thresh[84]) -loglik(x12[i-1],y12[i-1],v7[i],w7[i-1],dec4036.exceed,thresh[84])) aprobepsilon.xi7[i]<-min(1,(clikely7*((dnorm(canepsilon.xi7[i],a.xi7,1/(sqrt(zeta.xi8)))/(dnorm(w7[i-1],a.xi7/(sqrt(zeta.xi8)))))))) if(u=aprobepsilon.xi7[i]){ w7[i]<-w7[i-1] } #eta.one.one[i]<-x1[i]+v1[i] #xi.one.one[i]<-y1[i]+w1[i] #***SEVENTH UPDATE DONE!***# a.eta8<-rnorm(1,(c*b+(zeta.eta8*(v1[i]+v2[i]+v3[i]+v4[i]+v5[i]+v6[i]+v7[i]+v8[i-1]+v9[i-1]+v10[i-1]+v11[i-1]+v12[i-1])))/(c+12*zeta.eta8),1/(sqrt(c+12*zeta.eta8))) a.xi8<-rnorm(1,(c*b+(zeta.xi8*(w1[i]+w2[i]+w3[i]+w4[i]+w5[i]+w6[i]+w7[i]+w8[i-1]+w9[i-1]+w10[i-1]+w11[i-1]+w12[i-1])))/(c+12*zeta.xi8),1/(sqrt(c+12*zeta.xi8))) zeta.eta9<-rgamma(1,f+6,1/(g+0.5*(((v1[i]-a.eta8)^2)+((v2[i]-a.eta8)^2)+((v3[i]-a.eta8)^2)+((v4[i]-a.eta8)^2)+((v5[i]-a.eta8)^2)+((v6[i]-a.eta8)^2)+((v7[i]-a.eta8)^2)+((v8[i-1]-a.eta8)^2)+((v9[i-1]-a.eta8)^2)+((v10[i-1]-a.eta8)^2)+((v11[i-1]-a.eta8)^2)+((v12[i-1]-a.eta8)^2)))) zeta.xi9<-rgamma(1,f+6,1/(g+0.5*(((w1[i]-a.xi8)^2)+((w2[i]-a.xi8)^2)+((w3[i]-a.xi8)^2)+((w4[i]-a.xi8)^2)+((w5[i]-a.xi8)^2)+((w6[i]-a.xi8)^2)+((w7[i]-a.xi8)^2)+((w8[i-1]-a.xi8)^2)+((w9[i-1]-a.xi8)^2)+((w10[i-1]-a.xi8)^2)+((w11[i-1]-a.xi8)^2)+((w12[i-1]-a.xi8)^2)))) tau.eta9<-rgamma(1,d+6,1/(e+0.5*(((x1[i])^2)+((x2[i])^2)+((x3[i])^2)+((x4[i])^2)+((x5[i])^2)+((x6[i])^2)+((x7[i])^2)+((x8[i-1])^2)+((x9[i-1])^2)+((x10[i-1])^2)+((x11[i-1])^2)+((x12[i-1])^2)))) tau.xi9<-rgamma(1,d+6,1/(e+0.5*(((y1[i])^2)+((y2[i])^2)+((y3[i])^2)+((y4[i])^2)+((y5[i])^2)+((y6[i])^2)+((y7[i])^2)+((y8[i-1])^2)+((y9[i-1])^2)+((y10[i-1])^2)+((y11[i-1])^2)+((y12[i-1])^2)))) cangamma.eta8[i]<-x8[i-1]+rnorm(1,0,errgamma.eta8) likely8<-exp(loglik(cangamma.eta8[i],y8[i-1],v1[i],w1[i],augshef.exceed,thresh[8]) -loglik(x8[i-1],y8[i-1],v1[i],w1[i],augshef.exceed,thresh[8]) +loglik(cangamma.eta8[i],y8[i-1],v2[i],w2[i],augbrad.exceed,thresh[20]) -loglik(x8[i-1],y8[i-1],v2[i],w2[i],augbrad.exceed,thresh[20]) +loglik(cangamma.eta8[i],y8[i-1],v3[i],w3[i],aug2245.exceed,thresh[32]) -loglik(x8[i-1],y8[i-1],v3[i],w3[i],aug2245.exceed,thresh[32]) +loglik(cangamma.eta8[i],y8[i-1],v4[i],w4[i],aug2255.exceed,thresh[44]) -loglik(x8[i-1],y8[i-1],v4[i-1],w4[i-1],aug2255.exceed,thresh[44]) +loglik(cangamma.eta8[i],y8[i-1],v5[i],w5[i],aug2416.exceed,thresh[56]) -loglik(x8[i-1],y8[i-1],v5[i],w5[i],aug2416.exceed,thresh[56]) +loglik(cangamma.eta8[i],y8[i-1],v6[i],w6[i],aug2432.exceed,thresh[68]) -loglik(x8[i-1],y8[i-1],v6[i],w6[i],aug2432.exceed,thresh[68]) +loglik(cangamma.eta8[i],y8[i-1],v7[i],w7[i],aug4036.exceed,thresh[80]) -loglik(x8[i-1],y8[i-1],v7[i],w7[i],aug4036.exceed,thresh[80]) +loglik(cangamma.eta8[i],y8[i-1],v8[i-1],w8[i-1],aug4206.exceed,thresh[92]) -loglik(x8[i-1],y8[i-1],v8[i-1],w8[i-1],aug4206.exceed,thresh[92]) +loglik(cangamma.eta8[i],y8[i-1],v9[i-1],w9[i-1],aug4238.exceed,thresh[104]) -loglik(x8[i-1],y8[i-1],v9[i-1],w9[i-1],aug4238.exceed,thresh[104]) +loglik(cangamma.eta8[i],y8[i-1],v10[i-1],w10[i-1],aug4617.exceed,thresh[116]) -loglik(x8[i-1],y8[i-1],v10[i-1],w10[i-1],aug4617.exceed,thresh[116]) +loglik(cangamma.eta8[i],y8[i-1],v11[i-1],w11[i-1],aug7221.exceed,thresh[128]) -loglik(x8[i-1],y8[i-1],v11[i-1],w11[i-1],aug7221.exceed,thresh[128]) +loglik(cangamma.eta8[i],y8[i-1],v12[i-1],w12[i-1],aug7377.exceed,thresh[140]) -loglik(x8[i-1],y8[i-1],v12[i-1],w12[i-1],aug7377.exceed,thresh[140])) aprobgamma.eta8[i]<-min(1,(likely8*((dnorm(cangamma.eta8[i],0,1/(sqrt(tau.eta9)))/(dnorm(x8[i-1],0,1/(sqrt(tau.eta9)))))))) u<-runif(1) if(u=aprobgamma.eta8[i]){ x8[i]<-x8[i-1] } cangamma.xi8[i]<-y8[i-1]+rnorm(1,0,errgamma.xi8) alikely8<-exp(loglik(x8[i],cangamma.xi8[i],v1[i],w1[i],augshef.exceed,thresh[8]) -loglik(x8[i],y8[i-1],v1[i],w1[i],augshef.exceed,thresh[8]) +loglik(x8[i],cangamma.xi8[i],v2[i],w2[i],augbrad.exceed,thresh[20]) -loglik(x8[i],y8[i-1],v2[i],w2[i],augbrad.exceed,thresh[20]) +loglik(x8[i],cangamma.xi8[i],v3[i],w3[i],aug2245.exceed,thresh[32]) -loglik(x8[i],y8[i-1],v3[i],w3[i],aug2245.exceed,thresh[32]) +loglik(x8[i],cangamma.xi8[i],v4[i],w4[i],aug2255.exceed,thresh[44]) -loglik(x8[i],y8[i-1],v4[i],w4[i],aug2255.exceed,thresh[44]) +loglik(x8[i],cangamma.xi8[i],v5[i],w5[i],aug2416.exceed,thresh[56]) -loglik(x8[i],y8[i-1],v5[i],w5[i],aug2416.exceed,thresh[56]) +loglik(x8[i],cangamma.xi8[i],v6[i],w6[i],aug2432.exceed,thresh[68]) -loglik(x8[i],y8[i-1],v6[i],w6[i],aug2432.exceed,thresh[68]) +loglik(x8[i],cangamma.xi8[i],v7[i],w7[i],aug4036.exceed,thresh[80]) -loglik(x8[i],y8[i-1],v7[i],w7[i],aug4036.exceed,thresh[80]) +loglik(x8[i],cangamma.xi8[i],v8[i-1],w8[i-1],aug4206.exceed,thresh[92]) -loglik(x8[i],y8[i-1],v8[i-1],w8[i-1],aug4206.exceed,thresh[92]) +loglik(x8[i],cangamma.xi8[i],v9[i-1],w9[i-1],aug4238.exceed,thresh[104]) -loglik(x8[i],y8[i-1],v9[i-1],w9[i-1],aug4238.exceed,thresh[104]) +loglik(x8[i],cangamma.xi8[i],v10[i-1],w10[i-1],aug4617.exceed,thresh[116]) -loglik(x8[i],y8[i-1],v10[i-1],w10[i-1],aug4617.exceed,thresh[116]) +loglik(x8[i],cangamma.xi8[i],v11[i-1],w11[i-1],aug7221.exceed,thresh[128]) -loglik(x8[i],y8[i-1],v11[i-1],w11[i-1],aug7221.exceed,thresh[128]) +loglik(x8[i],cangamma.xi8[i],v12[i-1],w12[i-1],aug7377.exceed,thresh[140]) -loglik(x8[i],y8[i-1],v12[i-1],w12[i-1],aug7377.exceed,thresh[140])) aprobgamma.xi8[i]<-min(1,(alikely8*((dnorm(cangamma.xi8[i],0,1/(sqrt(tau.xi9)))/(dnorm(y8[i-1],0,1/(sqrt(tau.xi9)))))))) if(u=aprobgamma.xi8[i]){ y8[i]<-y8[i-1] } canepsilon.eta8[i]<-v8[i-1]+rnorm(1,0,errepsilon.eta8) blikely8<-exp(loglik(x1[i],y1[i],canepsilon.eta8[i],w8[i-1],jan4206.exceed,thresh[85]) -loglik(x1[i],y1[i],v8[i-1],w8[i-1],jan4206.exceed,thresh[85]) +loglik(x2[i],y2[i],canepsilon.eta8[i],w8[i-1],feb4206.exceed,thresh[86]) -loglik(x2[i],y2[i],v8[i-1],w8[i-1],feb4206.exceed,thresh[86]) +loglik(x3[i],y3[i],canepsilon.eta8[i],w8[i-1],mar4206.exceed,thresh[87]) -loglik(x3[i],y3[i],v8[i-1],w8[i-1],mar4206.exceed,thresh[87]) +loglik(x4[i],y4[i],canepsilon.eta8[i],w8[i-1],apr4206.exceed,thresh[88]) -loglik(x4[i],y4[i],v8[i-1],w8[i-1],apr4206.exceed,thresh[88]) +loglik(x5[i],y5[i],canepsilon.eta8[i],w8[i-1],may4206.exceed,thresh[89]) -loglik(x5[i],y5[i],v8[i-1],w8[i-1],may4206.exceed,thresh[89]) +loglik(x6[i],y6[i],canepsilon.eta8[i],w8[i-1],jun4206.exceed,thresh[90]) -loglik(x6[i],y6[i],v8[i-1],w8[i-1],jun4206.exceed,thresh[90]) +loglik(x7[i],y7[i],canepsilon.eta8[i],w8[i-1],jul4206.exceed,thresh[91]) -loglik(x7[i],y7[i],v8[i-1],w8[i-1],jul4206.exceed,thresh[91]) +loglik(x8[i],y8[i],canepsilon.eta8[i],w8[i-1],aug4206.exceed,thresh[92]) -loglik(x8[i],y8[i],v8[i-1],w8[i-1],aug4206.exceed,thresh[92]) +loglik(x9[i-1],y9[i-1],canepsilon.eta8[i],w8[i-1],sep4206.exceed,thresh[93]) -loglik(x9[i-1],y9[i-1],v8[i-1],w8[i-1],sep4206.exceed,thresh[93]) +loglik(x10[i-1],y10[i-1],canepsilon.eta8[i],w8[i-1],oct4206.exceed,thresh[94]) -loglik(x10[i-1],y10[i-1],v8[i-1],w8[i-1],oct4206.exceed,thresh[94]) +loglik(x11[i-1],y11[i-1],canepsilon.eta8[i],w8[i-1],nov4206.exceed,thresh[95]) -loglik(x11[i-1],y11[i-1],v8[i-1],w8[i-1],nov4206.exceed,thresh[95]) +loglik(x12[i-1],y12[i-1],canepsilon.eta8[i],w8[i-1],dec4206.exceed,thresh[96]) -loglik(x12[i-1],y12[i-1],v8[i-1],w8[i-1],dec4206.exceed,thresh[96])) aprobepsilon.eta8[i]<-min(1,(blikely8*((dnorm(canepsilon.eta8[i],a.eta8,1/(sqrt(zeta.eta9)))/(dnorm(v8[i-1],a.eta8,1/(sqrt(zeta.eta9)))))))) if(u=aprobepsilon.eta8[i]){ v8[i]<-v8[i-1] } canepsilon.xi8[i]<-w8[i-1]+rnorm(1,0,errepsilon.xi8) clikely8<-exp(loglik(x1[i],y1[i],v8[i],canepsilon.xi8[i],jan4206.exceed,thresh[85]) -loglik(x1[i],y1[i],v8[i],w8[i-1],jan4206.exceed,thresh[85]) +loglik(x2[i],y2[i],v8[i],canepsilon.xi8[i],feb4206.exceed,thresh[86]) -loglik(x2[i],y2[i],v8[i],w8[i-1],feb4206.exceed,thresh[86]) +loglik(x3[i],y3[i],v8[i],canepsilon.xi8[i],mar4206.exceed,thresh[87]) -loglik(x3[i],y3[i],v8[i],w8[i-1],mar4206.exceed,thresh[87]) +loglik(x4[i],y4[i],v8[i],canepsilon.xi8[i],apr4206.exceed,thresh[88]) -loglik(x4[i],y4[i],v8[i],w8[i-1],apr4206.exceed,thresh[88]) +loglik(x5[i],y5[i],v8[i],canepsilon.xi8[i],may4206.exceed,thresh[89]) -loglik(x5[i],y5[i],v8[i],w8[i-1],may4206.exceed,thresh[89]) +loglik(x6[i],y6[i],v8[i],canepsilon.xi8[i],jun4206.exceed,thresh[90]) -loglik(x6[i],y6[i],v8[i],w8[i-1],jun4206.exceed,thresh[90]) +loglik(x7[i],y7[i],v8[i],canepsilon.xi8[i],jul4206.exceed,thresh[91]) -loglik(x7[i],y7[i],v8[i],w8[i-1],jul4206.exceed,thresh[91]) +loglik(x8[i],y8[i],v8[i],canepsilon.xi8[i],aug4206.exceed,thresh[92]) -loglik(x8[i],y8[i],v8[i],w8[i-1],aug4206.exceed,thresh[92]) +loglik(x9[i-1],y9[i-1],v8[i],canepsilon.xi8[i],sep4206.exceed,thresh[93]) -loglik(x9[i-1],y9[i-1],v8[i],w8[i-1],sep4206.exceed,thresh[93]) +loglik(x10[i-1],y10[i-1],v8[i],canepsilon.xi8[i],oct4206.exceed,thresh[94]) -loglik(x10[i-1],y10[i-1],v8[i],w8[i-1],oct4206.exceed,thresh[94]) +loglik(x11[i-1],y11[i-1],v8[i],canepsilon.xi8[i],nov4206.exceed,thresh[95]) -loglik(x11[i-1],y11[i-1],v8[i],w8[i-1],nov4206.exceed,thresh[95]) +loglik(x12[i-1],y12[i-1],v8[i],canepsilon.xi8[i],dec4206.exceed,thresh[96]) -loglik(x12[i-1],y12[i-1],v8[i],w8[i-1],dec4206.exceed,thresh[96])) aprobepsilon.xi8[i]<-min(1,(clikely8*((dnorm(canepsilon.xi8[i],a.xi8,1/(sqrt(zeta.xi9)))/(dnorm(w8[i-1],a.xi8/(sqrt(zeta.xi9)))))))) if(u=aprobepsilon.xi8[i]){ w8[i]<-w8[i-1] } #eta.one.one[i]<-x1[i]+v1[i] #xi.one.one[i]<-y1[i]+w1[i] #***EIGHTH UPDATE DONE!***# a.eta9<-rnorm(1,(c*b+(zeta.eta9*(v1[i]+v2[i]+v3[i]+v4[i]+v5[i]+v6[i]+v7[i]+v8[i]+v9[i-1]+v10[i-1]+v11[i-1]+v12[i-1])))/(c+12*zeta.eta9),1/(sqrt(c+12*zeta.eta9))) a.xi9<-rnorm(1,(c*b+(zeta.xi9*(w1[i]+w2[i]+w3[i]+w4[i]+w5[i]+w6[i]+w7[i]+w8[i]+w9[i-1]+w10[i-1]+w11[i-1]+w12[i-1])))/(c+12*zeta.xi9),1/(sqrt(c+12*zeta.xi9))) zeta.eta10<-rgamma(1,f+6,1/(g+0.5*(((v1[i]-a.eta9)^2)+((v2[i]-a.eta9)^2)+((v3[i]-a.eta9)^2)+((v4[i]-a.eta9)^2)+((v5[i]-a.eta9)^2)+((v6[i]-a.eta9)^2)+((v7[i]-a.eta9)^2)+((v8[i]-a.eta9)^2)+((v9[i-1]-a.eta9)^2)+((v10[i-1]-a.eta9)^2)+((v11[i-1]-a.eta9)^2)+((v12[i-1]-a.eta9)^2)))) zeta.xi10<-rgamma(1,f+6,1/(g+0.5*(((w1[i]-a.xi9)^2)+((w2[i]-a.xi9)^2)+((w3[i]-a.xi9)^2)+((w4[i]-a.xi9)^2)+((w5[i]-a.xi9)^2)+((w6[i]-a.xi9)^2)+((w7[i]-a.xi9)^2)+((w8[i]-a.xi9)^2)+((w9[i-1]-a.xi9)^2)+((w10[i-1]-a.xi9)^2)+((w11[i-1]-a.xi9)^2)+((w12[i-1]-a.xi9)^2)))) tau.eta10<-rgamma(1,d+6,1/(e+0.5*(((x1[i])^2)+((x2[i])^2)+((x3[i])^2)+((x4[i])^2)+((x5[i])^2)+((x6[i])^2)+((x7[i])^2)+((x8[i])^2)+((x9[i-1])^2)+((x10[i-1])^2)+((x11[i-1])^2)+((x12[i-1])^2)))) tau.xi10<-rgamma(1,d+6,1/(e+0.5*(((y1[i])^2)+((y2[i])^2)+((y3[i])^2)+((y4[i])^2)+((y5[i])^2)+((y6[i])^2)+((y7[i])^2)+((y8[i])^2)+((y9[i-1])^2)+((y10[i-1])^2)+((y11[i-1])^2)+((y12[i-1])^2)))) cangamma.eta9[i]<-x9[i-1]+rnorm(1,0,errgamma.eta9) likely9<-exp(loglik(cangamma.eta9[i],y9[i-1],v1[i],w1[i],sepshef.exceed,thresh[9]) -loglik(x9[i-1],y9[i-1],v1[i],w1[i],sepshef.exceed,thresh[9]) +loglik(cangamma.eta9[i],y9[i-1],v2[i],w2[i],sepbrad.exceed,thresh[21]) -loglik(x9[i-1],y9[i-1],v2[i],w2[i],sepbrad.exceed,thresh[21]) +loglik(cangamma.eta9[i],y9[i-1],v3[i],w3[i],sep2245.exceed,thresh[33]) -loglik(x9[i-1],y9[i-1],v3[i],w3[i],sep2245.exceed,thresh[33]) +loglik(cangamma.eta9[i],y9[i-1],v4[i],w4[i],sep2255.exceed,thresh[45]) -loglik(x9[i-1],y9[i-1],v4[i-1],w4[i-1],sep2255.exceed,thresh[45]) +loglik(cangamma.eta9[i],y9[i-1],v5[i],w5[i],sep2416.exceed,thresh[57]) -loglik(x9[i-1],y9[i-1],v5[i],w5[i],sep2416.exceed,thresh[57]) +loglik(cangamma.eta9[i],y9[i-1],v6[i],w6[i],sep2432.exceed,thresh[69]) -loglik(x9[i-1],y9[i-1],v6[i],w6[i],sep2432.exceed,thresh[69]) +loglik(cangamma.eta9[i],y9[i-1],v7[i],w7[i],sep4036.exceed,thresh[81]) -loglik(x9[i-1],y9[i-1],v7[i],w7[i],sep4036.exceed,thresh[81]) +loglik(cangamma.eta9[i],y9[i-1],v8[i],w8[i],sep4206.exceed,thresh[93]) -loglik(x9[i-1],y9[i-1],v8[i],w8[i],sep4206.exceed,thresh[93]) +loglik(cangamma.eta9[i],y9[i-1],v9[i-1],w9[i-1],sep4238.exceed,thresh[105]) -loglik(x9[i-1],y9[i-1],v9[i-1],w9[i-1],sep4238.exceed,thresh[105]) +loglik(cangamma.eta9[i],y9[i-1],v10[i-1],w10[i-1],sep4617.exceed,thresh[117]) -loglik(x9[i-1],y9[i-1],v10[i-1],w10[i-1],sep4617.exceed,thresh[117]) +loglik(cangamma.eta9[i],y9[i-1],v11[i-1],w11[i-1],sep7221.exceed,thresh[129]) -loglik(x9[i-1],y9[i-1],v11[i-1],w11[i-1],sep7221.exceed,thresh[129]) +loglik(cangamma.eta9[i],y9[i-1],v12[i-1],w12[i-1],sep7377.exceed,thresh[141]) -loglik(x9[i-1],y9[i-1],v12[i-1],w12[i-1],sep7377.exceed,thresh[141])) aprobgamma.eta9[i]<-min(1,(likely9*((dnorm(cangamma.eta9[i],0,1/(sqrt(tau.eta10)))/(dnorm(x9[i-1],0,1/(sqrt(tau.eta10)))))))) u<-runif(1) if(u=aprobgamma.eta9[i]){ x9[i]<-x9[i-1] } cangamma.xi9[i]<-y9[i-1]+rnorm(1,0,errgamma.xi9) alikely9<-exp(loglik(x9[i],cangamma.xi9[i],v1[i],w1[i],sepshef.exceed,thresh[9]) -loglik(x9[i],y9[i-1],v1[i],w1[i],sepshef.exceed,thresh[9]) +loglik(x9[i],cangamma.xi9[i],v2[i],w2[i],sepbrad.exceed,thresh[21]) -loglik(x9[i],y9[i-1],v2[i],w2[i],sepbrad.exceed,thresh[21]) +loglik(x9[i],cangamma.xi9[i],v3[i],w3[i],sep2245.exceed,thresh[33]) -loglik(x9[i],y9[i-1],v3[i],w3[i],sep2245.exceed,thresh[33]) +loglik(x9[i],cangamma.xi9[i],v4[i],w4[i],sep2255.exceed,thresh[45]) -loglik(x9[i],y9[i-1],v4[i],w4[i],sep2255.exceed,thresh[45]) +loglik(x9[i],cangamma.xi9[i],v5[i],w5[i],sep2416.exceed,thresh[57]) -loglik(x9[i],y9[i-1],v5[i],w5[i],sep2416.exceed,thresh[57]) +loglik(x9[i],cangamma.xi9[i],v6[i],w6[i],sep2432.exceed,thresh[69]) -loglik(x9[i],y9[i-1],v6[i],w6[i],sep2432.exceed,thresh[69]) +loglik(x9[i],cangamma.xi9[i],v7[i],w7[i],sep4036.exceed,thresh[81]) -loglik(x9[i],y9[i-1],v7[i],w7[i],sep4036.exceed,thresh[81]) +loglik(x9[i],cangamma.xi9[i],v8[i],w8[i],sep4206.exceed,thresh[93]) -loglik(x9[i],y9[i-1],v8[i],w8[i],sep4206.exceed,thresh[93]) +loglik(x9[i],cangamma.xi9[i],v9[i-1],w9[i-1],sep4238.exceed,thresh[105]) -loglik(x9[i],y9[i-1],v9[i-1],w9[i-1],sep4238.exceed,thresh[105]) +loglik(x9[i],cangamma.xi9[i],v10[i-1],w10[i-1],sep4617.exceed,thresh[117]) -loglik(x9[i],y9[i-1],v10[i-1],w10[i-1],sep4617.exceed,thresh[117]) +loglik(x9[i],cangamma.xi9[i],v11[i-1],w11[i-1],sep7221.exceed,thresh[129]) -loglik(x9[i],y9[i-1],v11[i-1],w11[i-1],sep7221.exceed,thresh[129]) +loglik(x9[i],cangamma.xi9[i],v12[i-1],w12[i-1],sep7377.exceed,thresh[141]) -loglik(x9[i],y9[i-1],v12[i-1],w12[i-1],sep7377.exceed,thresh[141])) aprobgamma.xi9[i]<-min(1,(alikely9*((dnorm(cangamma.xi9[i],0,1/(sqrt(tau.xi10)))/(dnorm(y9[i-1],0,1/(sqrt(tau.xi10)))))))) if(u=aprobgamma.xi9[i]){ y9[i]<-y9[i-1] } canepsilon.eta9[i]<-v9[i-1]+rnorm(1,0,errepsilon.eta9) blikely9<-exp(loglik(x1[i],y1[i],canepsilon.eta9[i],w9[i-1],jan4238.exceed,thresh[97]) -loglik(x1[i],y1[i],v9[i-1],w9[i-1],jan4238.exceed,thresh[97]) +loglik(x2[i],y2[i],canepsilon.eta9[i],w9[i-1],feb4238.exceed,thresh[98]) -loglik(x2[i],y2[i],v9[i-1],w9[i-1],feb4238.exceed,thresh[98]) +loglik(x3[i],y3[i],canepsilon.eta9[i],w9[i-1],mar4238.exceed,thresh[99]) -loglik(x3[i],y3[i],v9[i-1],w9[i-1],mar4238.exceed,thresh[99]) +loglik(x4[i],y4[i],canepsilon.eta9[i],w9[i-1],apr4238.exceed,thresh[100]) -loglik(x4[i],y4[i],v9[i-1],w9[i-1],apr4238.exceed,thresh[100]) +loglik(x5[i],y5[i],canepsilon.eta9[i],w9[i-1],may4238.exceed,thresh[101]) -loglik(x5[i],y5[i],v9[i-1],w9[i-1],may4238.exceed,thresh[101]) +loglik(x6[i],y6[i],canepsilon.eta9[i],w9[i-1],jun4238.exceed,thresh[102]) -loglik(x6[i],y6[i],v9[i-1],w9[i-1],jun4238.exceed,thresh[102]) +loglik(x7[i],y7[i],canepsilon.eta9[i],w9[i-1],jul4238.exceed,thresh[103]) -loglik(x7[i],y7[i],v9[i-1],w9[i-1],jul4238.exceed,thresh[103]) +loglik(x8[i],y8[i],canepsilon.eta9[i],w9[i-1],aug4238.exceed,thresh[104]) -loglik(x8[i],y8[i],v9[i-1],w9[i-1],aug4238.exceed,thresh[104]) +loglik(x9[i],y9[i],canepsilon.eta9[i],w9[i-1],sep4238.exceed,thresh[105]) -loglik(x9[i],y9[i],v9[i-1],w9[i-1],sep4238.exceed,thresh[105]) +loglik(x10[i-1],y10[i-1],canepsilon.eta9[i],w9[i-1],oct4238.exceed,thresh[106]) -loglik(x10[i-1],y10[i-1],v9[i-1],w9[i-1],oct4238.exceed,thresh[106]) +loglik(x11[i-1],y11[i-1],canepsilon.eta9[i],w9[i-1],nov4238.exceed,thresh[107]) -loglik(x11[i-1],y11[i-1],v9[i-1],w9[i-1],nov4238.exceed,thresh[107]) +loglik(x12[i-1],y12[i-1],canepsilon.eta9[i],w9[i-1],dec4238.exceed,thresh[108]) -loglik(x12[i-1],y12[i-1],v9[i-1],w9[i-1],dec4238.exceed,thresh[108])) aprobepsilon.eta9[i]<-min(1,(blikely9*((dnorm(canepsilon.eta9[i],a.eta9,1/(sqrt(zeta.eta10)))/(dnorm(v9[i-1],a.eta9,1/(sqrt(zeta.eta10)))))))) if(u=aprobepsilon.eta9[i]){ v9[i]<-v9[i-1] } canepsilon.xi9[i]<-w9[i-1]+rnorm(1,0,errepsilon.xi9) clikely9<-exp(loglik(x1[i],y1[i],v9[i],canepsilon.xi9[i],jan4238.exceed,thresh[97]) -loglik(x1[i],y1[i],v9[i],w9[i-1],jan4238.exceed,thresh[97]) +loglik(x2[i],y2[i],v9[i],canepsilon.xi9[i],feb4238.exceed,thresh[98]) -loglik(x2[i],y2[i],v9[i],w9[i-1],feb4238.exceed,thresh[98]) +loglik(x3[i],y3[i],v9[i],canepsilon.xi9[i],mar4238.exceed,thresh[99]) -loglik(x3[i],y3[i],v9[i],w9[i-1],mar4238.exceed,thresh[99]) +loglik(x4[i],y4[i],v9[i],canepsilon.xi9[i],apr4238.exceed,thresh[100]) -loglik(x4[i],y4[i],v9[i],w9[i-1],apr4238.exceed,thresh[100]) +loglik(x5[i],y5[i],v9[i],canepsilon.xi9[i],may4238.exceed,thresh[101]) -loglik(x5[i],y5[i],v9[i],w9[i-1],may4238.exceed,thresh[101]) +loglik(x6[i],y6[i],v9[i],canepsilon.xi9[i],jun4238.exceed,thresh[102]) -loglik(x6[i],y6[i],v9[i],w9[i-1],jun4238.exceed,thresh[102]) +loglik(x7[i],y7[i],v9[i],canepsilon.xi9[i],jul4238.exceed,thresh[103]) -loglik(x7[i],y7[i],v9[i],w9[i-1],jul4238.exceed,thresh[103]) +loglik(x8[i],y8[i],v9[i],canepsilon.xi9[i],aug4238.exceed,thresh[104]) -loglik(x8[i],y8[i],v9[i],w9[i-1],aug4238.exceed,thresh[104]) +loglik(x9[i],y9[i],v9[i],canepsilon.xi9[i],sep4238.exceed,thresh[105]) -loglik(x9[i],y9[i],v9[i],w9[i-1],sep4238.exceed,thresh[105]) +loglik(x10[i-1],y10[i-1],v9[i],canepsilon.xi9[i],oct4238.exceed,thresh[106]) -loglik(x10[i-1],y10[i-1],v9[i],w9[i-1],oct4238.exceed,thresh[106]) +loglik(x11[i-1],y11[i-1],v9[i],canepsilon.xi9[i],nov4238.exceed,thresh[107]) -loglik(x11[i-1],y11[i-1],v9[i],w9[i-1],nov4238.exceed,thresh[107]) +loglik(x12[i-1],y12[i-1],v9[i],canepsilon.xi9[i],dec4238.exceed,thresh[108]) -loglik(x12[i-1],y12[i-1],v9[i],w9[i-1],dec4238.exceed,thresh[108])) aprobepsilon.xi9[i]<-min(1,(clikely9*((dnorm(canepsilon.xi9[i],a.xi9,1/(sqrt(zeta.xi10)))/(dnorm(w9[i-1],a.xi9/(sqrt(zeta.xi10)))))))) if(u=aprobepsilon.xi9[i]){ w9[i]<-w9[i-1] } #eta.one.one[i]<-x1[i]+v1[i] #xi.one.one[i]<-y1[i]+w1[i] #***NINTH UPDATE DONE!***# a.eta10<-rnorm(1,(c*b+(zeta.eta10*(v1[i]+v2[i]+v3[i]+v4[i]+v5[i]+v6[i]+v7[i]+v8[i]+v9[i]+v10[i-1]+v11[i-1]+v12[i-1])))/(c+12*zeta.eta10),1/(sqrt(c+12*zeta.eta10))) a.xi10<-rnorm(1,(c*b+(zeta.xi10*(w1[i]+w2[i]+w3[i]+w4[i]+w5[i]+w6[i]+w7[i]+w8[i]+w9[i]+w10[i-1]+w11[i-1]+w12[i-1])))/(c+12*zeta.xi10),1/(sqrt(c+12*zeta.xi10))) zeta.eta11<-rgamma(1,f+6,1/(g+0.5*(((v1[i]-a.eta10)^2)+((v2[i]-a.eta10)^2)+((v3[i]-a.eta10)^2)+((v4[i]-a.eta10)^2)+((v5[i]-a.eta10)^2)+((v6[i]-a.eta10)^2)+((v7[i]-a.eta10)^2)+((v8[i]-a.eta10)^2)+((v9[i]-a.eta10)^2)+((v10[i-1]-a.eta10)^2)+((v11[i-1]-a.eta10)^2)+((v12[i-1]-a.eta10)^2)))) zeta.xi11<-rgamma(1,f+6,1/(g+0.5*(((w1[i]-a.xi10)^2)+((w2[i]-a.xi10)^2)+((w3[i]-a.xi10)^2)+((w4[i]-a.xi10)^2)+((w5[i]-a.xi10)^2)+((w6[i]-a.xi10)^2)+((w7[i]-a.xi10)^2)+((w8[i]-a.xi10)^2)+((w9[i]-a.xi10)^2)+((w10[i-1]-a.xi10)^2)+((w11[i-1]-a.xi10)^2)+((w12[i-1]-a.xi10)^2)))) tau.eta11<-rgamma(1,d+6,1/(e+0.5*(((x1[i])^2)+((x2[i])^2)+((x3[i])^2)+((x4[i])^2)+((x5[i])^2)+((x6[i])^2)+((x7[i])^2)+((x8[i])^2)+((x9[i])^2)+((x10[i-1])^2)+((x11[i-1])^2)+((x12[i-1])^2)))) tau.xi11<-rgamma(1,d+6,1/(e+0.5*(((y1[i])^2)+((y2[i])^2)+((y3[i])^2)+((y4[i])^2)+((y5[i])^2)+((y6[i])^2)+((y7[i])^2)+((y8[i])^2)+((y9[i])^2)+((y10[i-1])^2)+((y11[i-1])^2)+((y12[i-1])^2)))) cangamma.eta10[i]<-x10[i-1]+rnorm(1,0,errgamma.eta10) likely10<-exp(loglik(cangamma.eta10[i],y10[i-1],v1[i],w1[i],octshef.exceed,thresh[10]) -loglik(x10[i-1],y10[i-1],v1[i],w1[i],octshef.exceed,thresh[10]) +loglik(cangamma.eta10[i],y10[i-1],v2[i],w2[i],octbrad.exceed,thresh[22]) -loglik(x10[i-1],y10[i-1],v2[i],w2[i],octbrad.exceed,thresh[22]) +loglik(cangamma.eta10[i],y10[i-1],v3[i],w3[i],oct2245.exceed,thresh[34]) -loglik(x10[i-1],y10[i-1],v3[i],w3[i],oct2245.exceed,thresh[34]) +loglik(cangamma.eta10[i],y10[i-1],v4[i],w4[i],oct2255.exceed,thresh[46]) -loglik(x10[i-1],y10[i-1],v4[i-1],w4[i-1],oct2255.exceed,thresh[46]) +loglik(cangamma.eta10[i],y10[i-1],v5[i],w5[i],oct2416.exceed,thresh[58]) -loglik(x10[i-1],y10[i-1],v5[i],w5[i],oct2416.exceed,thresh[58]) +loglik(cangamma.eta10[i],y10[i-1],v6[i],w6[i],oct2432.exceed,thresh[70]) -loglik(x10[i-1],y10[i-1],v6[i],w6[i],oct2432.exceed,thresh[70]) +loglik(cangamma.eta10[i],y10[i-1],v7[i],w7[i],oct4036.exceed,thresh[82]) -loglik(x10[i-1],y10[i-1],v7[i],w7[i],oct4036.exceed,thresh[82]) +loglik(cangamma.eta10[i],y10[i-1],v8[i],w8[i],oct4206.exceed,thresh[94]) -loglik(x10[i-1],y10[i-1],v8[i],w8[i],oct4206.exceed,thresh[94]) +loglik(cangamma.eta10[i],y10[i-1],v9[i],w9[i],oct4238.exceed,thresh[106]) -loglik(x10[i-1],y10[i-1],v9[i],w9[i],oct4238.exceed,thresh[106]) +loglik(cangamma.eta10[i],y10[i-1],v10[i-1],w10[i-1],oct4617.exceed,thresh[118]) -loglik(x10[i-1],y10[i-1],v10[i-1],w10[i-1],oct4617.exceed,thresh[118]) +loglik(cangamma.eta10[i],y10[i-1],v11[i-1],w11[i-1],oct7221.exceed,thresh[130]) -loglik(x10[i-1],y10[i-1],v11[i-1],w11[i-1],oct7221.exceed,thresh[130]) +loglik(cangamma.eta10[i],y10[i-1],v12[i-1],w12[i-1],oct7377.exceed,thresh[142]) -loglik(x10[i-1],y10[i-1],v12[i-1],w12[i-1],oct7377.exceed,thresh[142])) aprobgamma.eta10[i]<-min(1,(likely10*((dnorm(cangamma.eta10[i],0,1/(sqrt(tau.eta11)))/(dnorm(x10[i-1],0,1/(sqrt(tau.eta11)))))))) u<-runif(1) if(u=aprobgamma.eta10[i]){ x10[i]<-x10[i-1] } cangamma.xi10[i]<-y10[i-1]+rnorm(1,0,errgamma.xi10) alikely10<-exp(loglik(x10[i],cangamma.xi10[i],v1[i],w1[i],octshef.exceed,thresh[10]) -loglik(x10[i],y10[i-1],v1[i],w1[i],octshef.exceed,thresh[10]) +loglik(x10[i],cangamma.xi10[i],v2[i],w2[i],octbrad.exceed,thresh[22]) -loglik(x10[i],y10[i-1],v2[i],w2[i],octbrad.exceed,thresh[22]) +loglik(x10[i],cangamma.xi10[i],v3[i],w3[i],oct2245.exceed,thresh[34]) -loglik(x10[i],y10[i-1],v3[i],w3[i],oct2245.exceed,thresh[34]) +loglik(x10[i],cangamma.xi10[i],v4[i],w4[i],oct2255.exceed,thresh[46]) -loglik(x10[i],y10[i-1],v4[i],w4[i],oct2255.exceed,thresh[46]) +loglik(x10[i],cangamma.xi10[i],v5[i],w5[i],oct2416.exceed,thresh[58]) -loglik(x10[i],y10[i-1],v5[i],w5[i],oct2416.exceed,thresh[58]) +loglik(x10[i],cangamma.xi10[i],v6[i],w6[i],oct2432.exceed,thresh[70]) -loglik(x10[i],y10[i-1],v6[i],w6[i],oct2432.exceed,thresh[70]) +loglik(x10[i],cangamma.xi10[i],v7[i],w7[i],oct4036.exceed,thresh[82]) -loglik(x10[i],y10[i-1],v7[i],w7[i],oct4036.exceed,thresh[82]) +loglik(x10[i],cangamma.xi10[i],v8[i],w8[i],oct4206.exceed,thresh[94]) -loglik(x10[i],y10[i-1],v8[i],w8[i],oct4206.exceed,thresh[94]) +loglik(x10[i],cangamma.xi10[i],v9[i],w9[i],oct4238.exceed,thresh[106]) -loglik(x10[i],y10[i-1],v9[i],w9[i],oct4238.exceed,thresh[106]) +loglik(x10[i],cangamma.xi10[i],v10[i-1],w10[i-1],oct4617.exceed,thresh[118]) -loglik(x10[i],y10[i-1],v10[i-1],w10[i-1],oct4617.exceed,thresh[118]) +loglik(x10[i],cangamma.xi10[i],v11[i-1],w11[i-1],oct7221.exceed,thresh[130]) -loglik(x10[i],y10[i-1],v11[i-1],w11[i-1],oct7221.exceed,thresh[130]) +loglik(x10[i],cangamma.xi10[i],v12[i-1],w12[i-1],oct7377.exceed,thresh[142]) -loglik(x10[i],y10[i-1],v12[i-1],w12[i-1],oct7377.exceed,thresh[142])) aprobgamma.xi10[i]<-min(1,(alikely10*((dnorm(cangamma.xi10[i],0,1/(sqrt(tau.xi11)))/(dnorm(y10[i-1],0,1/(sqrt(tau.xi11)))))))) if(u=aprobgamma.xi10[i]){ y10[i]<-y10[i-1] } canepsilon.eta10[i]<-v10[i-1]+rnorm(1,0,errepsilon.eta10) blikely10<-exp(loglik(x1[i],y1[i],canepsilon.eta10[i],w10[i-1],jan4617.exceed,thresh[109]) -loglik(x1[i],y1[i],v10[i-1],w10[i-1],jan4617.exceed,thresh[109]) +loglik(x2[i],y2[i],canepsilon.eta10[i],w10[i-1],feb4617.exceed,thresh[110]) -loglik(x2[i],y2[i],v10[i-1],w10[i-1],feb4617.exceed,thresh[110]) +loglik(x3[i],y3[i],canepsilon.eta10[i],w10[i-1],mar4617.exceed,thresh[111]) -loglik(x3[i],y3[i],v10[i-1],w10[i-1],mar4617.exceed,thresh[111]) +loglik(x4[i],y4[i],canepsilon.eta10[i],w10[i-1],apr4617.exceed,thresh[112]) -loglik(x4[i],y4[i],v10[i-1],w10[i-1],apr4617.exceed,thresh[112]) +loglik(x5[i],y5[i],canepsilon.eta10[i],w10[i-1],may4617.exceed,thresh[113]) -loglik(x5[i],y5[i],v10[i-1],w10[i-1],may4617.exceed,thresh[113]) +loglik(x6[i],y6[i],canepsilon.eta10[i],w10[i-1],jun4617.exceed,thresh[114]) -loglik(x6[i],y6[i],v10[i-1],w10[i-1],jun4617.exceed,thresh[114]) +loglik(x7[i],y7[i],canepsilon.eta10[i],w10[i-1],jul4617.exceed,thresh[115]) -loglik(x7[i],y7[i],v10[i-1],w10[i-1],jul4617.exceed,thresh[115]) +loglik(x8[i],y8[i],canepsilon.eta10[i],w10[i-1],aug4617.exceed,thresh[116]) -loglik(x8[i],y8[i],v10[i-1],w10[i-1],aug4617.exceed,thresh[116]) +loglik(x9[i],y9[i],canepsilon.eta10[i],w10[i-1],sep4617.exceed,thresh[117]) -loglik(x9[i],y9[i],v10[i-1],w10[i-1],sep4617.exceed,thresh[117]) +loglik(x10[i],y10[i],canepsilon.eta10[i],w10[i-1],oct4617.exceed,thresh[118]) -loglik(x10[i],y10[i],v10[i-1],w10[i-1],oct4617.exceed,thresh[118]) +loglik(x11[i-1],y11[i-1],canepsilon.eta10[i],w10[i-1],nov4617.exceed,thresh[119]) -loglik(x11[i-1],y11[i-1],v10[i-1],w10[i-1],nov4617.exceed,thresh[119]) +loglik(x12[i-1],y12[i-1],canepsilon.eta10[i],w10[i-1],dec4617.exceed,thresh[120]) -loglik(x12[i-1],y12[i-1],v10[i-1],w10[i-1],dec4617.exceed,thresh[120])) aprobepsilon.eta10[i]<-min(1,(blikely10*((dnorm(canepsilon.eta10[i],a.eta10,1/(sqrt(zeta.eta11)))/(dnorm(v10[i-1],a.eta10,1/(sqrt(zeta.eta11)))))))) if(u=aprobepsilon.eta10[i]){ v10[i]<-v10[i-1] } canepsilon.xi10[i]<-w10[i-1]+rnorm(1,0,errepsilon.xi10) clikely10<-exp(loglik(x1[i],y1[i],v10[i],canepsilon.xi10[i],jan4617.exceed,thresh[109]) -loglik(x1[i],y1[i],v10[i],w10[i-1],jan4617.exceed,thresh[109]) +loglik(x2[i],y2[i],v10[i],canepsilon.xi10[i],feb4617.exceed,thresh[110]) -loglik(x2[i],y2[i],v10[i],w10[i-1],feb4617.exceed,thresh[110]) +loglik(x3[i],y3[i],v10[i],canepsilon.xi10[i],mar4617.exceed,thresh[111]) -loglik(x3[i],y3[i],v10[i],w10[i-1],mar4617.exceed,thresh[111]) +loglik(x4[i],y4[i],v10[i],canepsilon.xi10[i],apr4617.exceed,thresh[112]) -loglik(x4[i],y4[i],v10[i],w10[i-1],apr4617.exceed,thresh[112]) +loglik(x5[i],y5[i],v10[i],canepsilon.xi10[i],may4617.exceed,thresh[113]) -loglik(x5[i],y5[i],v10[i],w10[i-1],may4617.exceed,thresh[113]) +loglik(x6[i],y6[i],v10[i],canepsilon.xi10[i],jun4617.exceed,thresh[114]) -loglik(x6[i],y6[i],v10[i],w10[i-1],jun4617.exceed,thresh[114]) +loglik(x7[i],y7[i],v10[i],canepsilon.xi10[i],jul4617.exceed,thresh[115]) -loglik(x7[i],y7[i],v10[i],w10[i-1],jul4617.exceed,thresh[115]) +loglik(x8[i],y8[i],v10[i],canepsilon.xi10[i],aug4617.exceed,thresh[116]) -loglik(x8[i],y8[i],v10[i],w10[i-1],aug4617.exceed,thresh[116]) +loglik(x9[i],y9[i],v10[i],canepsilon.xi10[i],sep4617.exceed,thresh[117]) -loglik(x9[i],y9[i],v10[i],w10[i-1],sep4617.exceed,thresh[117]) +loglik(x10[i],y10[i],v10[i],canepsilon.xi10[i],oct4617.exceed,thresh[118]) -loglik(x10[i],y10[i],v10[i],w10[i-1],oct4617.exceed,thresh[118]) +loglik(x11[i-1],y11[i-1],v10[i],canepsilon.xi10[i],nov4617.exceed,thresh[119]) -loglik(x11[i-1],y11[i-1],v10[i],w10[i-1],nov4617.exceed,thresh[119]) +loglik(x12[i-1],y12[i-1],v10[i],canepsilon.xi10[i],dec4617.exceed,thresh[120]) -loglik(x12[i-1],y12[i-1],v10[i],w10[i-1],dec4617.exceed,thresh[120])) aprobepsilon.xi10[i]<-min(1,(clikely10*((dnorm(canepsilon.xi10[i],a.xi10,1/(sqrt(zeta.xi11)))/(dnorm(w10[i-1],a.xi10/(sqrt(zeta.xi11)))))))) if(u=aprobepsilon.xi10[i]){ w10[i]<-w10[i-1] } #eta.one.one[i]<-x1[i]+v1[i] #xi.one.one[i]<-y1[i]+w1[i] #***TENTH UPDATE DONE!***# a.eta11<-rnorm(1,(c*b+(zeta.eta11*(v1[i]+v2[i]+v3[i]+v4[i]+v5[i]+v6[i]+v7[i]+v8[i]+v9[i]+v10[i]+v11[i-1]+v12[i-1])))/(c+12*zeta.eta11),1/(sqrt(c+12*zeta.eta11))) a.xi11<-rnorm(1,(c*b+(zeta.xi11*(w1[i]+w2[i]+w3[i]+w4[i]+w5[i]+w6[i]+w7[i]+w8[i]+w9[i]+w10[i]+w11[i-1]+w12[i-1])))/(c+12*zeta.xi11),1/(sqrt(c+12*zeta.xi11))) zeta.eta12<-rgamma(1,f+6,1/(g+0.5*(((v1[i]-a.eta11)^2)+((v2[i]-a.eta11)^2)+((v3[i]-a.eta11)^2)+((v4[i]-a.eta11)^2)+((v5[i]-a.eta11)^2)+((v6[i]-a.eta11)^2)+((v7[i]-a.eta11)^2)+((v8[i]-a.eta11)^2)+((v9[i]-a.eta11)^2)+((v10[i]-a.eta11)^2)+((v11[i-1]-a.eta11)^2)+((v12[i-1]-a.eta11)^2)))) zeta.xi12<-rgamma(1,f+6,1/(g+0.5*(((w1[i]-a.xi11)^2)+((w2[i]-a.xi11)^2)+((w3[i]-a.xi11)^2)+((w4[i]-a.xi11)^2)+((w5[i]-a.xi11)^2)+((w6[i]-a.xi11)^2)+((w7[i]-a.xi11)^2)+((w8[i]-a.xi11)^2)+((w9[i]-a.xi11)^2)+((w10[i]-a.xi11)^2)+((w11[i-1]-a.xi11)^2)+((w12[i-1]-a.xi11)^2)))) tau.eta12<-rgamma(1,d+6,1/(e+0.5*(((x1[i])^2)+((x2[i])^2)+((x3[i])^2)+((x4[i])^2)+((x5[i])^2)+((x6[i])^2)+((x7[i])^2)+((x8[i])^2)+((x9[i])^2)+((x10[i])^2)+((x11[i-1])^2)+((x12[i-1])^2)))) tau.xi12<-rgamma(1,d+6,1/(e+0.5*(((y1[i])^2)+((y2[i])^2)+((y3[i])^2)+((y4[i])^2)+((y5[i])^2)+((y6[i])^2)+((y7[i])^2)+((y8[i])^2)+((y9[i])^2)+((y10[i])^2)+((y11[i-1])^2)+((y12[i-1])^2)))) cangamma.eta11[i]<-x11[i-1]+rnorm(1,0,errgamma.eta11) likely11<-exp(loglik(cangamma.eta11[i],y11[i-1],v1[i],w1[i],novshef.exceed,thresh[11]) -loglik(x11[i-1],y11[i-1],v1[i],w1[i],novshef.exceed,thresh[11]) +loglik(cangamma.eta11[i],y11[i-1],v2[i],w2[i],novbrad.exceed,thresh[23]) -loglik(x11[i-1],y11[i-1],v2[i],w2[i],novbrad.exceed,thresh[23]) +loglik(cangamma.eta11[i],y11[i-1],v3[i],w3[i],nov2245.exceed,thresh[35]) -loglik(x11[i-1],y11[i-1],v3[i],w3[i],nov2245.exceed,thresh[35]) +loglik(cangamma.eta11[i],y11[i-1],v4[i],w4[i],nov2255.exceed,thresh[47]) -loglik(x11[i-1],y11[i-1],v4[i-1],w4[i-1],nov2255.exceed,thresh[47]) +loglik(cangamma.eta11[i],y11[i-1],v5[i],w5[i],nov2416.exceed,thresh[59]) -loglik(x11[i-1],y11[i-1],v5[i],w5[i],nov2416.exceed,thresh[59]) +loglik(cangamma.eta11[i],y11[i-1],v6[i],w6[i],nov2432.exceed,thresh[71]) -loglik(x11[i-1],y11[i-1],v6[i],w6[i],nov2432.exceed,thresh[71]) +loglik(cangamma.eta11[i],y11[i-1],v7[i],w7[i],nov4036.exceed,thresh[83]) -loglik(x11[i-1],y11[i-1],v7[i],w7[i],nov4036.exceed,thresh[83]) +loglik(cangamma.eta11[i],y11[i-1],v8[i],w8[i],nov4206.exceed,thresh[95]) -loglik(x11[i-1],y11[i-1],v8[i],w8[i],nov4206.exceed,thresh[95]) +loglik(cangamma.eta11[i],y11[i-1],v9[i],w9[i],nov4238.exceed,thresh[107]) -loglik(x11[i-1],y11[i-1],v9[i],w9[i],nov4238.exceed,thresh[107]) +loglik(cangamma.eta11[i],y11[i-1],v10[i],w10[i],nov4617.exceed,thresh[119]) -loglik(x11[i-1],y11[i-1],v10[i],w10[i],nov4617.exceed,thresh[119]) +loglik(cangamma.eta11[i],y11[i-1],v11[i-1],w11[i-1],nov7221.exceed,thresh[131]) -loglik(x11[i-1],y11[i-1],v11[i-1],w11[i-1],nov7221.exceed,thresh[131]) +loglik(cangamma.eta11[i],y11[i-1],v12[i-1],w12[i-1],nov7377.exceed,thresh[143]) -loglik(x11[i-1],y11[i-1],v12[i-1],w12[i-1],nov7377.exceed,thresh[143])) aprobgamma.eta11[i]<-min(1,(likely11*((dnorm(cangamma.eta11[i],0,1/(sqrt(tau.eta12)))/(dnorm(x11[i-1],0,1/(sqrt(tau.eta12)))))))) u<-runif(1) if(u=aprobgamma.eta11[i]){ x11[i]<-x11[i-1] } cangamma.xi11[i]<-y11[i-1]+rnorm(1,0,errgamma.xi11) alikely11<-exp(loglik(x11[i],cangamma.xi11[i],v1[i],w1[i],novshef.exceed,thresh[11]) -loglik(x11[i],y11[i-1],v1[i],w1[i],novshef.exceed,thresh[11]) +loglik(x11[i],cangamma.xi11[i],v2[i],w2[i],novbrad.exceed,thresh[23]) -loglik(x11[i],y11[i-1],v2[i],w2[i],novbrad.exceed,thresh[23]) +loglik(x11[i],cangamma.xi11[i],v3[i],w3[i],nov2245.exceed,thresh[35]) -loglik(x11[i],y11[i-1],v3[i],w3[i],nov2245.exceed,thresh[35]) +loglik(x11[i],cangamma.xi11[i],v4[i],w4[i],nov2255.exceed,thresh[47]) -loglik(x11[i],y11[i-1],v4[i],w4[i],nov2255.exceed,thresh[47]) +loglik(x11[i],cangamma.xi11[i],v5[i],w5[i],nov2416.exceed,thresh[59]) -loglik(x11[i],y11[i-1],v5[i],w5[i],nov2416.exceed,thresh[59]) +loglik(x11[i],cangamma.xi11[i],v6[i],w6[i],nov2432.exceed,thresh[71]) -loglik(x11[i],y11[i-1],v6[i],w6[i],nov2432.exceed,thresh[71]) +loglik(x11[i],cangamma.xi11[i],v7[i],w7[i],nov4036.exceed,thresh[83]) -loglik(x11[i],y11[i-1],v7[i],w7[i],nov4036.exceed,thresh[83]) +loglik(x11[i],cangamma.xi11[i],v8[i],w8[i],nov4206.exceed,thresh[95]) -loglik(x11[i],y11[i-1],v8[i],w8[i],nov4206.exceed,thresh[95]) +loglik(x11[i],cangamma.xi11[i],v9[i],w9[i],nov4238.exceed,thresh[107]) -loglik(x11[i],y11[i-1],v9[i],w9[i],nov4238.exceed,thresh[107]) +loglik(x11[i],cangamma.xi11[i],v10[i],w10[i],nov4617.exceed,thresh[119]) -loglik(x11[i],y11[i-1],v10[i],w10[i],nov4617.exceed,thresh[119]) +loglik(x11[i],cangamma.xi11[i],v11[i-1],w11[i-1],nov7221.exceed,thresh[131]) -loglik(x11[i],y11[i-1],v11[i-1],w11[i-1],nov7221.exceed,thresh[131]) +loglik(x11[i],cangamma.xi11[i],v12[i-1],w12[i-1],nov7377.exceed,thresh[143]) -loglik(x11[i],y11[i-1],v12[i-1],w12[i-1],nov7377.exceed,thresh[143])) aprobgamma.xi11[i]<-min(1,(alikely11*((dnorm(cangamma.xi11[i],0,1/(sqrt(tau.xi12)))/(dnorm(y11[i-1],0,1/(sqrt(tau.xi12)))))))) if(u=aprobgamma.xi11[i]){ y11[i]<-y11[i-1] } canepsilon.eta11[i]<-v11[i-1]+rnorm(1,0,errepsilon.eta11) blikely11<-exp(loglik(x1[i],y1[i],canepsilon.eta11[i],w11[i-1],jan7221.exceed,thresh[121]) -loglik(x1[i],y1[i],v11[i-1],w11[i-1],jan7221.exceed,thresh[121]) +loglik(x2[i],y2[i],canepsilon.eta11[i],w11[i-1],feb7221.exceed,thresh[122]) -loglik(x2[i],y2[i],v11[i-1],w11[i-1],feb7221.exceed,thresh[122]) +loglik(x3[i],y3[i],canepsilon.eta11[i],w11[i-1],mar7221.exceed,thresh[123]) -loglik(x3[i],y3[i],v11[i-1],w11[i-1],mar7221.exceed,thresh[123]) +loglik(x4[i],y4[i],canepsilon.eta11[i],w11[i-1],apr7221.exceed,thresh[124]) -loglik(x4[i],y4[i],v11[i-1],w11[i-1],apr7221.exceed,thresh[124]) +loglik(x5[i],y5[i],canepsilon.eta11[i],w11[i-1],may7221.exceed,thresh[125]) -loglik(x5[i],y5[i],v11[i-1],w11[i-1],may7221.exceed,thresh[125]) +loglik(x6[i],y6[i],canepsilon.eta11[i],w11[i-1],jun7221.exceed,thresh[126]) -loglik(x6[i],y6[i],v11[i-1],w11[i-1],jun7221.exceed,thresh[126]) +loglik(x7[i],y7[i],canepsilon.eta11[i],w11[i-1],jul7221.exceed,thresh[127]) -loglik(x7[i],y7[i],v11[i-1],w11[i-1],jul7221.exceed,thresh[127]) +loglik(x8[i],y8[i],canepsilon.eta11[i],w11[i-1],aug7221.exceed,thresh[128]) -loglik(x8[i],y8[i],v11[i-1],w11[i-1],aug7221.exceed,thresh[128]) +loglik(x9[i],y9[i],canepsilon.eta11[i],w11[i-1],sep7221.exceed,thresh[129]) -loglik(x9[i],y9[i],v11[i-1],w11[i-1],sep7221.exceed,thresh[129]) +loglik(x10[i],y10[i],canepsilon.eta11[i],w11[i-1],oct7221.exceed,thresh[130]) -loglik(x10[i],y10[i],v11[i-1],w11[i-1],oct7221.exceed,thresh[130]) +loglik(x11[i],y11[i],canepsilon.eta11[i],w11[i-1],nov7221.exceed,thresh[131]) -loglik(x11[i],y11[i],v11[i-1],w11[i-1],nov7221.exceed,thresh[131]) +loglik(x12[i-1],y12[i-1],canepsilon.eta11[i],w11[i-1],dec7221.exceed,thresh[132]) -loglik(x12[i-1],y12[i-1],v11[i-1],w11[i-1],dec7221.exceed,thresh[132])) aprobepsilon.eta11[i]<-min(1,(blikely11*((dnorm(canepsilon.eta11[i],a.eta11,1/(sqrt(zeta.eta12)))/(dnorm(v11[i-1],a.eta11,1/(sqrt(zeta.eta12)))))))) if(u=aprobepsilon.eta11[i]){ v11[i]<-v11[i-1] } canepsilon.xi11[i]<-w11[i-1]+rnorm(1,0,errepsilon.xi11) clikely11<-exp(loglik(x1[i],y1[i],v11[i],canepsilon.xi11[i],jan7221.exceed,thresh[121]) -loglik(x1[i],y1[i],v11[i],w11[i-1],jan7221.exceed,thresh[121]) +loglik(x2[i],y2[i],v11[i],canepsilon.xi11[i],feb7221.exceed,thresh[122]) -loglik(x2[i],y2[i],v11[i],w11[i-1],feb7221.exceed,thresh[122]) +loglik(x3[i],y3[i],v11[i],canepsilon.xi11[i],mar7221.exceed,thresh[123]) -loglik(x3[i],y3[i],v11[i],w11[i-1],mar7221.exceed,thresh[123]) +loglik(x4[i],y4[i],v11[i],canepsilon.xi11[i],apr7221.exceed,thresh[124]) -loglik(x4[i],y4[i],v11[i],w11[i-1],apr7221.exceed,thresh[124]) +loglik(x5[i],y5[i],v11[i],canepsilon.xi11[i],may7221.exceed,thresh[125]) -loglik(x5[i],y5[i],v11[i],w11[i-1],may7221.exceed,thresh[125]) +loglik(x6[i],y6[i],v11[i],canepsilon.xi11[i],jun7221.exceed,thresh[126]) -loglik(x6[i],y6[i],v11[i],w11[i-1],jun7221.exceed,thresh[126]) +loglik(x7[i],y7[i],v11[i],canepsilon.xi11[i],jul7221.exceed,thresh[127]) -loglik(x7[i],y7[i],v11[i],w11[i-1],jul7221.exceed,thresh[127]) +loglik(x8[i],y8[i],v11[i],canepsilon.xi11[i],aug7221.exceed,thresh[128]) -loglik(x8[i],y8[i],v11[i],w11[i-1],aug7221.exceed,thresh[128]) +loglik(x9[i],y9[i],v11[i],canepsilon.xi11[i],sep7221.exceed,thresh[129]) -loglik(x9[i],y9[i],v11[i],w11[i-1],sep7221.exceed,thresh[129]) +loglik(x10[i],y10[i],v11[i],canepsilon.xi11[i],oct7221.exceed,thresh[130]) -loglik(x10[i],y10[i],v11[i],w11[i-1],oct7221.exceed,thresh[130]) +loglik(x11[i],y11[i],v11[i],canepsilon.xi11[i],nov7221.exceed,thresh[131]) -loglik(x11[i],y11[i],v11[i],w11[i-1],nov7221.exceed,thresh[131]) +loglik(x12[i-1],y12[i-1],v11[i],canepsilon.xi11[i],dec7221.exceed,thresh[132]) -loglik(x12[i-1],y12[i-1],v11[i],w11[i-1],dec7221.exceed,thresh[132])) aprobepsilon.xi11[i]<-min(1,(clikely11*((dnorm(canepsilon.xi11[i],a.xi11,1/(sqrt(zeta.xi12)))/(dnorm(w11[i-1],a.xi11/(sqrt(zeta.xi12)))))))) if(u=aprobepsilon.xi11[i]){ w11[i]<-w11[i-1] } #eta.one.one[i]<-x1[i]+v1[i] #xi.one.one[i]<-y1[i]+w1[i] #***ELEVENTH UPDATE DONE!***# a.eta12<-rnorm(1,(c*b+(zeta.eta12*(v1[i]+v2[i]+v3[i]+v4[i]+v5[i]+v6[i]+v7[i]+v8[i]+v9[i]+v10[i]+v11[i]+v12[i-1])))/(c+12*zeta.eta12),1/(sqrt(c+12*zeta.eta12))) a.xi12<-rnorm(1,(c*b+(zeta.xi12*(w1[i]+w2[i]+w3[i]+w4[i]+w5[i]+w6[i]+w7[i]+w8[i]+w9[i]+w10[i]+w11[i]+w12[i-1])))/(c+12*zeta.xi12),1/(sqrt(c+12*zeta.xi12))) zeta.eta13<-rgamma(1,f+6,1/(g+0.5*(((v1[i]-a.eta12)^2)+((v2[i]-a.eta12)^2)+((v3[i]-a.eta12)^2)+((v4[i]-a.eta12)^2)+((v5[i]-a.eta12)^2)+((v6[i]-a.eta12)^2)+((v7[i]-a.eta12)^2)+((v8[i]-a.eta12)^2)+((v9[i]-a.eta12)^2)+((v10[i]-a.eta12)^2)+((v11[i]-a.eta12)^2)+((v12[i-1]-a.eta12)^2)))) zeta.xi13<-rgamma(1,f+6,1/(g+0.5*(((w1[i]-a.xi12)^2)+((w2[i]-a.xi12)^2)+((w3[i]-a.xi12)^2)+((w4[i]-a.xi12)^2)+((w5[i]-a.xi12)^2)+((w6[i]-a.xi12)^2)+((w7[i]-a.xi12)^2)+((w8[i]-a.xi12)^2)+((w9[i]-a.xi12)^2)+((w10[i]-a.xi12)^2)+((w11[i]-a.xi12)^2)+((w12[i-1]-a.xi12)^2)))) tau.eta13<-rgamma(1,d+6,1/(e+0.5*(((x1[i])^2)+((x2[i])^2)+((x3[i])^2)+((x4[i])^2)+((x5[i])^2)+((x6[i])^2)+((x7[i])^2)+((x8[i])^2)+((x9[i])^2)+((x10[i])^2)+((x11[i])^2)+((x12[i-1])^2)))) tau.xi13<-rgamma(1,d+6,1/(e+0.5*(((y1[i])^2)+((y2[i])^2)+((y3[i])^2)+((y4[i])^2)+((y5[i])^2)+((y6[i])^2)+((y7[i])^2)+((y8[i])^2)+((y9[i])^2)+((y10[i])^2)+((y11[i])^2)+((y12[i-1])^2)))) cangamma.eta12[i]<-x12[i-1]+rnorm(1,0,errgamma.eta12) likely12<-exp(loglik(cangamma.eta12[i],y12[i-1],v1[i],w1[i],decshef.exceed,thresh[12]) -loglik(x12[i-1],y12[i-1],v1[i],w1[i],decshef.exceed,thresh[12]) +loglik(cangamma.eta12[i],y12[i-1],v2[i],w2[i],decbrad.exceed,thresh[24]) -loglik(x12[i-1],y12[i-1],v2[i],w2[i],decbrad.exceed,thresh[24]) +loglik(cangamma.eta12[i],y12[i-1],v3[i],w3[i],dec2245.exceed,thresh[36]) -loglik(x12[i-1],y12[i-1],v3[i],w3[i],dec2245.exceed,thresh[36]) +loglik(cangamma.eta12[i],y12[i-1],v4[i],w4[i],dec2255.exceed,thresh[48]) -loglik(x12[i-1],y12[i-1],v4[i-1],w4[i-1],dec2255.exceed,thresh[48]) +loglik(cangamma.eta12[i],y12[i-1],v5[i],w5[i],dec2416.exceed,thresh[60]) -loglik(x12[i-1],y12[i-1],v5[i],w5[i],dec2416.exceed,thresh[60]) +loglik(cangamma.eta12[i],y12[i-1],v6[i],w6[i],dec2432.exceed,thresh[72]) -loglik(x12[i-1],y12[i-1],v6[i],w6[i],dec2432.exceed,thresh[72]) +loglik(cangamma.eta12[i],y12[i-1],v7[i],w7[i],dec4036.exceed,thresh[84]) -loglik(x12[i-1],y12[i-1],v7[i],w7[i],dec4036.exceed,thresh[84]) +loglik(cangamma.eta12[i],y12[i-1],v8[i],w8[i],dec4206.exceed,thresh[96]) -loglik(x12[i-1],y12[i-1],v8[i],w8[i],dec4206.exceed,thresh[96]) +loglik(cangamma.eta12[i],y12[i-1],v9[i],w9[i],dec4238.exceed,thresh[108]) -loglik(x12[i-1],y12[i-1],v9[i],w9[i],dec4238.exceed,thresh[108]) +loglik(cangamma.eta12[i],y12[i-1],v10[i],w10[i],dec4617.exceed,thresh[120]) -loglik(x12[i-1],y12[i-1],v10[i],w10[i],dec4617.exceed,thresh[120]) +loglik(cangamma.eta12[i],y12[i-1],v11[i],w11[i],dec7221.exceed,thresh[132]) -loglik(x12[i-1],y12[i-1],v11[i],w11[i],dec7221.exceed,thresh[132]) +loglik(cangamma.eta12[i],y12[i-1],v12[i-1],w12[i-1],dec7377.exceed,thresh[144]) -loglik(x12[i-1],y12[i-1],v12[i-1],w12[i-1],dec7377.exceed,thresh[144])) aprobgamma.eta12[i]<-min(1,(likely12*((dnorm(cangamma.eta12[i],0,1/(sqrt(tau.eta13)))/(dnorm(x12[i-1],0,1/(sqrt(tau.eta13)))))))) u<-runif(1) if(u=aprobgamma.eta12[i]){ x12[i]<-x12[i-1] } cangamma.xi12[i]<-y12[i-1]+rnorm(1,0,errgamma.xi12) alikely12<-exp(loglik(x12[i],cangamma.xi12[i],v1[i],w1[i],decshef.exceed,thresh[12]) -loglik(x12[i],y12[i-1],v1[i],w1[i],decshef.exceed,thresh[12]) +loglik(x12[i],cangamma.xi12[i],v2[i],w2[i],decbrad.exceed,thresh[24]) -loglik(x12[i],y12[i-1],v2[i],w2[i],decbrad.exceed,thresh[24]) +loglik(x12[i],cangamma.xi12[i],v3[i],w3[i],dec2245.exceed,thresh[36]) -loglik(x12[i],y12[i-1],v3[i],w3[i],dec2245.exceed,thresh[36]) +loglik(x12[i],cangamma.xi12[i],v4[i],w4[i],dec2255.exceed,thresh[48]) -loglik(x12[i],y12[i-1],v4[i],w4[i],dec2255.exceed,thresh[48]) +loglik(x12[i],cangamma.xi12[i],v5[i],w5[i],dec2416.exceed,thresh[60]) -loglik(x12[i],y12[i-1],v5[i],w5[i],dec2416.exceed,thresh[60]) +loglik(x12[i],cangamma.xi12[i],v6[i],w6[i],dec2432.exceed,thresh[72]) -loglik(x12[i],y12[i-1],v6[i],w6[i],dec2432.exceed,thresh[72]) +loglik(x12[i],cangamma.xi12[i],v7[i],w7[i],dec4036.exceed,thresh[84]) -loglik(x12[i],y12[i-1],v7[i],w7[i],dec4036.exceed,thresh[84]) +loglik(x12[i],cangamma.xi12[i],v8[i],w8[i],dec4206.exceed,thresh[96]) -loglik(x12[i],y12[i-1],v8[i],w8[i],dec4206.exceed,thresh[96]) +loglik(x12[i],cangamma.xi12[i],v9[i],w9[i],dec4238.exceed,thresh[108]) -loglik(x12[i],y12[i-1],v9[i],w9[i],dec4238.exceed,thresh[108]) +loglik(x12[i],cangamma.xi12[i],v10[i],w10[i],dec4617.exceed,thresh[120]) -loglik(x12[i],y12[i-1],v10[i],w10[i],dec4617.exceed,thresh[120]) +loglik(x12[i],cangamma.xi12[i],v11[i],w11[i],dec7221.exceed,thresh[132]) -loglik(x12[i],y12[i-1],v11[i],w11[i],dec7221.exceed,thresh[132]) +loglik(x12[i],cangamma.xi12[i],v12[i-1],w12[i-1],dec7377.exceed,thresh[144]) -loglik(x12[i],y12[i-1],v12[i-1],w12[i-1],dec7377.exceed,thresh[144])) aprobgamma.xi12[i]<-min(1,(alikely12*((dnorm(cangamma.xi12[i],0,1/(sqrt(tau.xi13)))/(dnorm(y12[i-1],0,1/(sqrt(tau.xi13)))))))) if(u=aprobgamma.xi12[i]){ y12[i]<-y12[i-1] } canepsilon.eta12[i]<-v12[i-1]+rnorm(1,0,errepsilon.eta12) blikely12<-exp(loglik(x1[i],y1[i],canepsilon.eta12[i],w12[i-1],jan7377.exceed,thresh[133]) -loglik(x1[i],y1[i],v12[i-1],w12[i-1],jan7377.exceed,thresh[133]) +loglik(x2[i],y2[i],canepsilon.eta12[i],w12[i-1],feb7377.exceed,thresh[134]) -loglik(x2[i],y2[i],v12[i-1],w12[i-1],feb7377.exceed,thresh[134]) +loglik(x3[i],y3[i],canepsilon.eta12[i],w12[i-1],mar7377.exceed,thresh[135]) -loglik(x3[i],y3[i],v12[i-1],w12[i-1],mar7377.exceed,thresh[135]) +loglik(x4[i],y4[i],canepsilon.eta12[i],w12[i-1],apr7377.exceed,thresh[136]) -loglik(x4[i],y4[i],v12[i-1],w12[i-1],apr7377.exceed,thresh[136]) +loglik(x5[i],y5[i],canepsilon.eta12[i],w12[i-1],may7377.exceed,thresh[137]) -loglik(x5[i],y5[i],v12[i-1],w12[i-1],may7377.exceed,thresh[137]) +loglik(x6[i],y6[i],canepsilon.eta12[i],w12[i-1],jun7377.exceed,thresh[138]) -loglik(x6[i],y6[i],v12[i-1],w12[i-1],jun7377.exceed,thresh[138]) +loglik(x7[i],y7[i],canepsilon.eta12[i],w12[i-1],jul7377.exceed,thresh[139]) -loglik(x7[i],y7[i],v12[i-1],w12[i-1],jul7377.exceed,thresh[139]) +loglik(x8[i],y8[i],canepsilon.eta12[i],w12[i-1],aug7377.exceed,thresh[140]) -loglik(x8[i],y8[i],v12[i-1],w12[i-1],aug7377.exceed,thresh[140]) +loglik(x9[i],y9[i],canepsilon.eta12[i],w12[i-1],sep7377.exceed,thresh[141]) -loglik(x9[i],y9[i],v12[i-1],w12[i-1],sep7377.exceed,thresh[141]) +loglik(x10[i],y10[i],canepsilon.eta12[i],w12[i-1],oct7377.exceed,thresh[142]) -loglik(x10[i],y10[i],v12[i-1],w12[i-1],oct7377.exceed,thresh[142]) +loglik(x11[i],y11[i],canepsilon.eta12[i],w12[i-1],nov7377.exceed,thresh[143]) -loglik(x11[i],y11[i],v12[i-1],w12[i-1],nov7377.exceed,thresh[143]) +loglik(x12[i],y12[i],canepsilon.eta12[i],w12[i-1],dec7377.exceed,thresh[144]) -loglik(x12[i],y12[i],v12[i-1],w12[i-1],dec7377.exceed,thresh[144])) aprobepsilon.eta12[i]<-min(1,(blikely12*((dnorm(canepsilon.eta12[i],a.eta12,1/(sqrt(zeta.eta13)))/(dnorm(v12[i-1],a.eta12,1/(sqrt(zeta.eta13)))))))) if(u=aprobepsilon.eta12[i]){ v12[i]<-v12[i-1] } canepsilon.xi12[i]<-w12[i-1]+rnorm(1,0,errepsilon.xi12) clikely12<-exp(loglik(x1[i],y1[i],v12[i],canepsilon.xi12[i],jan7377.exceed,thresh[133]) -loglik(x1[i],y1[i],v12[i],w12[i-1],jan7377.exceed,thresh[133]) +loglik(x2[i],y2[i],v12[i],canepsilon.xi12[i],feb7377.exceed,thresh[134]) -loglik(x2[i],y2[i],v12[i],w12[i-1],feb7377.exceed,thresh[134]) +loglik(x3[i],y3[i],v12[i],canepsilon.xi12[i],mar7377.exceed,thresh[135]) -loglik(x3[i],y3[i],v12[i],w12[i-1],mar7377.exceed,thresh[135]) +loglik(x4[i],y4[i],v12[i],canepsilon.xi12[i],apr7377.exceed,thresh[136]) -loglik(x4[i],y4[i],v12[i],w12[i-1],apr7377.exceed,thresh[136]) +loglik(x5[i],y5[i],v12[i],canepsilon.xi12[i],may7377.exceed,thresh[137]) -loglik(x5[i],y5[i],v12[i],w12[i-1],may7377.exceed,thresh[137]) +loglik(x6[i],y6[i],v12[i],canepsilon.xi12[i],jun7377.exceed,thresh[138]) -loglik(x6[i],y6[i],v12[i],w12[i-1],jun7377.exceed,thresh[138]) +loglik(x7[i],y7[i],v12[i],canepsilon.xi12[i],jul7377.exceed,thresh[139]) -loglik(x7[i],y7[i],v12[i],w12[i-1],jul7377.exceed,thresh[139]) +loglik(x8[i],y8[i],v12[i],canepsilon.xi12[i],aug7377.exceed,thresh[140]) -loglik(x8[i],y8[i],v12[i],w12[i-1],aug7377.exceed,thresh[140]) +loglik(x9[i],y9[i],v12[i],canepsilon.xi12[i],sep7377.exceed,thresh[141]) -loglik(x9[i],y9[i],v12[i],w12[i-1],sep7377.exceed,thresh[141]) +loglik(x10[i],y10[i],v12[i],canepsilon.xi12[i],oct7377.exceed,thresh[142]) -loglik(x10[i],y10[i],v12[i],w12[i-1],oct7377.exceed,thresh[142]) +loglik(x11[i],y11[i],v12[i],canepsilon.xi12[i],nov7377.exceed,thresh[143]) -loglik(x11[i],y11[i],v12[i],w12[i-1],nov7377.exceed,thresh[143]) +loglik(x12[i],y12[i],v12[i],canepsilon.xi12[i],dec7377.exceed,thresh[144]) -loglik(x12[i],y12[i],v12[i],w12[i-1],dec7377.exceed,thresh[144])) aprobepsilon.xi12[i]<-min(1,(clikely12*((dnorm(canepsilon.xi12[i],a.xi12,1/(sqrt(zeta.xi13)))/(dnorm(w12[i-1],a.xi12/(sqrt(zeta.xi13)))))))) if(u=aprobepsilon.xi12[i]){ w12[i]<-w12[i-1] } #eta.one.one[i]<-x1[i]+v1[i] #xi.one.one[i]<-y1[i]+w1[i] #***TWELVTH UPDATE DONE!***# aprobgamma.eta1<-aprobgamma.eta1[!is.na(aprobgamma.eta1)] aprobgamma.eta2<-aprobgamma.eta2[!is.na(aprobgamma.eta2)] aprobgamma.eta3<-aprobgamma.eta3[!is.na(aprobgamma.eta3)] aprobgamma.eta4<-aprobgamma.eta4[!is.na(aprobgamma.eta4)] aprobgamma.eta5<-aprobgamma.eta5[!is.na(aprobgamma.eta5)] aprobgamma.eta6<-aprobgamma.eta6[!is.na(aprobgamma.eta6)] aprobgamma.eta7<-aprobgamma.eta7[!is.na(aprobgamma.eta7)] aprobgamma.eta8<-aprobgamma.eta8[!is.na(aprobgamma.eta8)] aprobgamma.eta9<-aprobgamma.eta9[!is.na(aprobgamma.eta9)] aprobgamma.eta10<-aprobgamma.eta10[!is.na(aprobgamma.eta10)] aprobgamma.eta11<-aprobgamma.eta11[!is.na(aprobgamma.eta11)] aprobgamma.eta12<-aprobgamma.eta12[!is.na(aprobgamma.eta12)] aprobgamma.xi1<-aprobgamma.xi1[!is.na(aprobgamma.xi1)] aprobgamma.xi2<-aprobgamma.xi2[!is.na(aprobgamma.xi2)] aprobgamma.xi3<-aprobgamma.xi3[!is.na(aprobgamma.xi3)] aprobgamma.xi4<-aprobgamma.xi4[!is.na(aprobgamma.xi4)] aprobgamma.xi5<-aprobgamma.xi5[!is.na(aprobgamma.xi5)] aprobgamma.xi6<-aprobgamma.xi6[!is.na(aprobgamma.xi6)] aprobgamma.xi7<-aprobgamma.xi7[!is.na(aprobgamma.xi7)] aprobgamma.xi8<-aprobgamma.xi8[!is.na(aprobgamma.xi8)] aprobgamma.xi9<-aprobgamma.xi9[!is.na(aprobgamma.xi9)] aprobgamma.xi10<-aprobgamma.xi10[!is.na(aprobgamma.xi10)] aprobgamma.xi11<-aprobgamma.xi11[!is.na(aprobgamma.xi11)] aprobgamma.xi12<-aprobgamma.xi12[!is.na(aprobgamma.xi12)] aprobepsilon.eta1<-aprobepsilon.eta1[!is.na(aprobepsilon.eta1)] aprobepsilon.eta2<-aprobepsilon.eta2[!is.na(aprobepsilon.eta2)] aprobepsilon.eta3<-aprobepsilon.eta3[!is.na(aprobepsilon.eta3)] aprobepsilon.eta4<-aprobepsilon.eta4[!is.na(aprobepsilon.eta4)] aprobepsilon.eta5<-aprobepsilon.eta5[!is.na(aprobepsilon.eta5)] aprobepsilon.eta6<-aprobepsilon.eta6[!is.na(aprobepsilon.eta6)] aprobepsilon.eta7<-aprobepsilon.eta7[!is.na(aprobepsilon.eta7)] aprobepsilon.eta8<-aprobepsilon.eta8[!is.na(aprobepsilon.eta8)] aprobepsilon.eta9<-aprobepsilon.eta9[!is.na(aprobepsilon.eta9)] aprobepsilon.eta10<-aprobepsilon.eta10[!is.na(aprobepsilon.eta10)] aprobepsilon.eta11<-aprobepsilon.eta11[!is.na(aprobepsilon.eta11)] aprobepsilon.eta12<-aprobepsilon.eta12[!is.na(aprobepsilon.eta12)] aprobepsilon.xi1<-aprobepsilon.xi1[!is.na(aprobepsilon.xi1)] aprobepsilon.xi2<-aprobepsilon.xi2[!is.na(aprobepsilon.xi2)] aprobepsilon.xi3<-aprobepsilon.xi3[!is.na(aprobepsilon.xi3)] aprobepsilon.xi4<-aprobepsilon.xi4[!is.na(aprobepsilon.xi4)] aprobepsilon.xi5<-aprobepsilon.xi5[!is.na(aprobepsilon.xi5)] aprobepsilon.xi6<-aprobepsilon.xi6[!is.na(aprobepsilon.xi6)] aprobepsilon.xi7<-aprobepsilon.xi7[!is.na(aprobepsilon.xi7)] aprobepsilon.xi8<-aprobepsilon.xi8[!is.na(aprobepsilon.xi8)] aprobepsilon.xi9<-aprobepsilon.xi9[!is.na(aprobepsilon.xi9)] aprobepsilon.xi10<-aprobepsilon.xi10[!is.na(aprobepsilon.xi10)] aprobepsilon.xi11<-aprobepsilon.xi11[!is.na(aprobepsilon.xi11)] aprobepsilon.xi12<-aprobepsilon.xi12[!is.na(aprobepsilon.xi12)] #Section 2.6# return(aprobgamma.eta1,aprobgamma.eta2,aprobgamma.eta3,aprobgamma.eta4,aprobgamma.eta5,aprobgamma.eta6,aprobgamma.eta7,aprobgamma.eta8,aprobgamma.eta9,aprobgamma.eta10,aprobgamma.eta11,aprobgamma.eta12, aprobgamma.xi1,aprobgamma.xi2,aprobgamma.xi3,aprobgamma.xi4,aprobgamma.xi5,aprobgamma.xi6,aprobgamma.xi7,aprobgamma.xi8,aprobgamma.xi9,aprobgamma.xi10,aprobgamma.xi11,aprobgamma.xi12,aprobepsilon.eta1,aprobepsilon.eta2,aprobepsilon.eta3,aprobepsilon.eta4,aprobepsilon.eta5,aprobepsilon.eta6,aprobepsilon.eta7,aprobepsilon.eta8,aprobepsilon.eta9,aprobepsilon.eta10,aprobepsilon.eta11,aprobepsilon.eta12, aprobepsilon.xi1,aprobepsilon.xi2,aprobepsilon.xi3,aprobepsilon.xi4,aprobepsilon.xi5,aprobepsilon.xi6,aprobepsilon.xi7,aprobepsilon.xi8,aprobepsilon.xi9,aprobepsilon.xi10,aprobepsilon.xi11,aprobepsilon.xi12,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,y11,y12,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,w1,w2,w3,w4,w5,w6,w7,w8,w9,w10,w11,w12) }