mable (version 4.1.1)

optimable: mable with degree selected by the method of moment and method of mode

Description

Maximum Approximate Bernstein/Beta Likelihood Estimation with an optimal model degree estimated by the Method of Moment

Usage

optimable(
  x,
  interval,
  m = NULL,
  mu = NULL,
  lam = NULL,
  modes = NULL,
  nmod = 1,
  ushaped = FALSE,
  maxit = 50L
)

Value

A class "mable" object with components

• m the given or a selected degree by method of change-point

• p the estimated vector of mixture proportions \(p = (p_0, \ldots, p_m)\) with the selected/given optimal degree m

• mloglik the maximum log-likelihood at degree m

• interval support/truncation interval (a,b)

• convergence An integer code. 0 indicates successful completion (all the EM iterations are convergent and an optimal degree is successfully selected in M). Possible error codes are

  • 1, indicates that the iteration limit maxit had been reached in at least one EM iteration;

  • 2, the search did not finish before m1.

• delta the convergence criterion delta value

Arguments

x

a univariate sample data in interval

interval

a closed interval c(a,b), default is [0,1]

m

initial degree, default is 2 times the number of modes nmod.

mu

a vector of component means of multimodal mixture density, default is NULL for unimodal or unknown

lam

a vector of mixture proportions of same length of mu

modes

a vector of the locations of modes, if it is NULL (default) and multimode::locmodes()

nmod

the number of modes, if nmod=0, the lower bound for m is estimated based on mean and variance only.

ushaped

logical, whether or not the density is clearly U-shaped including J- and L-shaped with mode occurs at the endpoint of the support.

maxit

maximum iterations

Author

Zhong Guan <[e-mail]>

Details

If the data show a clear uni- or multi-modal distribution, then give the value of nmod as the number of modes. Otherwise nmod=0. The degree is estimated by the iterative method of moment with an initial degree estimated by the method of mode. For multimodal density, if useful estimates of the component means mu and proportions lam are available then they can be used to give an initial degree. If the distribution is clearly U-, J-, or L-shaped, i.e., the mode occurs at the endpoint of interval, then set ushaped=TRUE. In this case the degree is estimated by the method of mode.

Examples

# \donttest{
## Old Faithful Data
x<-faithful
x1<-faithful[,1]
x2<-faithful[,2]
a<-c(0, 40); b<-c(7, 110)
mu<-(apply(x,2,mean)-a)/(b-a)
s2<-apply(x,2,var)/(b-a)^2
# mixing proportions
lambda<-c(mean(x1<3),mean(x2<65))
# guess component mean
mu1<-(c(mean(x1[x1<3]), mean(x2[x2<65]))-a)/(b-a)  
mu2<-(c(mean(x1[x1>=3]), mean(x2[x2>=65]))-a)/(b-a)  
# estimate lower bound for m
mb<-ceiling((mu*(1-mu)-s2)/(s2-lambda*(1-lambda)*(mu1-mu2)^2)-2)
mb
m1<-optimable(x1, interval=c(a[1],b[1]), nmod=2, modes=c(2,4.5))$m 
m2<-optimable(x2, interval=c(a[2],b[2]), nmod=2, modes=c(52.5,80))$m 
m1;m2
erupt1<-mable(x1, M=mb[1], interval=c(a[1],b[1]))
erupt2<-mable(x1, M=m1, interval=c(a[1],b[1])) 
wait1<-mable(x2, M=mb[2],interval=c(a[2],b[2])) 
wait2<-mable(x2, M=m2,interval=c(a[2],b[2])) 
ans1<- mable.mvar(faithful, M = mb, search =FALSE, method="em",  interval = cbind(a,b))
ans2<- mable.mvar(faithful, M = c(m1,m2), search =FALSE, method="em", interval = cbind(a,b))
op<-par(mfrow=c(1,2), cex=0.8)
hist(x1, probability = TRUE, col="grey", border="white", main="", 
      xlab="Eruptions", ylim=c(0,.65), las=1)
plot(erupt1, add=TRUE,"density")
plot(erupt2, add=TRUE,"density",lty=2,col=2)
legend("topleft", lty=c(1,2),col=1:2, bty="n", cex=.7,  
      c(expression(paste("m = ", m[b])),expression(paste("m = ", hat(m)))))
hist(x2, probability = TRUE, col="grey", border="white", main="", 
      xlab="Waiting", las=1)
plot(wait1, add=TRUE,"density")
plot(wait2, add=TRUE,"density",lty=2,col=2)
legend("topleft", lty=c(1,2),col=1:2, bty="n", cex=.7,  
      c(expression(paste("m = ", m[b])),expression(paste("m = ", hat(m)))))
par(op)
op<-par(mfrow=c(1,2), cex=0.7)
plot(ans1, which="density", contour=TRUE) 
plot(ans2, which="density", contour=TRUE, add=TRUE, lty=2, col=2) 
plot(ans1, which="cumulative", contour=TRUE)       
plot(ans2, which="cumulative", contour=TRUE, add=TRUE, lty=2, col=2)
par(op)
# }