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

In mable: Maximum Approximate Bernstein/Beta Likelihood Estimation

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
)

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

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.

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

Author(s)

Zhong Guan <[e-mail]>

Examples

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

mable documentation built on Oct. 1, 2024, 9:06 a.m.