This package contains Python and C code for fitting the basis-function
mixture model of Newman (2025).

Files:
  emlpa.c: C code for the EM algorithm
  blpa.c: C code for the MCMC algorithm
  Makefile: Makefile for compiling the code on Unix-style systems such as Linux
  pylpa.py: Python module for using the code in Python
  em_wine.py: Example code reproducing the wine analysis from the paper
  em_pancreatic.py: Example code for EM algorithm on biomarkers data set
  mc_pancreatic.py: Example code for MCMC algorithm on biomarkers data set
  mc_big5.py: Example code for MCMC algorithm on personality test data set
  wine.txt: The wine data set, processed for use with this code
  pancreatic:txt: The pancreatic biomarkers data set, processed suitably
  big5.txt: The Big 5 personality test data, processed suitably
  license.txt: Licensing
  readme.txt: This readme file

Installation:

The code that actually performs the calculations is written in the C
programming language for speed.  (C code is typically about 100 times
faster than equivalent code in Python.)  To install this package you must
compile the code, which requires that you have a C compiler, such as gcc on
Linux or Visual C on Windows.  On Linux, compilation can be done with the
following commands, executed in the directory containing the code:

gcc -O2 -shared -fPIC -o libemlpa.so emlpa.c -lm -lgsl -lgslcblas
gcc -O2 -shared -fPIC -o libblpa.so blpa.c -lm -lgsl -lgslcblas

This will produce files libemlpa.so and libblpa.so, which should be placed
in the directory where you will be running the code.  A Linux make file is
also included for automated compilation.  You will also need to have the GNU Scientific Library (GSL) installed on your system; it can be freely downloaded from the web if you don't already have it.

On Windows or Mac there are equivalent compilation commands, but the author
does not currently have access to these systems, so the code has not been
tested there.


Use:

The software takes a data set in the form of an array of N rows and M
columns.  In the nomenclature of the paper, each row corresponds to an
"observation" and each column corresponds to an "item".  For instance, for
the personality test data in "big5.txt" each row is a respondent and each
column is one of the fifty questions on the test.

Alternatively, the software can take processed data in the form of the
quantities phi[i,j,t] described in the paper.  This allows advanced users
to employ their own basis function sets, without needing to modify the code
in any way.

The package provides a number of Python classes in the pylpa package as
described below.  Most users will probably want to use the LPA_Bernstein()
and LPA_Gamma() classes, which are described toward the end.

class State(k=None, g=None, h=None):

    Objects in this class represent assignments of components and "slots".
    The object has the following variables, whose values are initialized
    from the arguments of the same name:

    k: integer, number of components

    g: array of N integers with values 0...k-1, representing the
    component to which each observation belongs

    h: N-by-M array of integers with values 0...T-1, representing the
    slot assignment for each observation and each item

class LPA(phi, nt):

    This class provides methods for fitting the model and analyzing the
    results.

    phi: A 3d array of floats with dimension [N,M,tdim], where tdim is at
    least as great as the largest value in nt[].  This array contains the
    values of the quantities phi[i,j,t] defined in the paper.

    nt: An array of M integer elements giving the number T[j] of basis
    functions for each item j.


    The LPA object provides the following methods:

    LPA.em(k,maxiter=10000,tolerance=1e-10):

	Run the EM algorithm with k components.  Iteration of the algorithm
	ends when either the parameter values change by less than tolerance
	on each iteration or maxiter iterations have been completed,
	whichever comes sooner.

	return values:
	    q: array of estimated probabilities q[i,r] that observation i is
               a member of component r
            theta: array of estimated parameter values theta[r,j,t]

    LPA.burnin(nsweeps,initial=None):

	Run the Monte Carlo algorithm, performing the given number of Monte
	Carlo sweeps to "burn in" the calculation.  Markov chain Monte
	Carlo methods, such as the one used here, need to equilibrate,
	meaning the distribution of samples needs to converge from its
	initial state (which is usually a low-probability state) to the
	posterior distribution of the model.  This function does this.

	nsweeps: integer, the number of sweeps to perform to burn in the
	Monte Carlo.  A "sweep" is N Monte Carlo steps, where N is the
	number of observations.  A typical value for thorough burn-in would
        be nsweeps=2500.

	initial: a State() object containing the initial state for the
	Monte Carlo.  If this parameter is given, the function will use it
	as the initial state.  If it is not given, the Monte Carlo will
	start from a state with all observations in a single component and
	randomly assigned slots.  This is usually a perfectly good starting
	point, so it's safe to omit this parameter.

    LPA.sample(nsweeps,tsample,initial=None,slots=False,permute=True):

	Perform the given number of sweeps of the Monte Carlo, starting
	either from the given initial state or the prevailing current
	state, and record the component assignments and (optionally) the
	slot assignments for all observations at specified intervals.  You
	can do repeated calls to this function and each run will
	automatically continue from wherever the previous one left off, so
	you can easily pause in the middle of a calculation, or go back and
	take more samples if you don't have enough.

	The "current state" would normally be created by first running the
	function LPA.burnin() above: normal use is to call LPA.burnin()
	first, then LPA.sample().  However, an initial state can also be
	specified explicitly if preferred.

	nsweeps: integer, total number of sweeps to perform

	tsample: integer, sampling interval.  The program records the
	assignments every tsample sweeps.

	initial: a State() object containing the initial state for the
	Monte Carlo.  If this parameter is given, the function will use it
	as the initial state.  If it is not given, the Monte Carlo will
	start from the current state, which is usually the state from a
	previous call to LPA.burnin() or LPA.sample().

	slots: record the slot assignments as well as the component
	assignments.  Usually the slot assignments are not of interest, and
	they can take up a lot of space, so there's no point in recording
	them unless you need them.  The most obvious use case where you
	would need them is if you were going to calculate the mutual
	information using the LPA.mi() function below, which requires slot
	assignments.

	permute: bool, whether to permute component labels.  If True (the
        default) the component labels are randomly permuted before each
        sample is drawn.  This has no effect on the results, but eliminates
        the possibility of any apparent significance to specific component
        labels, which is formally impossible because of symmetry with
        respect to permutations (so-called label switching).  The
        permutation action is performed efficiently and incurs no
        significant performance penalty.  Omit this parameter unless you
        know what you are doing.

	Return value: numpy array of State() objects, one for each sample,
	each containing a number of components k, a component assignment
	array g, and (optionally) a slot assignment array h.  The length of
	the array of State() objects is, by definition, equal to
	nsweeps//tsamples.  For instance, if you ask for 100,000 sweeps and
	a sampling interval of 10, you will get an array of 10,000 samples.

    cmatrix(samples):

	Calculate the consensus matrix for a set of samples.  The consensus
	matrix is the matrix with elements C[i,j] equal to the fraction of
	times observations i and j were assigned to the same component.
	Return value is an N-by-N array of floats containing the matrix.

    consensus(samples,kvalue=0,cmatrix=None):

        Calculate a consensus clustering for a set of samples, using the L2
        method described in the paper.

	samples: a numpy array of State() objects containing the samples

	kvalue: If non-zero, the function will return a consensus
	clustering with exactly kvalue components.  If zero (the default),
	then the function will return the best consensus clustering it can
	find, regardless of the number of components.

	cmatrix: An N-by-N consensus matrix of floats.  If provided, the
	calculation will use this consensus matrix to compute the
	clustering.  If not provided, the code will calculate its own
	consensus matrix from the samples provided.

	return value: A State() object containing the sample state that
        best represents the entire set.

    mi(samples):

        Calculate the mutual information between component and slot
        assignments for a set of samples.  Returns a float array of mutual
        information values, one for each item in the data set.

	samples: a numpy array of State() objects containing the samples

    esttheta(sample,tolerance=1e-10):

        Estimate the theta parameters from a set of component assignments
        using the iterative method described in the paper.

	sample: a State() object containing the component assignments.  Any
	included slot assignments are ignored.

	tolerance: iteration stops once the parameter values change by less
	than this amount on each iteration.


class LPA_Bernstein(x,d):

    A child class of LPA() that implements the model using the Bernstein
    polynomial basis.

    x: Two-dimensional float array of the actual observational data, with
    values mapped to the interval [0,1].  If values fall outside this
    interval, behavior is undefined.

    d: Degree of polynomials to use for the basis functions.

    This class has the same methods as the parent LPA() class (see above).


class LPA_Gamma(x,d):

    A child class of LPA() that implements the model with the gamma basis.

    x: Two dimensional float array of the actual observational data,
    suitably scaled.  We find that scaling values for each item to have
    mean 0.5 works well.

    d: Maximum degree of polynomials to use for the basis functions.

    This class has the same methods as the parent LPA() class.


Seeding the random number generator:

The C code uses the GSL Mersenne twister random number generator, but this
generator is seeded from the standard Python random number generator.  So
if you want to seed the random numbers explicitly you should just call
random.seed() in your Python code as normal before the start of any
calculations.  Repeated runs seeded with the same seed will then produce
identical results, which can be useful for reproducibility, debugging, and
running on parallel architectures.  If no seed is specified, you will get
the default Python behavior, which on most systems means the program will
choose a random seed, for instance from the system clock or from some
operating-system-level source of randomness.


Example applications:

The files em_wine.py and em_pancreatic.py contain example code to reproduce
the EM analysis of the wide and pancreatic biomarker data sets from the
paper.  The files mc_pancreatic.py and mc_big5.py contain example code to
reproduce the MCMC analysis of biomarker and personality test data sets.