require("SQUAREM")
#####---------------------------------------------------------------###########
##For MM algorithm, we will illustrate the acceleration of EM by Squarem ######
##using logistics regression maximum likelihood estimation example. ###########
#####---------------------------------------------------------------###########

#####---------------------------------------------------------------###########
###########################data################################################
#####---------------------------------------------------------------###########
ld <- c(1, 0.80, 0.83, 0.66, 1.9, 1.100, 0.996,	1, 
        1, 0.90, 0.36, 0.32, 1.4, 0.740, 0.992,	1,
        1, 0.80, 0.88, 0.70, 0.8, 0.176, 0.982,	0,
        1, 1.00, 0.87, 0.87, 0.7, 1.053, 0.986,	0,
        1, 0.90, 0.75, 0.68, 1.3, 0.519, 0.980,	1,
        1, 1.00, 0.65, 0.65, 0.6, 0.519, 0.982,	0,
        1, 0.95, 0.97, 0.92, 1.0, 1.230, 0.992,	1,
        1, 0.95, 0.87, 0.83, 1.9, 1.354, 1.020,	0,
        1, 1.00, 0.45, 0.45, 0.8, 0.322, 0.999,	0,
        1, 0.95, 0.36, 0.34, 0.5, 0.000, 1.038,	0,
        1, 0.85, 0.39, 0.33, 0.7, 0.279, 0.988,	0,
        1, 0.70, 0.76, 0.53, 1.2, 0.146, 0.982,	0,
        1, 0.80, 0.46, 0.37, 0.4, 0.380, 1.006,	0,
        1, 0.20, 0.39, 0.08, 0.8, 0.114, 0.990,	0,
        1, 1.00, 0.90, 0.90, 1.1, 1.037, 0.990,	0,
        1, 1.00, 0.84, 0.84, 1.9, 2.064, 1.020,	1,
        1, 0.65, 0.42, 0.27, 0.5, 0.114, 1.014,	0,
        1, 1.00, 0.75, 0.75, 1.0, 1.322, 1.004,	0,
        1, 0.50, 0.44, 0.22, 0.6, 0.114, 0.990,	0,
        1, 1.00, 0.63, 0.63, 1.1, 1.072, 0.986,	1,
        1, 1.00, 0.33, 0.33, 0.4, 0.176, 1.010,	0,
        1, 0.90, 0.93, 0.84, 0.6, 1.591, 1.020,	0,
        1, 1.00, 0.58, 0.58, 1.0, 0.531, 1.002,	1,
        1, 0.95, 0.32, 0.30, 1.6, 0.886, 0.988,	0,
        1, 1.00, 0.60, 0.60, 1.7, 0.964, 0.990,	1,
        1, 1.00, 0.69, 0.69, 0.9, 0.398, 0.986,	1,
        1, 1.00, 0.73, 0.73, 0.7, 0.398, 0.986,	0)
ld <- matrix(ld, byrow = T, ncol = 8)
ld <- data.frame(ld)

#####---------------------------------------------------------------###########
######################starting value###########################################
#####---------------------------------------------------------------###########
Z <- as.matrix(ld[, 1:7])
y <- ld[, 8]
p0 <- rep(10, 7)


#####---------------------------------------------------------------###########
####The fixed point mapping giving MM algorithm--uniform bound#################
#####---------------------------------------------------------------###########
qmub.update <- function(par, Z, y){
        Zmat <- solve(crossprod(Z)) %*% t(Z)
        zb <- c(Z %*% par)
        pib <- 1 / (1 + exp(-zb))
        ub <-  pib - y
        par <- par - 4 * c(Zmat %*% ub)
        par
}

#####---------------------------------------------------------------###########
####The fixed point mapping giving MM algorithm--non uniform bound#############
#####---------------------------------------------------------------###########
qmvb.update <- function(par, Z, y){
        zb <- c(Z %*% par)
        pib <- 1 / (1 + exp(-zb))
        wmat <- diag((2 * pib - 1)/(2 * zb))
        ub <-  pib - y
        Zmat <- solve(t(Z) %*% wmat %*% Z) %*% t(Z)
        par <- par - c(Zmat %*% ub)
        par
}

#####---------------------------------------------------------------###########
####Objective function whose local minimum is a fixed point. ##################
####Here it is the negative log-likelihood of logistic regression.#############
#####---------------------------------------------------------------###########
binom.loglike <- function(par, Z, y){
        zb <- c(Z %*% par)
        pib <- 1 / (1 + exp(-zb))
        return(as.numeric(-t(y) %*% (Z %*% par) - sum(log(1 - pib))))
}


#####---------------------------------------------------------------###########
#############################uniform-bound MM Algorithm########################
#####---------------------------------------------------------------###########
system.time(ans1 <- fpiter(par = p0, fixptfn = qmub.update, 
                           objfn = binom.loglike, 
                           control = list(maxiter = 20000), 
                           Z = Z, y = y))
ans1

#####---------------------------------------------------------------###########
############Squarem to accelerate uniform-bound MM Algorithm###################
#####---------------------------------------------------------------###########
system.time(ans2 <- squarem(par = p0, fixptfn = qmub.update, 
                            objfn = binom.loglike, 
                            Z = Z, y = y))
ans2

#####---------------------------------------------------------------###########
#########################non-uniform-bound MM Algorithm########################
#####---------------------------------------------------------------###########
system.time(ans3 <- fpiter(par = p0, fixptfn = qmvb.update, 
                           objfn = binom.loglike, 
                           control = list(maxiter = 20000),
                           Z = Z, y = y))
ans3


#####---------------------------------------------------------------###########
############Squarem to accelerate non uniform-bound MM Algorithm###############
#####---------------------------------------------------------------###########
system.time(ans4 <- squarem(par = p0, fixptfn = qmvb.update, 
                            objfn = binom.loglike, 
                            Z = Z, y = y))
ans4