Before attempting to run this Julia code, make sure your computer is properly set up to run this code by following the setup instructions in the README of the git repository.
We begin by loading the Distributions, LowRankApprox and other packages, as well as some function definitions used in the code chunks below.
using Pkg
using Printf
using Random
using Distributions
using LinearAlgebra
using SparseArrays
using LowRankApprox
include("../code/datasim.jl");
include("../code/likelihood.jl");
include("../code/mixSQP.jl");
Next, initialize the sequence of pseudorandom numbers.
Random.seed!(1);
Let's start with a smaller example with 50,000 samples.
z = normtmixdatasim(round(Int,5e4));
Compute the $n \times k$ likelihood matrix for a mixture of zero-centered normals, with $k = 20$. Note that the rows of the likelihood matrix are normalized by default.
sd = autoselectmixsd(z,nv = 20);
L = normlikmatrix(z,sd = sd);
size(L)
First we run the mix-SQP algorithm once to precompile the function.
out = mixSQP(L,maxiter = 10,verbose = false);
Observe that only a small number of iterations is needed to converge to the solution of the constrained optimization problem.
k = size(L,2);
x0 = ones(k)/k;
@time out = mixSQP(L,x = x0);
Next, let's see what happens when we use the SQP algorithm to fit a mixture model to a much larger data set.
Random.seed!(1);
z = normtmixdatasim(round(Int,1e5));
As before, we compute the $n \times k$ likelihood matrix for a mixture of zero-centered normals. This time, we use a finer grid of $k = 40$ normal densities.
sd = autoselectmixsd(z,nv = 40);
L = normlikmatrix(z,sd = sd);
size(L)
Even on this much larger data set, only a small number of iterations is needed to compute the solution.
k = size(L,2);
x0 = ones(k)/k;
@time out = mixSQP(L,x = x0);
With no low-rank approximation (lowrank = "none"), the algorithm still converges rapidly, even when using a very small correction factor eps = 1e-12.
@time out = mixSQP(L,x = x0,lowrank = "none",eps = 1e-12);
The section gives information about the computing environment used to generate the results contained in this notebook, including the version of Julia, and the versions of the Julia packages used here.
Pkg.status("Distributions");
Pkg.status("LowRankApprox");
versioninfo();