# mable with degree selected by the method of moment and method of mode Maximum Approximate Bernstein/Beta Likelihood Estimation with an optimal model degree estimated by the Method of Moment ```r 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 ## Returns 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 ## 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 ```r ## 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) ``` ## Author(s) Zhong Guan

optimable function