# Split testing with RJMCMC

Hi all, this post is going to be short.

Last time I posted an evidence based split test system; this time I have augmented it to use the RJMCMC algorithm. The value of this method is that now the data is represented as a Zero-inflated Weibull, which is a more accurate way of describing ARPI.

The RJMCMC code samples all three Zero-inflated Weibull parameters; the percent conversion, ARPI and ARPPI can be evaluated form the posterior samples. The trans-d steps allow any of the three parameters to be different from each other for the two groups.

```The R code is:

g.LTGGD <-function(x = seq(0, 1, length.out = 10001), mu = 0.1, alpha= 0.01, beta = 2, P0 = 0.05  )
{
n = length(x)
LP = 1:n
LP[x == 0] = log(P0)
P1 = log(beta)
P2 = - log(2*alpha) - lgamma(1/beta)
P3 = -(abs(x-mu)/alpha)^beta
P4 = - log(1 - (0.5 - 0.5*pgamma(q= (mu/alpha)^beta,  shape =1/beta, rate=1 ) )   )
P5 = log(1-P0)
LP[x > 0] = P1 + P2 + P3[x > 0] + P4 + P5
return(LP)
}

g.ZDW <-function(x = seq(0, 1, length.out = 10001), kappa = 0.1,lambda= 0.01, P0 = 0.05  )
{
n = length(x)
LP = 1:n
LP[x == 0] = log(1-P0) # prob of zero
LP[x > 0 ] = dweibull(x[x > 0], shape = kappa, scale = lambda, log = T) + log(P0) # prob non_zro
return(LP)
}

g.data.maker<-function(N = 500, kappa1 = 2, lambda1= 9, P01 = 0.02, kappa2 = 2, lambda2= 9, P02 = 0.02)
{
Group = c(rep(1, N), rep(2, N))
x1    = rweibull(N, shape = kappa1, scale = lambda1)
x2    = rweibull(N, shape = kappa2, scale = lambda2)
r1 =  runif(N)
r2 =  runif(N)
x1[r1 >= P01] = 0
x2[r2 >= P02] = 0
x  = c(x1,x2)
x[x>=200] = 200
out = data.frame(Group, x)
return(out)
}
hist(g.data.maker(N=1000, kappa1 = 2, lambda1= 9, P01 = 1, kappa2 = 2, lambda2= 9, P02 = 1)[,2], breaks=1001, xlim=c(0,25))

g.BD.step <- function(m.old, LL.old, DP, d, sd.j, LB, UB)
{
order = sample(1:3, 3)
for ( j in order)
{
m.new      = m.old
m.new[j]   = abs(m.new[j] - 1)
m.new[j+6] = 0
#print(c(LB[j], UB[j]))
if(m.new[j] == 1) {m.new[j+6] = runif(1, min = LB[j+6], max = UB[j+6])}

LL.new.1 = sum(g.ZDW(x=d[d[,1]==1,2],  kappa = m.new[4]                   , lambda= m.new[5]                   , P0 = m.new[6]                     )  )
LL.new.2 = sum(g.ZDW(x=d[d[,1]==2,2],  kappa = m.new[4]+ m.new[1]*m.new[7], lambda= m.new[5]+ m.new[2]*m.new[8], P0 = m.new[6]+ m.new[3]*m.new[9]  )  )
LL.new   = LL.new.1 + LL.new.2

r  = runif(1)

MH.transd = 0.5 # birth
if (m.old[j] == 1 ) {MH.transd =2 } # death

MH = exp(LL.new - LL.old) * MH.transd
if (is.na(MH )) {MH = 0}

#print(MH)

if( (r <= MH) & ((DP !=2 ) | (m.old[j] ==0) )  ) # mh criteria + must not be perterm or birth, i.e., if perterb then must be birth
{
#print("accept")
m.old  = m.new
LL.old = LL.new
}
}

return(c(LL.old, m.old, sd.j))
}

g.P.step<- function(m.old, LL.old, DP, BURN=1, d, sd.j, LB, UB)
{
order = sample(4:9, (6))
for ( j in order)
{

#print("--------------------------------------------")
#print(j)

#if (i == 160) {return(REC[1:160,])}

m.new = m.old
step = 0

if( j <= 6)                         {step = rnorm(1, mean=0, sd = sd.j[j]   )}
if( (j ==  7) & (m.new[1] == 1 ) )  {step = rnorm(1, mean=0, sd = sd.j[j]   )}
if( (j ==  8) & (m.new[2] == 1 ) )  {step = rnorm(1, mean=0, sd = sd.j[j]   )}
if( (j ==  9) & (m.new[3] == 1 ) )  {step = rnorm(1, mean=0, sd = sd.j[j]   )}

m.new[j] = m.new[j] + step
#if(j == 9) {print(rbind(m.new, m.old))
#            print(sd.j[j])
#            print(step)
#            print(m.new[3])}
#print(rbind(m.old, m.new))
#print(sd.j)
#print(LB)
#print(UB)
#print(m.new)

if ( (m.new[j] < UB[j]) & (m.new[j] > LB[j]) &

(m.new[4] + m.new[1]*m.new[7] > LB[4]) & (m.new[5] + m.new[2]*m.new[8] > LB[5]) & (m.new[6] + m.new[3]*m.new[9] > LB[6]) &
(m.new[4] + m.new[1]*m.new[7] < UB[4]) & (m.new[5] + m.new[2]*m.new[8] < UB[5]) & (m.new[6] + m.new[3]*m.new[9] < UB[6]) ) # check prior bound
{

#print(rbind(c(m.old[3], m.old[4], m.old[7]), c(m.new[3], m.new[4], m.new[7])))

LL.new.1 = sum(g.ZDW(x=d[d[,1]==1,2],  kappa = m.new[4]                   , lambda= m.new[5]                   , P0 = m.new[6]                     )  )
LL.new.2 = sum(g.ZDW(x=d[d[,1]==2,2],  kappa = m.new[4]+ m.new[1]*m.new[7], lambda= m.new[5]+ m.new[2]*m.new[8], P0 = m.new[6]+ m.new[3]*m.new[9]  )  )
LL.new = LL.new.1 + LL.new.2
#if(j == 9) {print(c(LL.new, LL.old))}

r  = runif(1)
MH = exp(LL.new - LL.old)
if (is.na(MH )) {MH = 0}

#print(MH)

if( (r <= MH) & ((j <= 6) | (DP == 2)) )
{
#print("accept")
m.old  = m.new
LL.old = LL.new
}
}

}
#print(c(LL.old, m.old, sd.j))
#print("------------------>>>>>>> gavin this line<<<<<<<<<<<<<<<--------------------------")
return(c(LL.old, m.old, sd.j))
}

g.RJMCMC <-function(d, N=500, NB = 500)
{
ND = 9
#          1     2       3   4      5         6       7         8        9
#         Imu, Ialpha, IP0, kappa1, lambda1,  P01,    Dkappa2, Dlambda2, DP02
m.old  = c( 1,   1,       1, 2,      9,       0.02,    0,      0,         0. )

LL.old = sum(g.ZDW(x=d[,2],  kappa = m.old[4], lambda= m.old[5], P0 = m.old[6] ) )
print(LL.old)

REC = matrix(0, nrow= NB, ncol=ND+1)
REC[1,] = c(LL.old, m.old)
Names = c("LL", "I.kappa", "I.lambda", "I.P0", "Kappa", "Lambda", "P0", "Delta.Kappa", "Delta.Lamda", "Detla.P0" )
colnames(REC) = Names
sd.j = rep(0.1, ND)

UB = c(1,1,1, 50,  50, 1,    25,  25,     1 )
LB = c(0,0, 0, 0,   0,  0,  -25, -25,  -  1 )

# burnin
for ( i in 1:NB)
{

DP   = sample(1:2, 1)
#print(c(LL.old, m.old, sd.j))
#print("DB before")
junk = g.BD.step(m.old, LL.old, DP, d,sd.j, LB, UB )
#print(junk)
#print("DB")
LL.old = junk[1]
m.old  = junk[2:10]
sd.j   = junk[(11):(19)]

#print(c(LL.old, m.old, sd.j))
#print("P before")
junk = g.P.step(m.old, LL.old, DP, BURN =1,d, sd.j, LB, UB)
#print(junk)
#print("P")
LL.old = junk[1]
m.old  = junk[2:10]
sd.j   = junk[(11):(19)]
#print(i)

if ( i >= 150 )
{ for (j in 4:9)
{
sd.j[j] = sd(REC[(1):i, j+1])/sqrt(6)
sd.j[sd.j <= 0.0001] = 0.00001
}
}

REC[i,] = c(LL.old, m.old)
}

print(sd.j)

#first sample
for ( i in 1:NB)
{

DP   = sample(1:2, 1)
junk = g.BD.step(m.old, LL.old, DP, d,sd.j, LB, UB )
LL.old = junk[1]
m.old  = junk[2:10]
sd.j   = junk[(11):(19)]

junk = g.P.step(m.old, LL.old, DP, BURN = 1,d, sd.j, LB, UB )
LL.old = junk[1]
m.old  = junk[2:10]
sd.j   = junk[(11):(19)]

REC[i,] = c(LL.old, m.old)

if ( i >= 150 )
{ for (j in 4:9)
{
sd.j[j] = sd(REC[(1):i, j+1])/sqrt(6)
sd.j[sd.j <= 0.0001] = 0.00001
}
}

REC[i,] = c(LL.old, m.old)

}

alpha = 1
REC = data.frame(REC, index = 1:NB)

print(sd.j)
par(mfrow = c(3,3))
k = 0
for ( i in 1:N)
{
k = k + 1

DP   = sample(1:2, 1)
junk = g.BD.step(m.old, LL.old, DP, d,sd.j, LB, UB )
LL.old = junk[1]
m.old  = junk[2:10]
sd.j   = junk[(11):(19)]

junk = g.P.step(m.old, LL.old, DP, BURN =0,d, sd.j, LB, UB)
LL.old = junk[1]
m.old  = junk[2:10]
sd.j   = junk[(11):(19)]

alpha = alpha*NB/(alpha + NB)
#print(alpha)

r = runif(1)
if ( r <= alpha)
{
#print("accept")
row = sample(1:NB, 1 )
REC[row,] = c(LL.old, m.old, i+NB)

}

if ( k == 250)
{
k= 0
plot(REC[,11], REC[,2] , xlab ="Sample index", ylab= Names[2 ])
plot(REC[,11], REC[,3] , xlab ="Sample index", ylab= Names[3 ])
plot(REC[,11], REC[,4] , xlab ="Sample index", ylab= Names[4 ])
plot(REC[,11], REC[,5] , xlab ="Sample index", ylab= Names[5 ])
plot(REC[,11], REC[,6] , xlab ="Sample index", ylab= Names[6 ])
plot(REC[,11], REC[,7] , xlab ="Sample index", ylab= Names[7 ])
plot(REC[,11], REC[,8] , xlab ="Sample index", ylab= Names[8 ])
plot(REC[,11], REC[,9] , xlab ="Sample index", ylab= Names[9 ])
plot(REC[,11], REC[,10], xlab ="Sample index", ylab= Names[10 ])

MU = 1:ND
for ( i in 1:ND ) {MU[i] = mean(REC[,i+1])}

print(MU)
#Sigma = cor(REC[,2:10])
#print(Sigma)

#print(sd.j)

}

}

plot(REC[,11], REC[,2] , xlab ="Sample index", ylab= Names[2 ])
plot(REC[,11], REC[,3] , xlab ="Sample index", ylab= Names[3 ])
plot(REC[,11], REC[,4] , xlab ="Sample index", ylab= Names[4 ])
plot(REC[,11], REC[,5] , xlab ="Sample index", ylab= Names[5 ])
plot(REC[,11], REC[,6] , xlab ="Sample index", ylab= Names[6 ])
plot(REC[,11], REC[,7] , xlab ="Sample index", ylab= Names[7 ])
plot(REC[,11], REC[,8] , xlab ="Sample index", ylab= Names[8 ])
plot(REC[,11], REC[,9] , xlab ="Sample index", ylab= Names[9 ])
plot(REC[,11], REC[,10], xlab ="Sample index", ylab= Names[10 ])

}

data = g.data.maker(N=10000, kappa1 = 2, lambda1= 10, P01 = 0.01, kappa2 = 2, lambda2= 10, P02 = 0.015)

g.RJMCMC(d = data, N=50000, NB = 1000)```