-
Notifications
You must be signed in to change notification settings - Fork 1
/
Example_SimA.Rmd
93 lines (77 loc) · 2.31 KB
/
Example_SimA.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
---
title: "Example Script Simulation (a)"
output:
pdf_document: default
html_document:
df_print: paged
---
#### Load functions
```{r message=FALSE, warning=FALSE, include=FALSE, paged.print=FALSE}
gc(verbose=FALSE)
source("MAIN_FUNCTION.R",verbose=FALSE)
```
#### Data Generation from a mixture model (See Section 3 and Table 1(a))
```{r}
c=40
g_full <- graph.lattice(length=c,dim=2)
n=c^2
#V(g_full)$index=1:n
net=get.adjacency(g_full,attr=NULL)
K=5
rho=5
label_T0<-pottshm(ncol=K,niter=1,c,m=c,beta=rho)
label_T1<-pottshm(ncol=K,niter=1,c,m=c,beta=rho)
cc=c/4
#####Partion the image
label_T0[1:cc,]<-1
label_T0[(cc+1):(cc*2),]<-2
label_T0[(cc*2+1):(cc*3),]<-3
label_T0[(cc*3+1):c,]<-4
Matrix_list=lapply(1:K, function(k) rwish(30, (k+1)*diag(3)/30))
Matrix_list[[K]]<-rwish(30, 1.5*diag(3)/50)
label_T1[1:cc,]<-1
label_T1[(cc+1):(cc*2),]<-2
label_T1[(cc*2+1):(cc*3),]<-3
label_T1[(cc*3+1):c,]<-4
#####Indicate the reg
label_T1[(cc+1):(cc*2),(cc+1):(cc*2)]<-5
Data=list()
N_T0=5
N_T1=5
T_index=c(rep(1,N_T0),rep(2,N_T1))
Data_T0<-lapply(1:N_T0, function(i) lapply(1:c^2,function(s) riwish(5,Matrix_list[[label_T0[s]]])))
Data_T1<-lapply(1:N_T1, function(i) lapply(1:c^2,function(s) riwish(5,Matrix_list[[label_T1[s]]])))
Data<-append(Data_T0,Data_T1)
par(mfrow=c(1,2))
image(label_T0)
image(label_T1)
p=3
n=c^2
K=10
N=N_T0+N_T1
no.T=2
Ca=3
Cb=3
adj_list<- apply(net==1,1,which)
S=diag(2)*0.1
rho_alpha=0.3
rho_beta=2
```
#### Codes implementation
```{r message=FALSE, warning=FALSE, paged.print=FALSE}
result<-Potts_Bayesian_Semi_Multi(Data,net,n,p=3,N,
K=10,T_index,iter=10000,
LU_alpha=c(0,20),
LU_beta=c(0,20),
LU_xi=c(0,1),
rho_alpha=2,
rho_beta=4,
rho_xi=0.1)
```
#### Posterior Mean of $mathcal{I}(h_{0v}\not=h_{1v})$
```{r message=FALSE, warning=FALSE, include=FALSE, paged.print=FALSE}
dd=rbind(rep(0,40),rep(0,40),(1*(matrix(colMeans(result$Diff[101:500,])>0.5,40,40)))[3:40,])
```
```{r}
image(matrix(colMeans(result$Diff[101:500,]),40,40))
```