### Install and get started with mixsqp in R Source: https://stephenslab.github.io/mixsqp/index.html Steps to install and explore the mixsqp package. First install from CRAN, then load the library contact, run a built-in example, and open the help page. ```R install.packages("mixsqp") ``` ```R library(mixsqp) ``` ```R example("mixsqp") ``` ```R help("mixsqp") ``` -------------------------------- ### Build and test mixsqp package with R CMD commands Source: https://stephenslab.github.io/mixsqp/index.html Build, install, and check the mixsqp package from the command line in a shell. Requires dependencies to be installed. ```shell R CMD build mixsqp R CMD INSTALL mixsqp_0.3-18.tar.gz R CMD check --as-cran mixsqp_0.3-18.tar.gz ``` -------------------------------- ### Install mixsqp from GitHub with devtools Source: https://stephenslab.github.io/mixsqp/index.html Install the latest mixsqp version from GitHub using devtools. This command automatically installs required packages. ```R install.packages("devtools") library(devtools) install_github("stephenslab/mixsqp",build_vignettes = TRUE) ``` -------------------------------- ### Install local mixsqp repository Source: https://stephenslab.github.io/mixsqp/index.html Install mixsqp from a local clone of the repository using the install_local function from devtools ```R library(devtools) install_local("mixsqp") ``` -------------------------------- ### Build and view mixsqp vignette Source: https://stephenslab.github.io/mixsqp/index.html Build the vignette and view it locally after installing REBayes and Rmosek for the optional 'mixkwdual' solver. ```R devtools::install_github("stephenslab/mixsqp",build_vignettes = TRUE) library(mixsqp) vignette("mixsqp-intro") ``` -------------------------------- ### Fit IP algorithm with mixkwdual Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html Computes the solution using the IP algorithm from the REBayes package. Requires REBayes to be installed; otherwise this step can be skipped. ```R fit.ip <- mixkwdual(L) ``` -------------------------------- ### tacks Source: https://stephenslab.github.io/mixsqp/reference/index.html Beckett & Diaconis tack rolling example. ```APIDOC ## tacks ### Description Beckett & Diaconis tack rolling example. ``` -------------------------------- ### Generate likelihood matrix with simulatemixdata Source: https://stephenslab.github.io/mixsqp/reference/simulatemixdata.html Example usage of `simulatemixdata` to generate a likelihood matrix for 1,000 samples and 20 mixture components. The result is stored in `dat`. ```R # Generate the likelihood matrix for a data set with 1,000 samples # and a nonparametric Empirical Bayes model with 20 mixture # components. dat <- simulatemixdata(1000,20) ``` -------------------------------- ### Fit mixture proportions with mixsqp() and evaluate with mixobjective() Source: https://stephenslab.github.io/mixsqp/reference/mixsqp.html Basic call to mixsqp on a simulated 100000 x 10 matrix, then evaluates the objective with mixobjective. The run prints algorithm settings and a per-iteration table, ending when convergence criteria are met in about 0.5 seconds. The optional comparison with mixkwdual requires MOSEK and REBayes, so it is wrapped in if (FALSE). ```R __ set.seed(1) n <- 1e5 m <- 10 w <- rep(1,n)/n L <- simulatemixdata(n,m)$L out.mixsqp <- mixsqp(L,w) #> Running mix-SQP algorithm 0.3-17 on 100000 x 10 matrix #> convergence tol. (SQP): 1.0e-08 #> conv. tol. (active-set): 1.0e-10 #> zero threshold (solution): 1.0e-08 #> zero thresh. (search dir.): 1.0e-10 #> l.s. sufficient decrease: 1.0e-02 #> step size reduction factor: 7.5e-01 #> minimum step size: 1.0e-08 #> max. iter (SQP): 1000 #> max. iter (active-set): 11 #> number of EM iterations: 10 #> Computing SVD of 100000 x 10 matrix. #> Matrix is not low-rank; falling back to full matrix. #> iter objective max(rdual) nnz stepsize max.diff nqp nls #> 1 +6.649406794e-01 -- EM -- 10 1.00e+00 7.79e-02 -- -- #> 2 +6.336404241e-01 -- EM -- 10 1.00e+00 2.61e-02 -- -- #> 3 +6.201431305e-01 -- EM -- 10 1.00e+00 2.41e-02 -- -- #> 4 +6.131770842e-01 -- EM -- 10 1.00e+00 2.06e-02 -- -- #> 5 +6.091788165e-01 -- EM -- 10 1.00e+00 1.71e-02 -- -- #> 6 +6.066948308e-01 -- EM -- 10 1.00e+00 1.42e-02 -- -- #> 7 +6.050509636e-01 -- EM -- 10 1.00e+00 1.18e-02 -- -- #> 8 +6.039062025e-01 -- EM -- 10 1.00e+00 9.95e-03 -- -- #> 9 +6.030755992e-01 -- EM -- 10 1.00e+00 8.45e-03 -- -- #> 10 +6.024525768e-01 -- EM -- 10 1.00e+00 7.25e-03 -- -- #> 1 +6.024525768e-01 +2.365e-02 10 ------ ------ -- -- #> 2 +6.058450302e-01 +4.073e+00 4 1.00e+00 3.91e-01 11 1 #> 3 +5.999690521e-01 +1.894e+00 5 1.00e+00 3.21e-02 10 1 #> 4 +5.997818561e-01 +7.929e-01 5 1.00e+00 2.95e-04 2 1 #> 5 +5.997500584e-01 +2.749e-01 5 1.00e+00 1.98e-03 2 1 #> 6 +5.996811464e-01 +1.482e-01 5 1.00e+00 5.64e-03 2 1 #> 7 +5.994718234e-01 +6.121e-02 5 1.00e+00 9.54e-03 2 1 #> 8 +5.991459366e-01 +1.505e-02 5 1.00e+00 8.23e-03 2 1 #> 9 +5.989768704e-01 +1.193e-03 5 1.00e+00 2.88e-03 2 1 #> 10 +5.989591648e-01 +8.215e-06 5 1.00e+00 2.42e-04 2 1 #> 11 +5.989590364e-01 +3.904e-10 5 1.00e+00 1.63e-06 2 1 #> Optimization took 0.54 seconds. #> Convergence criteria met---optimal solution found. f <- mixobjective(L,out.mixsqp$x,w) print(f,digits = 16) #> [1] 0.5989590363728419 # We can also compare this result with solution found from an # interior-point approach called via the "KWDual" function from the # REBayes package. (This requires installation of the MOSEK # optimization library as well as the REBayes package, so we have # made this step optional.) if (FALSE) { out.kwdual <- mixkwdual(L,w) print(mixobjective(L,out.kwdual$x,w),digits = 16) } ``` -------------------------------- ### Build pkgdown site for mixsqp Source: https://stephenslab.github.io/mixsqp/index.html Build the pkgdown website for mixsqp while connected to the Internet. Version 1.4.1 of pkgdown was used. ```R pkgdown::build_site() ``` -------------------------------- ### Load mixsqp package Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html Loads the mixsqp package so its functions can be used. ```R library(mixsqp) ``` -------------------------------- ### Compare SQP and IP objective values Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html Prints the objective values from the SQP and IP solutions and the difference between them, showing they are nearly identical. ```R cat(sprintf("Objective at SQP solution: %0.16f\n",fit.sqp$value,digits = 16)) cat(sprintf("Objective at IP solution: %0.16f\n",fit.ip$value,digits = 16)) cat(sprintf("Difference in objectives: %0.4e\n",fit.sqp$value - fit.ip$value)) ``` -------------------------------- ### Run IP algorithm with mixkwdual Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html Runs the IP method using the mixkwdual function from the REBayes package, timing the computation and printing the elapsed time. ```R timing <- system.time(fit.ip <- mixkwdual(L)) cat(sprintf("Computation took %0.2f seconds\n",timing["elapsed"])) ``` -------------------------------- ### mixsqp-package Source: https://stephenslab.github.io/mixsqp/reference/index.html mixsqp: Sequential Quadratic Programming for Fast Maximum-Likelihood Estimation of Mixture Proportions. ```APIDOC ## mixsqp-package ### Description mixsqp: Sequential Quadratic Programming for Fast Maximum-Likelihood Estimation of Mixture Proportions. ``` -------------------------------- ### Fit mixture model with mixsqp Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html Runs the SQP algorithm on the likelihood matrix. The output shows default settings and progress, including objective, dual residual, and number of nonzeros. Convergence is achieved when the dual residual is below conv.tol (default 1e-8). ```R fit.sqp <- mixsqp(L) # Running mix-SQP algorithm 0.3-17 on 1000 x 20 matrix # convergence tol. (SQP): 1.0e-08 # conv. tol. (active-set): 1.0e-10 # zero threshold (solution): 1.0e-08 # zero thresh. (search dir.): 1.0e-10 # l.s. sufficient decrease: 1.0e-02 # step size reduction factor: 7.5e-01 # minimum step size: 1.0e-08 # max. iter (SQP): 1000 # max. iter (active-set): 20 # number of EM iterations: 10 # Computing SVD of 1000 x 20 matrix. # Matrix is not low-rank; falling back to full matrix. # iter objective max(rdual) nnz stepsize max.diff nqp nls # 1 +6.825854577e-01 -- EM -- 20 1.00e+00 3.43e-02 -- -- # 2 +6.608901451e-01 -- EM -- 20 1.00e+00 1.12e-02 -- -- # 3 +6.501638059e-01 -- EM -- 20 1.00e+00 8.83e-03 -- -- # 4 +6.441429919e-01 -- EM -- 20 1.00e+00 7.64e-03 -- -- # 5 +6.405380236e-01 -- EM -- 20 1.00e+00 6.44e-03 -- -- # 6 +6.382624101e-01 -- EM -- 20 1.00e+00 5.36e-03 -- -- # 7 +6.367521107e-01 -- EM -- 20 1.00e+00 4.46e-03 -- -- # 8 +6.357010188e-01 -- EM -- 20 1.00e+00 3.75e-03 -- -- # 9 +6.349367202e-01 -- EM -- 20 1.00e+00 3.18e-03 -- -- # 10 +6.343585101e-01 -- EM -- 20 1.00e+00 2.73e-03 -- -- # 1 +6.343585101e-01 +2.161e-02 20 ------ ------ -- -- # 2 +6.331717425e-01 +8.779e-03 4 1.00e+00 4.36e-01 20 1 # 3 +6.282364710e-01 +8.715e-05 4 1.00e+00 3.39e-03 2 1 # 4 +6.281978492e-01 +9.817e-09 4 1.00e+00 4.13e-05 2 1 # Optimization took 0.01 seconds. # Convergence criteria met---optimal solution found. ``` -------------------------------- ### mixsqp_control_default() Source: https://stephenslab.github.io/mixsqp/reference/index.html Default control parameters for mixsqp. ```APIDOC ## mixsqp_control_default() ### Description Default control parameters for mixsqp. ### Usage mixsqp_control_default() ``` -------------------------------- ### Plot optimal solution for tacks data Source: https://stephenslab.github.io/mixsqp/reference/tacks.html Loads the tacks dataset and plots the optimal solution, which is noted for being extremely sparse. ```R data(tacks) plot(tacks$x,type = "l",col = "royalblue") ``` -------------------------------- ### mixsqp Source: https://stephenslab.github.io/mixsqp/reference/mixsqp.html Solves the convex optimization problem of minimizing a weighted sum of negative log-likelihoods subject to a simplex constraint, using a sequential quadratic programming (SQP) algorithm. ```APIDOC ## mixsqp ### Description Solves the convex optimization problem of minimizing a weighted sum of negative log-likelihoods subject to a simplex constraint, using a sequential quadratic programming (SQP) algorithm. ### Method Function call ### Endpoint mixsqp(L, w, control = list()) ### Parameters #### Arguments - **L** (matrix) - Required - A matrix with n rows and m columns containing only non-negative entries. - **w** (numeric vector) - Required - A vector of non-negative weights of length n. - **control** (list) - Optional - A list of control parameters that override the default optimization algorithm settings. See Details for available components. ### Return Value A list with the following elements: - **x** (numeric vector) - The solution to the convex optimization problem if the SQP algorithm converges; otherwise, the best estimate of the solution achieved by the algorithm. Note that if the algorithm fails to converge, the entries of `x` may not sum to 1. - **value** (numeric) - The value of the objective function, f(x), at `x`. - **status** (character) - A character string describing the status of the algorithm upon termination. - **progress** (data.frame) - A data frame containing detailed information about the algorithm's progress, with one row per SQP iteration. See the description of the `verbose` control parameter in Details for an explanation of the columns. ### Details The `mixsqp` function solves the following optimization problem. Let L be a matrix with n rows and m columns containing only non-negative entries, and let w=(w1,…,wn) be a vector of non-negative "weights". `mixsqp` computes the value of vector x=(x1,…,xm) minimizing the objective function f(x)=−∑j=1n wj log(∑k=1m Ljk xk), subject to the constraint that x lie within the simplex; that is, the entries of x are non-negative and sum to 1. Implicitly, there is an additional constraint L∗x>0 to ensure the objective has a finite value. If all weights are equal, solving this optimization problem corresponds to finding the maximum-likelihood estimate of the mixture proportions x given n independent data points drawn from a mixture model with m components. In this case, Ljk is the likelihood for mixture component k and data point j. The Expectation Maximization (EM) algorithm can be used to solve this optimization problem, but it is intolerably slow in many interesting cases, and mixsqp is much faster. A special feature of this optimization problem is that the gradient of the objective does not change with re-scaling; for example, if all the entries of matrix `L` are multiplied by 100, the gradient does not change. A practical benefit of this property is that the optimization algorithm will behave similarly irrespective of the scale of `L`; for example, the same value for the convergence tolerance `convtol.sqp` will have the same effect at different scales. A related feature is that the solution to the optimization problem is invariant to rescaling the rows of `L`; for example, the solution will remain the same after all the entries in a row of `L` are multiplied by 10. A simple normalization scheme divides each row by the largest entry in the row so that all entries of `L` are at most 1: `L <- L / apply(L,1,max)` Occasionally, it can be helpful to normalize the rows when some of the entries are unusually large or unusually small. This can help to avoid numerical overflow or underflow errors. The SQP algorithm is implemented using the Armadillo C++ linear algebra library, which can automatically take advantage of multithreaded matrix computations to speed up `mixsqp` for large `L` matrices, but only when R has been configured with a multithreaded BLAS/LAPACK library (e.g., OpenBLAS). The `control` argument is a list in which any of the following named components will override the default optimization algorithm settings (as they are defined by `mixsqp_control_default`): - **normalize.rows** (logical) - When `normalize.rows = TRUE`, the rows of the data matrix `L` are automatically scaled so that the largest entry in each row is 1. This is the recommended setting for better stability of the optimization. When `log = TRUE`, this setting is ignored because the rows are already normalized. Note that the objective is computed on the original (unnormalized) matrix to make the results easier to interpret. - **force.sparse.init** (logical) - The computational complexity of the SQP updates grows rapidly by the number of nonzeros in the current solution estimate. When `force.sparse.init = TRUE`, and there are more than 20 nonzeros in the solution before beginning the SQP updates, only the 20 largest entries are kept, the rest are forced to zero. - **tol.svd** (numeric) - Setting used to determine rank of truncated SVD approximation for L. The rank of the truncated singular value decomposition is determined by the number of singular values surpassing `tol.svd`. When `tol.svd = 0` or when `L` has 4 or fewer columns, all computations are performed using full L matrix. ``` -------------------------------- ### mixobjective() Source: https://stephenslab.github.io/mixsqp/reference/index.html Compute objective optimized by mixsqp. ```APIDOC ## mixobjective() ### Description Compute objective optimized by mixsqp. ### Usage mixobjective() ``` -------------------------------- ### Generate small data set Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html Generates a small data set using simulatemixdata and prints its dimensions. ```R L <- simulatemixdata(1000,20)$L dim(L) ``` -------------------------------- ### mixsqp() Source: https://stephenslab.github.io/mixsqp/reference/index.html Maximum-likelihood estimation of mixture proportions using SQP. ```APIDOC ## mixsqp() ### Description Maximum-likelihood estimation of mixture proportions using SQP. ### Usage mixsqp() ``` -------------------------------- ### Show session information Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html Displays the R session information, including the version of R, platform, and packages used to generate the vignette results. ```R sessionInfo() # R version 3.6.2 (2019-12-12) # Platform: x86_64-apple-darwin15.6.0 (64-bit) # Running under: macOS Catalina 10.15.2 # # Matrix products: default # BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib # LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib # # locale: # [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 # # attached base packages: # [1] stats graphics grDevices utils datasets methods base # # other attached packages: # [1] mixsqp_0.3-18 # # loaded via a namespace (and not attached): # [1] Rcpp_1.0.3 knitr_1.26 magrittr_1.5 REBayes_1.8 # [5] MASS_7.3-51.4 lattice_0.20-38 R6_2.4.1 rlang_0.4.2 # [9] stringr_1.4.0 tools_3.6.2 grid_3.6.2 xfun_0.11 # [13] irlba_2.3.3 htmltools_0.4.0 yaml_2.2.0 assertthat_0.2.1 # [17] digest_0.6.23 rprojroot_1.3-2 pkgdown_1.4.1 crayon_1.3.4 # [21] Matrix_1.2-18 fs_1.3.1 memoise_1.1.0 evaluate_0.14 # [25] rmarkdown_2.0 stringi_1.4.3 compiler_3.6.2 desc_1.2.0 # [29] backports_1.1.5 Rmosek_9.0.96 ``` -------------------------------- ### simulatemixdata() Source: https://stephenslab.github.io/mixsqp/reference/index.html Create likelihood matrix from simulated data set. ```APIDOC ## simulatemixdata() ### Description Create likelihood matrix from simulated data set. ### Usage simulatemixdata() ``` -------------------------------- ### Set random seed Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html Initializes the sequence of pseudorandom numbers for reproducible results. ```R set.seed(1) ``` -------------------------------- ### mixobjective(L, x, w = rep(1, nrow(L))) Source: https://stephenslab.github.io/mixsqp/reference/mixobjective.html Computes the value of the objective function optimized by the mix-SQP algorithm at the point x. The objective is defined by the matrix L and optional weights w. Returns Inf if any entry of L %*% x is non-positive. ```APIDOC ## mixobjective(L, x, w = rep(1, nrow(L))) ### Description Computes the value of the objective function optimized by the mix-SQP algorithm at the point x. The objective is defined by the matrix L and optional weights w. Returns Inf if any entry of L %*% x is non-positive. ### Usage mixobjective(L, x, w = rep(1, nrow(L))) ### Arguments - **L** (matrix) - Required - Matrix specifying the optimization problem. In the context of mixture-model fitting, L[j,k] should be the value of the kth mixture component density at the jth data point. L should be a numeric matrix with at least two columns, with all entries being non-negative and finite (and not missing). No column should be entirely zeros. For large matrices, it is preferable that the matrix is stored in double-precision; see storage.mode. - **x** (numeric) - Required - The point at which the objective is evaluated; see argument x0 in mixsqp for details. - **w** (numeric vector) - Optional - An optional numeric vector, with one entry for each row of L, specifying the "weights" associated with the rows of L. All weights must be finite, non-negative and not missing. Internally, the weights are normalized to sum to 1, which does not change the problem, but does change the value of the objective function reported. By default, all weights are equal. ### Value The value of the objective at x. If any entry of L %*% x is less than or equal to zero, Inf is returned. ### See Also mixsqp ``` -------------------------------- ### Plot objective progress from mixsqp fit Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html Uses the progress data stored in fit.sqp$progress to plot the objective value at each SQP iteration. ```R numiter <- nrow(fit.sqp$progress) plot(1:numiter,fit.sqp$progress$objective,type = "b", pch = 20,lwd = 2,xlab = "SQP iteration", ylab = "objective",xaxp = c(1,numiter,numiter - 1)) ``` -------------------------------- ### mixkwdual() Source: https://stephenslab.github.io/mixsqp/reference/index.html Maximum-likelihood estimation of mixture proportions using MOSEK. ```APIDOC ## mixkwdual() ### Description Maximum-likelihood estimation of mixture proportions using MOSEK. ### Usage mixkwdual() ``` -------------------------------- ### Simulate mixture data with simulatemixdata Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html Creates an n×m conditional likelihood matrix for a mixture of zero-centered normals. By default, rows are normalized so the maximum entry in each row is 1. ```R L <- simulatemixdata(1000,20)$L dim(L) # [1] 1000 20 ``` -------------------------------- ### Simulate larger mixture data set Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html Creates a larger likelihood matrix with 100,000 samples and 100 mixture components for performance comparison. ```R L <- simulatemixdata(1e5,100)$L dim(L) # [1] 100000 100 ``` -------------------------------- ### Time SQP fit on large data set Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html Measures the computation time for fitting the SQP model on the large data set and prints the elapsed time. ```R timing <- system.time(fit.sqp <- mixsqp(L)) cat(sprintf("Computation took %0.2f seconds\n",timing["elapsed"])) ``` -------------------------------- ### simulatemixdata function signature Source: https://stephenslab.github.io/mixsqp/reference/simulatemixdata.html The function signature for `simulatemixdata`, which simulates data and computes a likelihood matrix. It takes parameters for sample size, number of mixture components, simulation type, log transformation, and row normalization. ```R simulatemixdata( n, m, simtype = c("n", "nt"), log = FALSE, normalize.rows = !log ) ``` === COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.