SVGD

phasic.svgd.SVGD(
    model,
    observed_data,
    prior=None,
    n_particles=None,
    n_iterations=1000,
    learning_rate=None,
    bandwidth='median_per_dim',
    theta_init=None,
    theta_dim=None,
    seed=None,
    verbose=False,
    progress=True,
    jit=None,
    parallel=None,
    n_devices=None,
    precompile=True,
    compilation_config=None,
    positive_params=True,
    param_transform=None,
    regularization=0.0,
    nr_moments=2,
    rewards=None,
    fixed=None,
    optimizer=None,
    preconditioner='auto',
    exposure=None,
    exposure_param_index=None,
    _validated=False,
)

Stein Variational Gradient Descent (SVGD) for Bayesian parameter inference.

This class provides an object-oriented interface for SVGD inference with automatic result storage and diagnostic plotting capabilities.

When to use `Graph.svgd` vs `SVGD` direct

Use :meth:phasic.Graph.svgd when you have a phasic Graph object — it handles model construction, reward validation, zero-inflation attachment for partial-coverage rewards, daisy-chain epoch wiring, tied parameters, and per-edge weight callbacks for you. The parameters that :meth:Graph.svgd and :class:SVGD share have identical names, defaults, and semantics.

Use SVGD(model=..., ...) direct only when you have a pre-built JAX-compatible model callable (e.g. an external PMF) and a fixed theta_dim. Parameters unique to :meth:Graph.svgd — discrete, joint_index, tied, callback, epoch_starts, daisy_chain_t_eval, daisy_chain_granularity, daisy_chain_probe_theta, daisy_chain_t_eval_tol, validate_rewards, return_history — are pre-model-construction concerns and have no meaning when SVGD is handed an opaque model callable.

Direct SVGD callers whose model carries a _source_graph attribute (stamped by the public PMF builders :meth:Graph.pmf_from_graph, :meth:Graph.pmf_and_moments_from_graph, :meth:Graph.pmf_and_moments_from_graph_multivariate, :meth:Graph.pmf_from_graph_joint_index) get the zero-inflated likelihood term auto-attached when rewards indicates partial coverage. Hand-built models without a source graph reference skip this safety net silently — see :mod:phasic.zero_inflation for the underlying helpers.

Parameters

model : callable: JAX-compatible parameterized model with signature: model(theta, data) -> values
observed_data : np.ndarray: Observed data points
prior : callable or list of Prior objects = None: Log prior function for parameters. Can be: - Single callable: prior(theta) -> scalar, applied to entire theta vector - List of Prior objects: One prior per parameter dimension. Use None for fixed parameters: prior=[GaussPrior(…), None, GaussPrior(…)] If None, uses standard normal prior: log p(theta) = -0.5 * sum(theta^2) With fixed parameters: When using a list of priors with the fixed parameter, you must provide None at indices corresponding to fixed parameters. This is validated at initialization. Example: prior=[GaussPrior(ci=[0,1]), None, GaussPrior(ci=[0,1])], fixed=[(1, 0.5)] # theta[1] fixed, prior[1] must be None
n_particles : int = is 20 times length of theta: Number of SVGD particles
n_iterations : int = 1000: Number of SVGD optimization steps
learning_rate : float, StepSizeSchedule, or None = None: SVGD step size. Can be: - None: Uses Adam optimizer with adaptive learning rates (default) - float: constant step size (uses fixed learning rate approach) - StepSizeSchedule object: dynamic step size schedule When None and no optimizer is provided, Adam is used as the default optimizer (unless regularization > 0, which falls back to fixed lr=0.001). Examples: ConstantStepSize(0.01), ExpStepSize(0.1, 0.01, 500.0)
bandwidth : str, float, or np.ndarray = 'median_per_dim': Kernel bandwidth selection. Can be: - ‘median_per_dim’: Per-dimension median heuristic (default). Uses a separate bandwidth for each parameter dimension, giving an anisotropic kernel that adapts to different parameter scales. - ‘median’: Scalar median heuristic (isotropic kernel) - float: Fixed scalar bandwidth value - np.ndarray: Fixed per-dimension bandwidth vector
theta_init : np.ndarray = None: Initial particle positions (n_particles, theta_dim)
theta_dim : int = None: Dimension of theta parameter vector (required if theta_init is None)
seed : int = 42: Random seed for reproducibility
verbose : bool = False: Print progress information
progress : bool = True: Display progress bar during optimization
jit : bool or None = None: Enable JIT compilation. If None, uses value from phasic.get_config().jit. If True, requires JAX to be available (raises PTDConfigError otherwise). JIT compilation provides significant speedup but adds initial compilation overhead.
parallel : str or None = None: Parallelization strategy: - ‘vmap’: Vectorize across particles (single device) - ‘pmap’: Parallelize across devices (uses multiple CPUs/GPUs) - ‘none’: No parallelization (sequential, useful for debugging) - None: Auto-select (pmap if multiple devices, vmap otherwise) Single-machine multi-CPU: Auto-selection uses pmap for multi-core parallelization. Multi-node SLURM: Call initialize_distributed() then set parallel=‘pmap’ explicitly.
n_devices : int or None = None: Number of devices to use for pmap. Only used when parallel=‘pmap’. If None, uses all available devices. Must be <= number of available JAX devices. See: jax.devices() to check available devices, or configure via PTDALG_CPUS environment variable before import.
precompile : bool = True: (Deprecated: use jit parameter instead) Precompile model and gradient functions for faster execution. Implies jit=True for backward compatibility. First run will take longer (compilation time) but subsequent iterations will be much faster. Compiled functions are cached in memory and on disk (~/.phasic_cache/).
compilation_config : CompilationConfig, dict, str, or Path = None: JAX compilation optimization configuration. Can be: - CompilationConfig object from phasic.CompilationConfig - dict with CompilationConfig parameters - str/Path to JSON config file - None (uses default balanced configuration) The configuration controls JAX/XLA compilation behavior including: - Persistent cache directory for cross-session caching - Optimization level (0-3) - Parallel compilation settings Examples: - Use preset: CompilationConfig.fast_compile() - Load from file: ‘my_config.json’ - Custom dict: {‘optimization_level’: 2, ‘cache_dir’: ‘/tmp/cache’}
positive_params : bool = True: If True, constrains parameters to positive domain using softplus transformation. SVGD operates in unconstrained space (can be negative) but results are transformed to positive values. DEFAULT=True because phase-type distribution parameters (rates) must be positive. Set to False only if you have a specific reason to allow negative parameters (e.g., regression coefficients, log-space parameterization).
param_transform : callable = None: Custom parameter transformation function. Overrides positive_params if provided. Should map unconstrained space to constrained space (e.g., lambda x: jax.nn.sigmoid(x) for parameters in [0,1]). Cannot be used together with positive_params=True.
regularization : float or RegularizationSchedule = 0.0: Moment-based regularization strength. Can be: - float: constant regularization (0.0 = no regularization, >0.0 = regularized SVGD) - RegularizationSchedule object: dynamic regularization schedule Examples: ConstantRegularization(1.0), ExpRegularization(5.0, 0.1, 500.0) If > 0.0, adds penalty term to match model moments to sample moments. Sample moments are computed from observed_data at initialization. Note: Using RegularizationSchedule disables gradient precompilation for flexibility, which may be slower than constant regularization but allows dynamic strategies.
nr_moments : int = 2: Number of moments to use for regularization. Only used if regularization > 0. Typical values: 2 (mean and variance) or 3 (mean, variance, skewness).
fixed : np.ndarray or list of tuples = None: Pin selected parameters at known constants so SVGD only optimises the learnable dimensions. Equivalent to a point-mass prior at value on the pinned slots, with those positions removed from the kernel and the gradient computation. Two accepted forms: (a) Index/value tuples (recommended) — [(idx, value), ...]. Each tuple pins parameter idx at value. Values may be any positive number (not just 1.0):: fixed=[(1, 0.01)] # theta[1] = 0.01, rest learned fixed=[(0, 2.5), (2, 0.1)] # theta[0]=2.5, theta[2]=0.1 (b) Binary mask (legacy) — a 1D array of length theta_dim where 1 pins the slot at the value 1.0 and 0 leaves it learnable. Use form (a) whenever the fix value is not 1.0:: fixed=[0, 1] # theta[1] pinned at 1.0 When to use. Fix a parameter when its value is known from prior data, when you want to test sensitivity to a single dimension while holding the rest at MLE, or when a slot is structurally unidentifiable from the observed data alone (e.g. a global rate scale that the model absorbs into another parameter). Interaction with prior. If prior is a per-slot list, the entries at fixed indices must be None; mismatches raise at construction time. A scalar prior callable is fine and is auto-masked at the fixed slots. Combination with tied. Tying across epochs (tied=[...]) is a daisy-chain-only feature owned by :meth:Graph.svgd; it is not accepted by this direct SVGD(model=...) constructor. When using Graph.svgd with both kwargs, fixed and tied are compatible as long as no (local_idx, epoch) slot is claimed by both (rule R20). To pin a parameter at the same value in every epoch use fixed=[(local_idx, value)]; to learn a single value shared across a subset of epochs, use tied. Performance. SVGD operates in reduced parameter space (learnable dims only). Kernel bandwidth is computed only over varying dimensions, improving convergence.
tied : not accepted on this path: Tying a parameter slot across multiple daisy-chain epochs (tied=[(local_idx, [epoch_a, epoch_b, ...]), ...]) is a pre-model-construction concern owned by :meth:Graph.svgd — it requires building the daisy-chain JSP model with an _apply_tying scatter wired into the forward pass, which cannot be retrofitted onto an opaque model callable handed to SVGD(model=...). The direct SVGD constructor therefore does not accept a tied= kwarg, and passing one raises TypeError. When to use. If you need to share a single learned value across two or more epochs of a time-inhomogeneous (daisy-chain) model — e.g. a baseline rate that is biologically constant across a subset of epochs while other parameters change — switch to the high-level interface:: graph.svgd( observed_data=…, epoch_starts=[0.0, t1, t2, …], tied=[(local_idx, [epoch_master, epoch_slave, …])], fixed=[…], # optional, see fixed= rules ) See :meth:phasic.Graph.svgd for the full tied semantics: master/slave selection, the rule R16 requirement that epoch_starts is also set, the rule R20 non-overlap with fixed, and the per-epoch local-index convention shared with fixed. Single-epoch alternative. To pin a parameter at a known constant (rather than learn a shared value), use fixed= instead — that is accepted on this path. tied only makes sense when the shared value should be learned, which requires the daisy-chain machinery.
optimizer : (Adam, SGDMomentum, RMSprop, Adagrad, OptaxOptimizer) = None: Optimizer for adaptive per-parameter learning rates. Default is None. Options include: - Adam: (default)Standard Adam optimizer - SGDMomentum: SGD with momentum - RMSprop: RMSprop optimizer - Adagrad: Adagrad optimizer When an optimizer is provided, the learning_rate parameter is ignored (the optimizer has its own learning rate). Example: >>> from phasic import SVGD, Adam >>> # Default: uses Adam >>> svgd = SVGD(model, data, theta_dim=2) >>> # Explicit optimizer >>> svgd = SVGD(model, data, theta_dim=2, optimizer=Adam(learning_rate=0.01))
exposure : float, array-like, or None = None: Per-observation exposure :math:\alpha_i — a known multiplicative scaling on a rate-typed component of :math:\boldsymbol{\theta}. For observation i the model is evaluated at :math:\boldsymbol{\theta}^{(i)} where :math:\theta^{(i)}_j = \theta_j for :math:j \neq k and :math:\theta^{(i)}_k = \theta_k \cdot \alpha_i, with :math:k = exposure_param_index. The exposed rate parameter and :math:\alpha_i jointly determine each observation’s expected event count (or hazard, or PMF, depending on the model). This is the GLM “exposure” / “offset” construct: it linearises the relationship between a rate parameter and an observation-specific outcome that scales with a known quantity. Concrete instances: - Coalescent-with-mutation: :math:\alpha_i = segment length :math:L_i in bases; :math:\theta_k is the per-base mutation rate. - Survival / failure-time: :math:\alpha_i = time-at-risk for unit :math:i; :math:\theta_k is the hazard rate. - Spatial Poisson: :math:\alpha_i = area or volume of region :math:i; :math:\theta_k is the intensity per unit area. - None (default): no exposure correction; existing behaviour. - scalar: same :math:\alpha applied to every observation. - 1D array of length n_observations: per-observation :math:\alpha. For dense 2D observed_data of shape (n_observations, n_features) the same :math:\alpha_i is shared across all features of observation :math:i. Requires exposure_param_index. Rejected for SparseObservations (raises NotImplementedError).
exposure_param_index : int or None = None: Index :math:k of the rate-typed parameter in :math:\boldsymbol{\theta} that exposure scales. Required when exposure is set. Must be in [0, param_length). Under daisy-chain (epoch_starts=[…]), :math:\boldsymbol{\theta} has flat layout (n_epochs * param_length,); exposure_param_index remains the local per-epoch index and is broadcast across every epoch internally.

Attributes

particles : array: Final posterior samples (n_particles, theta_dim)
theta_mean : array: Posterior mean estimate
theta_std : array: Posterior standard deviation
history : list of arrays, optional: Particle evolution over iterations (if fit was called with return_history=True)
is_fitted : bool: Whether fit() has been called

Examples

>>> # Basic usage with auto-configuration
>>> svgd = SVGD(model, observed_data, theta_dim=1)
>>> svgd.fit()
>>>
>>> # Explicit single-device configuration
>>> svgd = SVGD(model, observed_data, theta_dim=1, jit=True, parallel='vmap')
>>> svgd.fit()
>>>
>>> # Multi-device parallelization
>>> svgd = SVGD(model, observed_data, theta_dim=1,
...             jit=True, parallel='pmap', n_devices=8)
>>> svgd.fit()
>>>
>>> # No JIT (for debugging)
>>> svgd = SVGD(model, observed_data, theta_dim=1, jit=False, parallel='none')
>>> svgd.fit()
>>>
>>> # Multi-node SLURM (explicit distributed initialization)
>>> from phasic import initialize_distributed
>>> dist = initialize_distributed()  # Auto-detects SLURM environment
>>> svgd = SVGD(model, observed_data, theta_dim=1,
...             jit=True, parallel='pmap', n_devices=dist.num_processes)
>>> svgd.fit()
>>>
>>> # Using step size schedules to prevent divergence
>>> from phasic import ExpStepSize
>>> schedule = ExpStepSize(first_step=0.1, last_step=0.01, tau=500.0)
>>> svgd = SVGD(model, observed_data, theta_dim=1, learning_rate=schedule)
>>> svgd.fit()
>>>
>>> # Using adaptive step size based on particle spread
>>> from phasic import AdaptiveStepSize
>>> schedule = AdaptiveStepSize(base_step=0.01, kl_target=0.1, adjust_rate=0.1)
>>> svgd = SVGD(model, observed_data, theta_dim=1, learning_rate=schedule)
>>> svgd.fit()
>>>
>>> # Using regularization schedules for moment matching
>>> from phasic import ExpRegularization
>>> reg_schedule = ExpRegularization(first_reg=5.0, last_reg=0.1, tau=500.0)
>>> svgd = SVGD(model, observed_data, theta_dim=1,
...             regularization=reg_schedule, nr_moments=2)
>>> svgd.fit()  # Starts with strong regularization, gradually reduces
>>>
>>> # Using CDF-based regularization schedule (bidirectional)
>>>
>>> # Constant regularization (no schedule)
>>> svgd = SVGD(model, observed_data, theta_dim=1, regularization=1.0, nr_moments=2)
>>> svgd.fit()
>>>
>>> # Using custom bandwidth schedule
>>> from phasic import LocalAdaptiveBandwidth
>>> bandwidth = LocalAdaptiveBandwidth(alpha=0.9, k_frac=0.1)
>>> svgd = SVGD(model, observed_data, theta_dim=1, kernel=bandwidth)
>>> svgd.fit()
>>>
>>> # Access results
>>> print(svgd.theta_mean)
>>> print(svgd.theta_std)
>>>
>>> # Generate diagnostic plots
>>> svgd.plot_posterior()
>>> svgd.plot_trace()

Methods

Name	Description
analyze_trace	Analyze SVGD convergence and suggest parameter improvements.
animate	Create an animation showing the evolution of a single parameter’s distribution.
animate_pairwise	Create an animated pairwise scatter plot showing SVGD particle evolution.
animate_parameter_pairs	Animate multiple parameter pairs simultaneously.
check_convergence	Monitor convergence of SVGD by tracking statistics for n-dimensional parameters.
estimate_hdr	Estimate the Highest Density Region (HDR) from particles.
get_results	Get inference results as a dictionary.
log_likelihood	Maximum log-likelihood at the MAP (or at a supplied theta).
map_estimate_from_particles	Find the MAP estimate from a set of particles by finding the particle
map_estimate_with_optimization	Refine MAP estimate by starting from the best particle and performing
optimize	Run SVGD inference with optional moment-based regularization.
plot_ci	Plot mean parameter trajectory with credible interval ribbon.
plot_convergence	Plot convergence diagnostics showing mean and std over iterations.
plot_hdr	Plot 2D highest-density region with optional marginal HPD bands.
plot_pairwise	Plot pairwise scatter plots for all parameter pairs.
plot_posterior	Plot posterior distributions for each parameter.
plot_trace	Plot trace plots showing particle evolution over iterations.
summary	Print a summary of the inference results.

analyze_trace

phasic.svgd.SVGD.analyze_trace(burnin=None, verbose=True, return_dict=False)

Analyze SVGD convergence and suggest parameter improvements.

Computes convergence diagnostics, detects issues, and recommends parameter updates for better performance.

Parameters

burnin : int = None: Number of initial iterations to discard as burn-in. If None, auto-detects using convergence detection.
verbose : bool = True: Print detailed diagnostic report
return_dict : bool = False: Return full diagnostics dictionary

Returns

: dict or None: If return_dict=True, returns diagnostics dictionary with: - ‘converged’: bool - Whether SVGD converged - ‘convergence_point’: int or None - Iteration where converged - ‘diversity’: dict - Particle diversity metrics - ‘variance_collapse’: dict - Variance collapse diagnostics - ‘suggestions’: dict - Recommended parameter updates - ‘issues’: list - Detected problems

Raises

: RuntimeError: If fit() not called or history not available

Examples

>>> svgd = SVGD(model, data, theta_dim=1, n_iterations=700)
>>> svgd.fit(return_history=True)
>>> svgd.analyze_trace()  # Prints diagnostic report

>>> # Get full diagnostics
>>> diag = svgd.analyze_trace(return_dict=True, verbose=False)
>>> if not diag['converged']:
>>>     print("Need more iterations!")

animate

phasic.svgd.SVGD.animate(
    param_idx=0,
    true_theta=None,
    param_name=None,
    figsize=(8, 3),
    skip=0,
    thin=1,
    interval=100,
    duration=None,
    bins=30,
    show_particles=True,
    max_particles=20,
    save_as_gif=None,
    save_as_mp4=None,
    unconstrained=False,
)

Create an animation showing the evolution of a single parameter’s distribution.

This method creates a side-by-side animation with: - Left panel: Histogram of current parameter distribution - Right panel: Particle trajectories over time

Parameters

param_idx : int = 0: Index of the parameter to animate (0-indexed)
true_theta : np.ndarray = None: True parameter values. If provided, will overlay the true value for param_idx.
param_name : str = None: Name for the parameter (e.g., ‘jump rate’). If None, uses ‘θ_{param_idx}’.
figsize : tuple = (8, 3): Figure size (width, height)
skip : int = 0: Number of initial iterations to skip in the animation
thin : int, thin=1 = 1: Interval of interations to plot/annimate
interval : int = 100: Delay between frames in milliseconds
duration : int = None: Duration of the animation in seconds, overrides interval and thin if set
bins : int = 30: Number of histogram bins
show_particles : bool = True: If True, show individual particle trajectories in right panel
max_particles : int = 20: Maximum number of particle trajectories to show (for clarity)
save_as_gif : str = None: Path to save animation as GIF (requires pillow)
save_as_mp4 : str = None: Path to save animation as MP4 (requires ffmpeg)
unconstrained : bool = False: If False, show constrained (model-space) parameter values. If True, show unconstrained (optimization-space) values. Only relevant when using parameter transformations.

Returns

: IPython.display.HTML: Animation as HTML for Jupyter notebook display

Examples

>>> svgd = SVGD(model, data, theta_dim=3, n_iterations=70)
>>> svgd.fit(return_history=True)
>>> anim = svgd.animate(param_idx=0, true_theta=[2.0, 3.0, 2.0],
...                     param_name='jump rate')

animate_pairwise

phasic.svgd.SVGD.animate_pairwise(
    true_theta=None,
    param_names=None,
    figsize=None,
    skip=0,
    thin=1,
    interval=100,
    duration=None,
    save_as_gif=None,
    save_as_mp4=None,
    unconstrained=False,
)

Create an animated pairwise scatter plot showing SVGD particle evolution.

Parameters

true_theta : np.ndarray = None: True parameter values (if known) to overlay as red ‘x’ markers
param_names : list of str = None: Names for each parameter dimension (e.g., [‘jump’, ‘flood_left’, ‘flood_right’])
figsize : tuple = None: Figure size (width, height). Auto-sized based on parameter dimension if None.
skip : int = 0: Number of initial iterations to skip in the animation
thin : int, thin=1 = 1: Interval of interations to plot/annimate
interval : int = 100: Delay between frames in milliseconds
duration : int = None: Duration of the animation in seconds, overrides interval and thin if set
save_as_gif : str = None: Path to save animation as GIF (requires pillow)
save_as_mp4 : str = None: Path to save animation as MP4 (requires ffmpeg)
unconstrained : bool = False: If False, show constrained (model-space) parameter values. If True, show unconstrained (optimization-space) values. Only relevant when using parameter transformations.

Returns

: IPython.display.HTML: Animation as HTML for Jupyter notebook display

Raises

: RuntimeError: If fit() was not called with return_history=True
: ImportError: If matplotlib or required animation backend is not installed

Examples

>>> svgd = SVGD(model, data, theta_dim=3, n_iterations=70)
>>> svgd.fit(return_history=True)
>>> anim = svgd.animate_pairwise(
...     true_theta=[2.0, 3.0, 2.0],
...     param_names=['jump', 'flood_left', 'flood_right'],
...     save_as_gif='svgd_evolution.gif'
... )

animate_parameter_pairs

phasic.svgd.SVGD.animate_parameter_pairs(
    param_pairs=None,
    true_params=None,
    figsize=(15, 5),
    hide_fixed=True,
    save_as_gif=None,
)

Animate multiple parameter pairs simultaneously.

Parameters

param_pairs : list of tuple[int, int] = None: Parameter pairs to show, e.g. [(0, 1), (2, 3)]. Defaults to consecutive pairs.
true_params : np.ndarray = None: True parameter values for comparison.
figsize : tuple = (15, 5): Figure size (width, height).
hide_fixed : bool = True: If True, hide fixed parameters in the trace plot.
save_as_gif : str = None: Path to save animation as GIF.

check_convergence

phasic.svgd.SVGD.check_convergence(every=1, text=None, param_indices=None)

Monitor convergence of SVGD by tracking statistics for n-dimensional parameters.

estimate_hdr

phasic.svgd.SVGD.estimate_hdr(alpha=0.95)

Estimate the Highest Density Region (HDR) from particles.

Parameters

alpha : float = 0.95: Coverage probability (e.g., 0.95 for 95% HDR).

Returns

hdr_particles : array: Particles that are within the HDR.
threshold : float: The log probability threshold that defines the HDR.

get_results

phasic.svgd.SVGD.get_results()

Get inference results as a dictionary.

Returns

: dict: Dictionary containing: - ‘particles’: Final posterior samples (in constrained space if using transformation) - ‘theta_mean’: Posterior mean (in constrained space if using transformation) - ‘theta_std’: Posterior standard deviation (in constrained space if using transformation) - ‘history’: Particle evolution (if available, in constrained space if using transformation) - ‘particles_unconstrained’: Particles in unconstrained space (only if using transformation)

log_likelihood

phasic.svgd.SVGD.log_likelihood(theta=None, *, refine=False)

Maximum log-likelihood at the MAP (or at a supplied theta).

Returns the pure log-likelihood Σ log p(x_i | θ̂) with no prior and no moment-regularization penalty. This is the quantity consumed by AIC, BIC, and likelihood-ratio tests.

For partial-coverage reward models (Graph.svgd attaches a _zero_inflated_p_fn to the model when rewards don’t cover every absorbing trajectory), the returned value also includes the explicit zero-inflation point-mass term Σ_j n_zero_j · log(1 − p_j(θ)), matching the likelihood SVGD optimised against. Prior to v0.28.7 this term was silently dropped here, which biased AIC/BIC for partial-coverage fits.

Parameters

theta : array - like or None = None: If None, evaluate at the MAP particle (the one maximizing the log-posterior, found via :meth:map_estimate_from_particles). If refine=True, the MAP is refined via gradient ascent on the log-posterior before LL is evaluated there. The returned value is always pure LL — the prior is used only to locate the MAP, not added to the result. If an array is provided, it is interpreted in CONSTRAINED (model) space — i.e. the same space as theta_mean. The model is evaluated directly without re-applying param_transform. Must have shape (theta_dim,).
refine : bool = False: When theta is None, controls whether the MAP is refined via gradient ascent on the log-posterior (70 steps, step size 0.01) before LL is evaluated. Ignored otherwise.

Returns

: float: Σ log p(x_i | θ̂) over all valid observations.

Raises

: RuntimeError: If the SVGD instance has not been fitted yet.
: ValueError: If theta has wrong shape, or if the observed data has zero valid entries.

Examples

Maximum log-likelihood at the MAP, for AIC/BIC:

>>> svgd = graph.svgd(observed_data, theta_dim=2).optimize()
>>> ll_map = svgd.log_likelihood()
>>> ll_refined = svgd.log_likelihood(refine=True)  # closer to MLE

Evaluate LL at a specific θ (constrained space):

>>> import jax.numpy as jnp
>>> ll_at = svgd.log_likelihood(theta=jnp.array([1.5, 2.0]))

map_estimate_from_particles

phasic.svgd.SVGD.map_estimate_from_particles(unconstrained=False)

Find the MAP estimate from a set of particles by finding the particle with the highest log probability.

Parameters

unconstrained : bool = False: If False, return constrained (model-space) parameter values. If True, return unconstrained (optimization-space) values. Only relevant when using parameter transformations.

Returns

: tuple[list, float]: Tuple of (parameter values as list, log probability).

map_estimate_with_optimization

phasic.svgd.SVGD.map_estimate_with_optimization(
    n_steps=70,
    step_size=0.01,
    unconstrained=False,
)

Refine MAP estimate by starting from the best particle and performing gradient ascent on the log probability.

Parameters

n_steps : int = 70: Number of optimization steps.
step_size : float = 0.01: Step size for gradient ascent.
unconstrained : bool = False: If False, return constrained (model-space) parameter values. If True, return unconstrained (optimization-space) values. Only relevant when using parameter transformations.

Returns

: tuple[list, float]: Tuple of (refined parameter values as list, log probability).

optimize

phasic.svgd.SVGD.optimize(rewards=None, return_history=True)

Run SVGD inference with optional moment-based regularization.

Regularization settings are configured at SVGD initialization via regularization and nr_moments parameters.

Parameters

rewards : np.ndarray = None: Reward vector/matrix for multivariate phase-type distributions. - 1D array (n_vertices,): Single reward vector for univariate distribution - 2D array (n_features, n_vertices): Reward matrix for multivariate distribution where each row rewards[j, :] defines the reward vector for feature j Must be provided if model requires rewards (multivariate distributions). Each reward serves as multiplier of vertex value in trace.
return_history : bool = True: If True, store particle positions throughout optimization

Returns

: SVGD: Returns self for method chaining.

Raises

: NotImplementedError: If rewards parameter is provided (not yet implemented)

Examples

>>> # Standard SVGD (no regularization)
>>> model = Graph.pmf_and_moments_from_graph(graph)
>>> svgd = SVGD(model, observed_data, theta_dim=1, regularization=0.0)
>>> svgd.optimize()

>>> # SVGD with moment regularization
>>> model = Graph.pmf_and_moments_from_graph(graph, nr_moments=2)
>>> svgd = SVGD(model, observed_data, theta_dim=1, regularization=1.0, nr_moments=2)
>>> svgd.optimize()

>>> # With custom moments and strong regularization
>>> svgd = SVGD(model, observed_data, theta_dim=1, regularization=5.0, nr_moments=3)
>>> svgd.optimize()

Notes

Supports all JAX transformations: jit, grad, vmap, pmap
Supports multi-core parallelization via parallel=‘vmap’/‘pmap’
Supports multi-machine distribution via initialize_distributed()
Gradient compilation is cached (both memory and disk) for performance
All functionality from fit() and fit_regularized() is preserved

plot_ci

phasic.svgd.SVGD.plot_ci(
    figsize=(7, 3),
    save_path=None,
    skip=0,
    unconstrained=False,
    true_theta=None,
    ci=0.95,
    alpha=0.2,
    target=None,
    median=False,
    return_fig=False,
    ci_method='hpd',
    hide_fixed=True,
)

Plot mean parameter trajectory with credible interval ribbon.

Shows the posterior mean over iterations with a shaded region representing the specified credible interval (default 95%).

Parameters

figsize : tuple = (7, 3): Figure size (width, height)
save_path : str = None: Path to save the plot
skip : int = 0: Number of initial iterations to skip
unconstrained : bool = False: If False, show constrained (model-space) parameter values. If True, show unconstrained (optimization-space) values.
true_theta : np.ndarray = None: True parameter values to overlay as horizontal lines
ci : float = 0.95: Credible interval level (e.g., 0.95 for 95% CI)
alpha : float = 0.2: Transparency of the CI ribbon
target : np.ndarray = None: Target parameter value to overlay as horizontal line
median : bool = False: If True, plot the median trajectory as a dashed line
return_fig : bool = True: If True, return (fig, ax). If False, call plt.show() instead.
ci_method : str = 'hpd': Method for credible intervals. - 'hpd': Highest Posterior Density interval — the shortest interval containing ci fraction of posterior samples at each iteration. Computed by sorting the particles and sliding a window of ceil(n * ci) samples to find the narrowest span. - 'percentile': Equal-tailed percentile interval using symmetric quantiles.
hide_fixed : bool = True: If True, hide fixed parameters in the trace plot.

Returns

: tuple[matplotlib.figure.Figure, matplotlib.axes.Axes]: Matplotlib figure and axes (only if return_fig=True).

plot_convergence

phasic.svgd.SVGD.plot_convergence(
    figsize=(7, 3),
    save_path=None,
    skip=0,
    unconstrained=False,
    hide_fixed=True,
    return_fig=False,
)

Plot convergence diagnostics showing mean and std over iterations.

Requires fit() to have been called with return_history=True.

Parameters

figsize : tuple = (7, 4): Figure size (width, height)
save_path : str = None: Path to save the plot
skip : int = 0: Number of initial iterations to skip. Defaults to 0.
unconstrained : bool = False: If False, show constrained (model-space) parameter values. If True, show unconstrained (optimization-space) values. Only relevant when using parameter transformations.
hide_fixed : bool = True: If True, hide fixed parameters in the trace plot.
return_fig : bool = True: If True, return (fig, axes). If False, call plt.show() instead.

Returns

: tuple[matplotlib.figure.Figure, numpy.ndarray]: Matplotlib figure and axes array (only if return_fig=True).

plot_hdr

phasic.svgd.SVGD.plot_hdr(
    alphas=[0.95, 0.5],
    idx=None,
    figsize=(5, 4),
    hexgrid=True,
    trim=True,
    n=15,
    margin=0.1,
    xlim=None,
    ylim=None,
    palette='viridis',
    pad=2,
    unconstrained=False,
    return_fig=False,
    hide_fixed=True,
    show_hpd=False,
    hpd_alpha=0.95,
)

Plot 2D highest-density region with optional marginal HPD bands.

Displays a hex-grid log-likelihood heatmap and KDE-based HDR contours for two selected parameter dimensions.

Parameters

alphas : list of float = [0.95, 0.5]: HDR contour levels (fraction of mass enclosed).
idx : list of int = None: Indices of the two parameter dimensions to plot. Defaults to first two non-fixed parameters.
figsize : tuple = (5, 4): Figure size (width, height).
hexgrid : bool = True: Whether to show the log-likelihood hex-grid heatmap.
trim : bool = True: Clip axes to the hex-grid extent.
n : int = 15: Approximate number of hexagons along the shorter axis.
margin : float = 0.1: Fractional margin around particle range.
xlim : tuple = None: Manual axis limits.
ylim : tuple = None: Manual axis limits.
palette : str = 'viridis': Colormap for the hex-grid heatmap.
pad : int = 2: Extra hex rows/columns beyond the data range.
unconstrained : bool = False: Show unconstrained (optimization-space) values.
return_fig : bool = False: If True, return the axes object instead of calling plt.show().
show_hpd : bool = False: If True, overlay marginal HPD intervals as translucent vertical and horizontal bands. The HPD interval is the shortest contiguous interval containing hpd_alpha fraction of the marginal particle samples (computed by sorting the samples and sliding a window of ceil(n * hpd_alpha) elements to find the narrowest span).
hpd_alpha : float = 0.95: Credible level for the marginal HPD bands (only used when show_hpd=True).
hide_fixed : bool = True: If True, hide fixed parameters in the trace plot.

plot_pairwise

phasic.svgd.SVGD.plot_pairwise(
    true_theta=None,
    param_names=None,
    figsize=None,
    save_path=None,
    hide_fixed=True,
    unconstrained=False,
    return_fig=False,
)

Plot pairwise scatter plots for all parameter pairs.

Parameters

true_theta : np.ndarray = None: True parameter values (if known) to overlay on plot
param_names : list of str = None: Names for each parameter dimension
figsize : tuple = None: Figure size (width, height)
save_path : str = None: Path to save the plot
unconstrained : bool = False: If False, show constrained (model-space) parameter values. If True, show unconstrained (optimization-space) values. Only relevant when using parameter transformations.
hide_fixed : bool = True: If True, hide fixed parameters in the trace plot.
return_fig : bool = True: If True, return (fig, axes). If False, call plt.show() instead.

Returns

: tuple[matplotlib.figure.Figure, numpy.ndarray]: Matplotlib figure and axes array (only if return_fig=True).

plot_posterior

phasic.svgd.SVGD.plot_posterior(
    true_theta=None,
    param_names=None,
    bins=20,
    figsize=None,
    save_path=None,
    unconstrained=False,
    return_fig=False,
    ci_method='hpd',
    ci_level=0.95,
)

Plot posterior distributions for each parameter.

Parameters

true_theta : np.ndarray = None: True parameter values (if known) to overlay on plot
param_names : list of str = None: Names for each parameter dimension
bins : int = 20: Number of histogram bins
figsize : tuple = None: Figure size (width, height)
save_path : str = None: Path to save the plot
unconstrained : bool = False: If False, show constrained (model-space) parameter values. If True, show unconstrained (optimization-space) values. Only relevant when using parameter transformations.
return_fig : bool = True: If True, return (fig, axes). If False, call plt.show() instead.
ci_method : str = 'hpd': Method for credible intervals shown as a shaded region on each histogram. - 'hpd': Highest Posterior Density interval — the shortest interval containing ci_level fraction of posterior samples. Computed by sorting the particles and sliding a window of ceil(n * ci_level) samples to find the narrowest span. Better centred on the mode for skewed posteriors. - 'percentile': Equal-tailed percentile interval using symmetric quantiles.
ci_level : float = 0.95: Credible level for the shaded interval (e.g. 0.95 for 95%).

Returns

: tuple[matplotlib.figure.Figure, numpy.ndarray]: Matplotlib figure and axes array (only if return_fig=True).

plot_trace

phasic.svgd.SVGD.plot_trace(
    param_names=None,
    figsize=None,
    skip=0,
    max_particles=None,
    save_path=None,
    unconstrained=False,
    hide_fixed=True,
    return_fig=False,
)

Plot trace plots showing particle evolution over iterations.

Requires fit() to have been called with return_history=True.

Parameters

param_names : list of str = None: Names for each parameter dimension
figsize : tuple = None: Figure size (width, height)
skip : int = 0: Number of initial iterations to skip. Defaults to 0.
max_particles : int = None: Max number of particles to plot. Defaults to all particles.
save_path : str = None: Path to save the plot
unconstrained : bool = False: If False, show constrained (model-space) parameter values. If True, show unconstrained (optimization-space) values. Only relevant when using parameter transformations.
hide_fixed : bool = True: If True, hide fixed parameters in the trace plot.
return_fig : bool = True: If True, return (fig, axes). If False, call plt.show() instead.

Returns

: tuple[matplotlib.figure.Figure, numpy.ndarray]: Matplotlib figure and axes array (only if return_fig=True).

summary

phasic.svgd.SVGD.summary(ci_method='hpd', ci_level=0.95)

Print a summary of the inference results.

Parameters

ci_method : str = 'hpd': Method for credible intervals. - 'hpd': Highest Posterior Density interval — the shortest interval containing ci_level fraction of the posterior samples. Computed by sorting the particles and sliding a window of ceil(n * ci_level) samples to find the narrowest span. Better centred on the mode for skewed posteriors. - 'percentile': Equal-tailed percentile interval using the (1 - ci_level)/2 and (1 + ci_level)/2 quantiles.
ci_level : float = 0.95: Credible level (e.g. 0.95 for 95% intervals).

SVGD

When to use Graph.svgd vs SVGD direct

Parameters

Attributes

Examples

Methods

analyze_trace

Parameters

Returns

Raises

Examples

animate

Parameters

Returns

Examples

animate_pairwise

Parameters

Returns

Raises

Examples

animate_parameter_pairs

Parameters

check_convergence

estimate_hdr

Parameters

Returns

get_results

Returns

log_likelihood

Parameters

Returns

Raises

Examples

See Also

map_estimate_from_particles

Parameters

Returns

map_estimate_with_optimization

Parameters

Returns

optimize

Parameters

Returns

Raises

Examples

Notes

plot_ci

Parameters

Returns

plot_convergence

Parameters

Returns

plot_hdr

Parameters

plot_pairwise

Parameters

Returns

plot_posterior

Parameters

Returns

plot_trace

Parameters

Returns

summary

Parameters

When to use `Graph.svgd` vs `SVGD` direct