icons.R
In hyper2: The Hyperdirichlet Distribution, Mark 2

## -----------------------------------------------------------------------------
library("hyper2",quietly=TRUE)
M <- icons_table # saves typing
M

## -----------------------------------------------------------------------------
icons
icons == saffy(icons_table)  # should be TRUE

## -----------------------------------------------------------------------------
options("digits" = 4)
(mic <- maxp(icons))
dotchart(mic,pch=16)

## -----------------------------------------------------------------------------
L1 <- loglik(indep(mic),icons)
L1

## -----------------------------------------------------------------------------
equalp.test(icons)

## -----------------------------------------------------------------------------
specificp.test(icons,1)

## -----------------------------------------------------------------------------
o <- function(Ul,Cl,startp,give=FALSE){
    small <- 1e-4  #  ensure start at an interior point
    if(missing(startp)){startp <- small*(1:5)+rep(0.1,5)}			
    out <- maxp(icons, startp=small*(1:5)+rep(0.1,5), give=TRUE, fcm=Ul,fcv=Cl)
    if(give){
        return(out)
    }else{
        return(out$value)
    }
}

p2max <- o(c(-1, 1, 0, 0, 0), 0)  # p1 <= p2
p3max <- o(c(-1, 0, 1, 0, 0), 0)  # p1 <= p3
p4max <- o(c(-1, 0, 0, 1, 0), 0)  # p1 <= p4
p5max <- o(c(-1, 0, 0, 0, 1), 0)  # p1 <= p5
p6max <- o(c(-2,-1,-1,-1,-1),-1)  # p1 <= p6 (fillup)

## -----------------------------------------------------------------------------
likes <- c(p2max,p3max,p4max,p5max,p6max)
likes
ml <- max(likes) 
ml

## -----------------------------------------------------------------------------
L1-ml

## -----------------------------------------------------------------------------
o2 <- function(Ul,Cl){
  jj <-o(Ul,Cl,give=TRUE)
  out <- c(jj[[1]],1-sum(jj[[1]]),jj[[2]])
  names(out) <- c("p1","p2","p3","p4","p5","p6","support")
  return(out)
}
rbind(
o2(c(-1, 1, 0, 0, 0), 0),  # p1 <= p2
o2(c(-1, 0, 1, 0, 0), 0),  # p1 <= p3
o2(c(-1, 0, 0, 1, 0), 0),  # p1 <= p4
o2(c(-1, 0, 0, 0, 1), 0),  # p1 <= p5
o2(c(-2,-1,-1,-1,-1),-1)   # p1 <= p6
)

## -----------------------------------------------------------------------------
jj <- o(c(-1,-1,-1,-1,0) , -2/3, give=TRUE,start=indep((1:6)/21))$value
jj

## -----------------------------------------------------------------------------
L1-jj

## -----------------------------------------------------------------------------
small <- 1e-4
start <- indep(c(small,small,small,small,0.5-2*small,0.5-2*small))
jj <- c(
   o(c(-1, 0, 0, 0, 1), 0,start=start),
   o(c( 0,-1, 0, 0, 1), 0,start=start),
   o(c( 0, 0,-1, 0, 1), 0,start=start),
   o(c( 0, 0, 0,-1, 1), 0,start=start),

   o(c(-2,-1,-1,-1,-1),-1,start=start),
   o(c(-1,-2,-1,-1,-1),-1,start=start),
   o(c(-1,-1,-2,-1,-1),-1,start=start),
   o(c(-1,-1,-1,-2,-1),-1,start=start)
   )
jj
max(jj)

## -----------------------------------------------------------------------------
L1-max(jj)

## -----------------------------------------------------------------------------
   o(c( 0, 0, 0,-1, 1), 0,give=TRUE,start=start)