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)
Advertisement