User Guide

01Introduction #

What problem does DLSM solve?

Traditional learning analytics treats engagement as a static quantity: take a window of LMS activity, compute some counts, and predict an outcome. This works when learning is cohort-based and regular. It fails in self-paced learning, where:

• Students register at any time.
• LMS activity is irregular and sparse.
• Much learning happens offline (textbook reading, word processing).
• A student can pause for weeks and still return.

The Dynamic Latent-State Model (DLSM) instead treats engagement as a latent, dynamic, multi-dimensional state that evolves continuously and generates observations through multiple noisy channels. The library implements this model end-to-end: simulation, filtering, parameter estimation, and evaluation.

Who is this guide for?

• Learning analytics researchers who want to model engagement as a dynamic state.
• Educational data scientists who need Bayesian state-space inference on LMS data.
• Developers building engagement-aware interventions.

What you can do with the library

02Installation #

From PyPI

pip install dlsm

From source

git clone https://github.com/hongxueharriswang/dlsm-lib.git
cd dlsm-lib
pip install -e .

Optional extras

The core library needs only NumPy and SciPy. Visualization and progress bars are optional:

pip install "dlsm[viz]"        # matplotlib
pip install "dlsm[progress]"   # tqdm
pip install "dlsm[all]"        # everything

Verifying the installation

import dlsm
print(dlsm.__version__)          # 0.1.0

from dlsm import DLSM, DLSParams
model = DLSM(DLSParams.default())
print(model.params.alpha_O)      # 0.3

If that runs without error, you are ready.

Development environment

pip install -r requirements-dev.txt
pip install -r requirements-lock.txt   # exact pinned versions
pytest tests/ -v

03Your First Five Minutes #

A complete end-to-end run in twenty lines:

import numpy as np
from dlsm import DLSM, DLSParams, ParticleFilter, Simulator

# 1. Set up the model with default parameters
params = DLSParams.default()
model = DLSM(params)

# 2. Simulate a cohort of 10 students over 30 time steps
sim = Simulator(model, seed=42)
data = sim.simulate(n_students=10, n_timesteps=30)

# 3. Filter the first student's observations
pf = ParticleFilter(model, n_particles=500, seed=0)
result = pf.filter(
    observations=data.observations[0],
    events=data.events[0],
)

# 4. Inspect the result
print("Filtered states shape:", result.mean_states.shape)   # (30, 6)
print("Log-likelihood:      ", round(result.log_likelihood, 2))
print("Mean ESS:            ", round(result.ess.mean(), 1))
print("Final status:        ", result.mean_status[-1])

# 5. Compare with ground truth
true_engagement = data.states[0, :, :4]
est_engagement  = result.mean_states[:, :4]
rmse = np.sqrt(np.mean((true_engagement - est_engagement) ** 2))
print("Engagement RMSE:     ", round(rmse, 4))

04The Latent State #

State vector

At each time t, a student's latent state is:

Xₜ = ( Sₜ , Eₜ , Pₜ , Rₜ )

Engagement dimensions

These map onto established engagement theory (Fredricks et al., 2004; Reeve & Tseng, 2011). The hierarchical constraint eO ≥ eP ≥ eC reflects the observation that more intensive engagement typically requires a foundation of basic engagement.

Statuses

Status is an IntEnum:

from dlsm.core.state import Status

Status.REGISTERED    # 0
Status.PROGRESSING   # 1
Status.STOPPED_OUT   # 2
Status.WITHDRAWN     # 3   <- absorbing
Status.COMPLETED     # 4   <- absorbing

WITHDRAWN and COMPLETED are absorbing: once a student enters, they cannot leave.

Risk

R is a continuous latent variable representing the propensity for adverse outcomes. It is not the same as the discrete status:

• Status is observed (administrative records).
• Risk is latent and inferred from engagement, progress, and inactivity.
• Risk influences the timing of status transitions.
• A student can be "Progressing" with high risk, or "Stopped-Out" with low risk.

The ordered-cone constraint

The constraint eO ≥ eP ≥ eC is enforced by projection onto the set:

{ e ∈ [0,1]⁴ : eO ≥ eP ≥ eC }

You rarely need to call the projection directly, but it is exposed:

from dlsm.core.constraints import project_ordered_cone
import numpy as np

e = np.array([[0.5, 0.7, 0.3, 0.2]])   # violates the constraint
p = project_ordered_cone(e)
print(p)                                # [[0.6 0.6 0.3 0.2]]

The DLSMState class

from dlsm.core.state import DLSMState, Status
import numpy as np

s = DLSMState(
    status=Status.PROGRESSING,
    continuous=np.array([0.6, 0.4, 0.2, 0.5, 0.3, 0.15]),
)
print(s.engagement)   # [0.6 0.4 0.2 0.5]
print(s.progress)     # 0.3
print(s.risk)         # 0.15
print(s.e_O)          # 0.6

05Parameters and What They Mean #

Creating a parameter set

from dlsm import DLSParams

params = DLSParams.default()             # sensible defaults
params = DLSParams(alpha_O=0.4)          # override one field
params = DLSParams(
    alpha_O=0.35, lambda_O=0.12,
    gamma=0.25, omega_3=0.08,
)                                        # override several

Engagement growth rates (α)

Control how fast engagement rises when the corresponding event occurs.

All are strictly positive. Typical range: 0.1 – 0.5.

Engagement decay rates (λ)

Control how fast engagement decays in the absence of events. All are strictly positive. Typical range: 0.05 – 0.3.

Progress dynamics (β)

Risk dynamics (γ, ω)

Process noise (q_diag)

Length-6 vector of process noise variances, one per continuous component:

params.q_diag = np.array([0.02, 0.02, 0.02, 0.02, 0.01, 0.02])

Serializing parameters

v = params.as_vector()                 # 1-D numpy array of 17 floats
print(v.shape)                         # (17,)

params2 = DLSParams.from_vector(v, template=params)

The 17-vector contains only the rate parameters (α, λ, β, γ, ω). Non-rate fields are taken from the template. This is what PMCMC uses internally.

The parameter names

from dlsm.inference.pmcmc import PMCMC
print(PMCMC.PARAM_NAMES)
# ['alpha_O', 'alpha_P', 'alpha_C', 'alpha_SR',
#  'lambda_O', 'lambda_P', 'lambda_C', 'lambda_SR',
#  'beta_0', 'beta_1', 'beta_2', 'beta_3',
#  'gamma',
#  'omega_0', 'omega_1', 'omega_2', 'omega_3']

06Observations and Measurement Channels #

The Observation class

from dlsm.core.observations import Observation
import numpy as np

obs = Observation(
    counts=np.array([3, 1, 0, 2]),       # shape (K_count,)
    binaries=np.array([1]),              # shape (K_binary,)
    ordinals=np.array([2]),              # shape (K_ordinal,)
    events=np.array([1, 1, 0, 1]),       # shape (4,)
)

Any channel can be None, which means the observation is missing for that channel at that time step. The particle filter handles missing channels gracefully.

Observation channels

Poisson (counts)

Used for LMS activity counts. Rate is exponential in the latent state:

y_k ~ Poisson( exp( w_kᵀ X ) )

Bernoulli (binary)

Used for submission indicators, login flags, etc.:

y_k ~ Bernoulli( σ( w_kᵀ X ) )

Ordered probit (ordinal)

Used for self-reported engagement, Likert responses, etc.:

y_k ~ OrderedProbit( w_kᵀ X , θ_k )

The default channel layout

Missing data

obs_partial = Observation(
    counts=np.array([3, 1, 0, 2]),
    binaries=None,       # missing
    ordinals=None,       # missing
)

The particle filter propagates the state forward without updating on the missing channels.

07Simulating Data #

from dlsm import DLSM, DLSParams, Simulator

model = DLSM(DLSParams.default())
sim = Simulator(model, seed=42)
data = sim.simulate(
    n_students=50,
    n_timesteps=60,
    event_prob=0.3,
)

What you get

data.states           # (50, 60, 6)   continuous states
data.statuses         # (50, 60)      discrete statuses
data.observations     # list of 50 lists of 60 Observation objects
data.events           # (50, 60, 4)   binary event indicators

Custom initial state

import numpy as np

x0 = np.array([0.3, 0.2, 0.1, 0.4, 0.0, 0.5])
data = sim.simulate(n_students=20, n_timesteps=60, x0=x0)

Custom event process

rng = np.random.default_rng(0)
events = (rng.random((20, 60, 4)) < 0.4).astype(int)
events[:, :10, :] = 0    # first 10 steps: no events

states, statuses, obs = model.simulate_trajectory(
    n_timesteps=60, events=events[0], rng=rng,
)

Simulating a specific scenario

x0_engaged = np.array([0.8, 0.6, 0.4, 0.7, 0.3, 0.1])
data = sim.simulate(
    n_students=1, n_timesteps=90,
    event_prob=0.05, x0=x0_engaged,
)
# Engagement decays over time

08Filtering: Tracking Latent States #

from dlsm import DLSM, DLSParams, ParticleFilter

model = DLSM(DLSParams.default())
pf = ParticleFilter(model, n_particles=1000, seed=0)

result = pf.filter(
    observations=data.observations[0],
    events=data.events[0],
)

The FilterResult object

Computing credible intervals

import numpy as np

mean = result.mean_states          # (T, 6)
std  = np.sqrt(result.var_states)  # (T, 6)

lower = mean - 1.96 * std
upper = mean + 1.96 * std

Storing full particles

result = pf.filter(
    observations=data.observations[0],
    events=data.events[0],
    store_particles=True,
)

particles = result.particles      # (T, N, 6)
weights   = result.weights        # (T, N)

Choosing the number of particles

09Estimating Parameters with PMCMC #

from dlsm import DLSM, DLSParams, PMCMC

model = DLSM(DLSParams.default())
pmcmc = PMCMC(
    model,
    n_particles=2000,
    n_iterations=10000,
    burn_in=2000,
    thin=10,
    proposal_scale=0.02,
    seed=0,
)

result = pmcmc.fit(
    observations=data.observations[0],
    events=data.events[0],
)

The PMCMCResult object

Reading the output

from dlsm.evaluation import metrics

table = metrics.parameter_recovery_table(
    true_params=DLSParams.default().as_vector(),
    posterior_samples=result.chain,
    param_names=result.param_names,
)
for row in table:
    print(
        f"{row['parameter']:10s}  "
        f"true={row['true']:+.3f}  "
        f"post={row['posterior_mean']:+.3f}  "
        f"std={row['posterior_std']:.3f}  "
        f"covered={row['covered']}"
    )

Tuning the proposal scale

Target acceptance rate: 0.2 – 0.4.

Multi-chain runs

chains = []
for seed in (0, 1, 2):
    pmcmc = PMCMC(model, seed=seed, ...)
    chains.append(pmcmc.fit(obs, events).chain)

stacked = np.stack(chains)     # (3, n_keep, 17)

from dlsm.inference.diagnostics import gelman_rubin, effective_sample_size
rhat = gelman_rubin(stacked)
ess  = effective_sample_size(stacked)

print("Max R-hat:", rhat.max())          # target < 1.1
print("Min ESS:  ", ess.min())           # target > 500

10Working with Real Data #

Data requirements

To apply the DLSM to real data, you need per-student time series of:

1. Registration and status records
2. LMS activity counts (daily or weekly views, downloads, plays)
3. Assessment records (submission dates, scores)
4. Communication logs (tutor contacts, forum posts)
5. Extension requests (dates and outcomes)

Building Observation objects from a DataFrame

from dlsm.core.observations import Observation
import numpy as np
import pandas as pd

def df_to_observations(df: pd.DataFrame, freq="W") -> list[Observation]:
    df = df.set_index("date").sort_index()
    agg = df.resample(freq).agg({
        "page_views": "sum",
        "video_plays": "sum",
        "forum_reads": "sum",
        "forum_posts": "sum",
        "assignments_submitted": "sum",
        "self_report": "mean",
    }).fillna(0)

    obs_list = []
    for _, row in agg.iterrows():
        counts = np.array([
            int(row["page_views"] + row["video_plays"]),
            int(row["forum_reads"]),
            int(row["forum_posts"]),
            0,
        ])
        binaries = np.array([int(row["assignments_submitted"] > 0)])
        ordinals = np.array([int(row["self_report"])]) \
                    if row["self_report"] > 0 else None
        obs_list.append(Observation(
            counts=counts, binaries=binaries, ordinals=ordinals,
        ))
    return obs_list

Batch filtering a cohort

from joblib import Parallel, delayed

def filter_one(i):
    pf = ParticleFilter(model, n_particles=500, seed=1000 + i)
    return pf.filter(cohort_obs[i], events=events[i])

results = Parallel(n_jobs=-1)(
    delayed(filter_one)(i) for i in range(n_students)
)

Real-time risk estimation

pf = ParticleFilter(fitted_model, n_particles=1000, seed=0)
result = pf.filter(observations=new_student_obs)

latest_risk = result.mean_states[-1, 5]
latest_risk_std = np.sqrt(result.var_states[-1, 5])
print(f"Risk: {latest_risk:.2f} ± {latest_risk_std:.2f}")

11Custom Observation Channels #

Redefining the channel layout

import numpy as np
from dlsm import DLSParams

params = DLSParams.default()

# Rows: channels. Columns: [eO, eP, eC, eSR, P, R]
params.w_count = np.array([
    [1.0, 0.0, 0.0, 0.0, 0.0, 0.0],   # page_views -> eO
    [0.7, 0.5, 0.0, 0.0, 0.0, 0.0],   # video_plays -> eO + eP
])

params.w_binary = np.array([
    [0.0, 0.0, 0.0, 0.0, 1.0, 0.0],   # submitted -> P
])

params.w_ordinal = np.array([
    [0.0, 0.0, 0.0, 1.0, 0.0, 0.0],   # self-report -> eSR
])

params.ordinal_thresholds = np.array([-0.5, 0.5, 1.5])   # 4 levels

Adding a new channel type

from dlsm.core.observations import ObservationModel
import numpy as np

class ExtendedObservationModel(ObservationModel):

    def __init__(self, params, w_gauss=None, sigma_gauss=0.3):
        super().__init__(params)
        self.w_gauss = w_gauss if w_gauss is not None else np.zeros((1, 6))
        self.sigma_gauss = sigma_gauss

    def log_likelihood_gaussian(self, x, y):
        mu = x @ self.w_gauss.T
        resid = y - mu
        return -0.5 * (resid / self.sigma_gauss) ** 2 \
               - np.log(self.sigma_gauss) \
               - 0.5 * np.log(2 * np.pi)

    def log_likelihood(self, x, counts=None, binaries=None,
                       ordinals=None, gaussian=None):
        ll = super().log_likelihood(x, counts, binaries, ordinals)
        if gaussian is not None:
            ll = ll + self.log_likelihood_gaussian(x, gaussian).sum(axis=1)
        return ll

12Custom Drift Functions #

from dlsm.core.dynamics import Dynamics
import numpy as np

class SaturatedDynamics(Dynamics):

    def drift(self, x, events, dt=None):
        base = super().drift(x, events, dt=dt)
        # Add a saturating term that pulls eO back toward 0.6
        base[:, 0] += -0.1 * (x[:, 0] - 0.6)
        return base

model = DLSM(params)
model.dynamics = SaturatedDynamics(params)

Time-varying dt

The step method accepts a dt argument for irregularly-spaced observations:

dt = 5.0    # five days elapsed since last observation
x_next = dynamics.step(x, events, dt=dt, rng=rng)

The drift is scaled by dt and the process noise by √dt, preserving the continuous-time interpretation.

13Diagnostics and Model Checking #

Particle-filter diagnostics

result = pf.filter(obs, events=events)
print("Mean ESS:   ", result.ess.mean())
print("Min ESS:    ", result.ess.min())
print("Below 30%:  ", (result.ess < 0.3 * pf.n_particles).mean())

If mean(ess) < 0.3 * n_particles, increase n_particles or investigate the model.

PMCMC diagnostics

Trace plots

from dlsm.viz import plots

fig = plots.plot_posterior(
    result.chain,
    result.param_names,
    true_params=params.as_vector(),
)
fig.savefig("posteriors.png", dpi=150)

Gelman–Rubin

from dlsm.inference.diagnostics import gelman_rubin

# Requires ≥ 2 chains
rhat = gelman_rubin(np.stack([chain1, chain2, chain3]))
# R-hat < 1.1 for all parameters indicates convergence

Effective sample size

from dlsm.inference.diagnostics import effective_sample_size

ess = effective_sample_size(np.stack([chain1, chain2]))
# Aim for ESS > 500 per parameter

Posterior predictive checks

theta_hat = result.posterior_mean()
params_hat = DLSParams.from_vector(theta_hat, template=params)
model_hat = DLSM(params_hat)

from dlsm import Simulator
sim = Simulator(model_hat, seed=0)
sim_data = sim.simulate(n_students=50, n_timesteps=60)

print("Observed mean eO: ", data.states[..., 0].mean())
print("Simulated mean eO:", sim_data.states[..., 0].mean())

14Performance and Scaling #

Profiling a single filtering pass

import time

t0 = time.perf_counter()
pf.filter(obs, events=events)
elapsed = time.perf_counter() - t0
print(f"Filtering: {elapsed:.3f} s")

Parallelizing across students

from joblib import Parallel, delayed

def _filter_one(i, obs, events, model, seed):
    pf = ParticleFilter(model, n_particles=500, seed=seed + i)
    return pf.filter(obs, events=events)

results = Parallel(n_jobs=-1)(
    delayed(_filter_one)(i, data.observations[i], data.events[i], model, 0)
    for i in range(len(data.observations))
)

Memory considerations

For a trajectory of T = 100 with N_p = 1000 particles:

• Without store_particles: ~50 MB.
• With store_particles=True: ~480 MB.

Only request stored particles when you need them.

15Reproducibility #

Seeds everywhere

sim = Simulator(model, seed=42)
pf = ParticleFilter(model, n_particles=1000, seed=0)
pmcmc = PMCMC(model, seed=0, ...)

Same seed → identical results:

a = Simulator(model, seed=42).simulate(n_students=5, n_timesteps=20)
b = Simulator(model, seed=42).simulate(n_students=5, n_timesteps=20)
assert np.allclose(a.states, b.states)

Locking the environment

pip install -r requirements-lock.txt

Logging a run

import json
from datetime import datetime

run_metadata = {
    "date": datetime.now().isoformat(),
    "model_params": params.to_dict(),
    "n_particles": pf.n_particles,
    "n_students": len(data.observations),
    "n_timesteps": len(data.observations[0]),
}
with open("run_metadata.json", "w") as f:
    json.dump(run_metadata, f, indent=2)

16API Cheat Sheet #

Core

from dlsm import (
    DLSParams,        # parameter container
    DLSM,             # model
    DLSMState,        # single state
    Status,           # discrete status enum
    Dynamics,         # transition dynamics
    ObservationModel, # observation channels
    Observation,      # single observation
)

Inference

from dlsm import ParticleFilter, PMCMC
from dlsm.inference import (
    gelman_rubin,
    effective_sample_size,
)

Simulation

from dlsm import Simulator

Evaluation & visualization

from dlsm.evaluation import metrics
from dlsm.viz import plots

Common patterns

The 17-vector

v = params.as_vector()
# v[0:4]    α (alpha_O, alpha_P, alpha_C, alpha_SR)
# v[4:8]    λ (lambda_O, lambda_P, lambda_C, lambda_SR)
# v[8:12]   β (beta_0, beta_1, beta_2, beta_3)
# v[12]     γ (gamma)
# v[13:17]  ω (omega_0, omega_1, omega_2, omega_3)

The 6-vector (continuous state)

x = state.continuous
# x[0]  eO
# x[1]  eP
# x[2]  eC
# x[3]  eSR
# x[4]  P
# x[5]  R

17Troubleshooting & FAQ #

Q1: The particle filter ESS collapses to 1 immediately

Symptom: result.ess drops to 1 after the first few steps.

• Too few particles — increase n_particles (try 2000+).
• Observation model mismatch — rescale counts or adjust w_count.
• Process noise too small — increase q_diag.

Q2: PMCMC acceptance rate is near 0

Fix: decrease proposal_scale:

pmcmc = PMCMC(model, proposal_scale=0.005, ...)

Q3: PMCMC acceptance rate is near 1

Fix: increase proposal_scale:

pmcmc = PMCMC(model, proposal_scale=0.05, ...)

Q4: Parameter recovery is poor for some parameters

Diagnosis: check posterior correlations.

corr = np.corrcoef(result.chain, rowvar=False)
# Look for |corr| > 0.8 pairs

Q5: The ordered-cone projection changes my estimates too much

raw = result.mean_states[:, :3]
from dlsm.core.constraints import project_ordered_cone
projected = project_ordered_cone(raw)
violations = np.abs(raw - projected).sum(axis=1) > 1e-6
print(f"Projection active at {violations.mean():.1%} of time steps")

Q6: I have missing timestamps in my data

obs = Observation(
    counts=np.zeros(K_count, dtype=int),
    binaries=np.zeros(K_binary, dtype=int),
    ordinals=None,           # truly missing, not zero
)

Q7: ValueError: q_diag must have length 6

params.q_diag = np.array([0.02, 0.02, 0.02, 0.02, 0.01, 0.02])

Q8: ValueError: w_status must have shape (5, 5, 6)

from dlsm.core.parameters import _default_status_coefficients
params.w_status = _default_status_coefficients()

Q9: How do I know if my model is well-specified?

1. ESS stays high — no degeneracy.
2. Posterior predictive checks match — simulated data resembles observed.
3. Residuals are unstructured — no systematic patterns.
4. PMCMC converges — R-hat < 1.1.
5. Parameters are identified — no extreme posterior correlations.

18Best Practices #

Model development

1. Start with the default parameters. Fit them to your data before modifying the model.
2. Prototype with simulation. Verify the model recovers what you put in before applying to real data.
3. Inspect the posterior. Never trust a single point estimate; always look at the chain.
4. Document your channel mapping. Which observation maps to which latent dimension is a modeling choice; record it.

Inference

1. Always check ESS. If ESS is low, increase particles before drawing conclusions.
2. Use multiple chains for PMCMC. A single chain cannot diagnose convergence.
3. Tune the proposal scale. Target 20–40% acceptance.
4. Thin the chain. Retain every 10th sample after burn-in.

Evaluation

1. Report calibration, not just accuracy. Brier score and calibration plots reveal more than AUC alone.
2. Compare against static baselines. Logistic regression and random forest are strong baselines.
3. Held-out evaluation. Train on some students, test on others — never evaluate on training data.
4. Report uncertainty. Always report credible intervals alongside point estimates.

Deployment

1. Version your model. Save params.to_dict() alongside any deployed predictions.
2. Log the seed. Even for production, log the RNG seed for incident debugging.
3. Monitor ESS in production. If ESS collapses, downstream predictions are unreliable.
4. Warm-start with the last posterior. For real-time tracking, reuse the last particle cloud as the prior for the next step.