Graph

phasic.Graph(arg, ipv=None, cache_graph=False, **kwargs)

Attributes

Name	Description
cache_trace	Whether this graph caches the elimination trace.
hierarchical	Deprecated: use cache_trace instead.
trace_valid	Whether the cached trace is valid (not dirty).

Methods

Name	Description
absorbing_state_rewards	Get a reward vector that is 1 for states with edges to absorbing states.
accumulated_occupancy	Compute expected occupancy (visits or time) for each state.
accumulated_visiting_time	Compute expected time spent in each state (continuous only).
accumulated_visits	Compute expected number of visits to each state (discrete only).
as_matrices	Convert the graph to its matrix representation.
cdf	Compute cumulative distribution function.
clone	Create a deep copy of this graph.
compute_trace	Compute elimination trace with optional hierarchical caching.
copy	Returns a deep copy of the graph.
covariance	Compute covariance matrix for multivariate phase-type distributions.
discretize	Discretizes graph inplace and returns reward matrix for added auxiliary states.
distribution_context	Create a distribution context for efficient repeated sampling.
eliminate_to_dag	Perform symbolic graph elimination to create a reusable DAG structure.
expectation	Compute expected value (first moment) of the phase-type distribution.
expected_sojourn_time	Compute expected sojourn time (residence time) in each state.
expected_waiting_time	Compute expected waiting time until absorption.
extend	Extend the graph by continuing to visit unvisited vertices using a callback.
find_or_create_vertex	Find or create a vertex with the given state.
from_matrices	Construct a Graph from matrix representation.
from_serialized	Reconstruct Graph from serialize() output.
laplace_transform	Create a Laplace-transformed graph.
method_of_moments	Find parameter estimates by matching model moments to sample moments.
moments	Compute k-th moment of the phase-type distribution.
moments_from_graph	Convert a parameterized Graph to a JAX-compatible function that computes moments.
normalize	Normalize edge weights to make the graph a proper probability distribution.
pdf	Compute probability density/mass function using forward algorithm.
plot	Plots the graph using graphviz. See plot::plot_graph.py for more details.
pmf_and_moments_from_graph	Convert a parameterized Graph to a function that computes both PMF/PDF and moments.
pmf_and_moments_from_graph_multivariate	Create a multivariate phase-type model that handles 2D observations and rewards.
pmf_from_cpp	Load a phase-type model from a user’s C++ file and return a JAX-compatible function.
pmf_from_graph	Convert a Python-built Graph to a JAX-compatible function with full gradient support.
pmf_from_graph_joint_index	Create a JAX-compatible model for joint index distributions.
pmf_from_graph_parameterized	Convert a parameterized Python graph builder to a JAX-compatible function.
probability_matching	Find parameter estimates by matching model probabilities to empirical probabilities.
reward_transform	Apply reward transformation to create a new graph with modified rewards.
reward_transform_discrete	Apply reward transformation for discrete-time distributions.
sample	Generate random samples from the phase-type distribution.
serialize	Serialize graph to array representation for efficient computation.
state_probability	Compute probability of being in each state at a given time.
stop_probability	Compute probability of being in each state at a given time.
svgd	Run Stein Variational Gradient Descent (SVGD) inference for Bayesian parameter estimation.
update_weights	Update parameterized edge weights with given parameters.
variance	Compute variance of the phase-type distribution.

absorbing_state_rewards

phasic.Graph.absorbing_state_rewards()

Get a reward vector that is 1 for states with edges to absorbing states.

This reward vector is used for computing the Laplace transform value from a Laplace-transformed graph.

Returns

: np.ndarray: Reward vector of length n_vertices. Element i is 1.0 if vertex i has an edge to an absorbing state, 0.0 otherwise.

Examples

>>> graph = Graph(coalescent_callback)
>>> rewards = graph.absorbing_state_rewards()
>>> laplace_graph = graph.laplace_transform(0.5)
>>> L_0_5 = laplace_graph.expectation(rewards=rewards)

accumulated_occupancy

phasic.Graph.accumulated_occupancy(*args, **kwargs)

Compute expected occupancy (visits or time) for each state.

Parameters

*args : tuple = (): Additional positional arguments passed to C++ implementation.
****kwargs** : dict = {}: Additional keyword arguments passed to C++ implementation.

Returns

: np.ndarray: Expected visits (discrete) or time (continuous) in each state.

Notes

This is a convenience method that dispatches to either: - accumulated_visits() for discrete distributions - accumulated_visiting_time() for continuous distributions

accumulated_visiting_time

phasic.Graph.accumulated_visiting_time(*args, **kwargs)

Compute expected time spent in each state (continuous only).

Parameters

*args : tuple = (): Additional positional arguments passed to C++ implementation.
****kwargs** : dict = {}: Additional keyword arguments passed to C++ implementation.

Returns

: np.ndarray: Expected time spent in each state until absorption.

Notes

Only available for continuous-time phase-type distributions. For discrete distributions, use accumulated_visits() instead.

accumulated_visits

phasic.Graph.accumulated_visits(*args, **kwargs)

Compute expected number of visits to each state (discrete only).

Parameters

*args : tuple = (): Additional positional arguments passed to C++ implementation.
****kwargs** : dict = {}: Additional keyword arguments passed to C++ implementation.

Returns

: np.ndarray: Expected number of visits to each state until absorption.

Raises

: ValueError: If graph is not discrete.

Notes

Only available for discrete-time phase-type (DPH) distributions. The graph must be discretized via discretize() before calling this method.

as_matrices

phasic.Graph.as_matrices()

Convert the graph to its matrix representation.

Returns a NamedTuple containing the traditional phase-type distribution matrices and associated information.

Returns

: MatrixRepresentation: NamedTuple with the following attributes: - states: np.ndarray of shape (n_states, state_dim), dtype=int32 State vector for each vertex - sim: np.ndarray of shape (n_states, n_states), dtype=float64 Sub-intensity matrix - ipv: np.ndarray of shape (n_states,), dtype=float64 Initial probability vector - indices: np.ndarray of shape (n_states,), dtype=int32 1-based indices for vertices (for use with vertex_at())

Examples

>>> g = Graph(1)
>>> start = g.starting_vertex()
>>> v1 = g.find_or_create_vertex([1])
>>> v2 = g.find_or_create_vertex([2])
>>> start.add_edge(v1, 1.0)
>>> v1.add_edge(v2, 2.0)
>>> g.normalize()
>>>
>>> matrices = g.as_matrices()
>>> print(matrices.sim)  # Sub-intensity matrix (attribute access)
>>> print(matrices.ipv)  # Initial probability vector
>>> # Can also use index access like a tuple
>>> states, sim, ipv, indices = matrices

cdf

phasic.Graph.cdf(time, **kwargs)

Compute cumulative distribution function.

Parameters

time : float or ArrayLike: Time point(s) at which to evaluate the CDF.
****kwargs** : dict = {}: Additional keyword arguments passed to C++ implementation.

Returns

: float or np.ndarray: CDF value(s) P(T ≤ t) at the specified time point(s).

Notes

For continuous distributions: F(t) = P(T ≤ t) = 1 - α · exp(S·t) · 1 For discrete distributions: F(n) = P(N ≤ n) = sum of PMF up to n

clone

phasic.Graph.clone()

Create a deep copy of this graph.

The cloned graph preserves the cache_trace setting but starts with a fresh (invalid) trace cache.

Returns

: Graph: A new Graph instance with the same structure.

compute_trace

phasic.Graph.compute_trace(
    param_length=None,
    hierarchical=True,
    min_size=50,
    parallel='auto',
    verbose=False,
    force=False,
)

Compute elimination trace with optional hierarchical caching.

When Graph was created with cache_trace=True, the trace is cached on the instance for use by moments(), expectation(), etc. In this mode, the operation is NON-DESTRUCTIVE (graph is preserved via cloning).

When Graph was created with cache_trace=False (default), the operation is DESTRUCTIVE and will empty the graph during trace recording.

Parameters

param_length : int = None: Number of parameters (auto-detect if None)
hierarchical : bool = True: If True, use hierarchical SCC-based caching (recommended). If False, use direct trace recording without caching. Caching provides 10-100x speedup on repeated calls.
min_size : int = 50: Minimum vertices to subdivide (only used if hierarchical=True)
parallel : str = 'auto': Parallelization: ‘auto’, ‘vmap’, ‘pmap’, or ‘sequential’
verbose : bool = False: If True, show progress bars for major computation stages
force : bool = False: If True, recompute trace even if a cached trace exists and is valid. Only applicable when Graph was created with cache_trace=True.

Returns

: EliminationTrace: Elimination trace (from cache or computed)

Notes

Caching Mode (Graph created with cache_trace=True): - Trace is cached on the instance and reused by moments(), expectation(), etc. - Graph is cloned before trace recording, preserving the original - Subsequent calls return cached trace unless force=True or graph was modified

Non-Caching Mode (default): - The graph elimination algorithm is DESTRUCTIVE - vertices are eliminated - Use disk caching (hierarchical=True parameter) to avoid re-recording

Examples

>>> # Cache mode - non-destructive, cached on instance
>>> g = Graph(model, nr_samples=5, cache_trace=True)
>>> g.normalize()
>>> trace = g.compute_trace()  # Graph preserved, trace cached
>>> g.update_weights([1.0, 2.0])
>>> mean = g.expectation()  # Uses cached trace
>>>
>>> # Standard mode with disk caching
>>> g = Graph(model, nr_samples=5)
>>> trace = g.compute_trace()  # Graph emptied, trace returned

copy

phasic.Graph.copy()

Returns a deep copy of the graph.

Creates an independent copy with all vertices, edges, and metadata. Modifications to the copy will not affect the original graph.

Returns

: Graph: Deep copy of the graph

covariance

phasic.Graph.covariance(rewards1=[], rewards2=[], **kwargs)

Compute covariance matrix for multivariate phase-type distributions.

Parameters

rewards1 : list of float or ndarray = []: The first set of rewards, which should be applied to the phase-type distribution. Must have length equal to vertices_length().
rewards2 : list of float or ndarray = []: The second set of rewards, which should be applied to the phase-type distribution. Must have length equal to vertices_length().
****kwargs** : dict = {}: Additional keyword arguments passed to C++ implementation.

Returns

: np.ndarray: Covariance matrix for the multivariate distribution.

Notes

This method is for multivariate phase-type distributions with multiple marginals. For univariate distributions, use variance().

discretize

phasic.Graph.discretize(rate, skip_existing=False, **kwargs)

Discretizes graph inplace and returns reward matrix for added auxiliary states.

Parameters

rate : float | Callable: float or callable
skip_existing : bool = False: If True, skip vertices that already have auxiliary vertices, by default False

Returns

: NDArray[np.int64]: Reward matrix for added auxiliary states

distribution_context

phasic.Graph.distribution_context(*args, **kwargs)

Create a distribution context for efficient repeated sampling.

Parameters

*args : tuple = (): Additional positional arguments passed to C++ implementation.
****kwargs** : dict = {}: Additional keyword arguments passed to C++ implementation.

Returns

: object: Distribution context object that can be used for efficient sampling.

Notes

The distribution context precomputes data structures needed for sampling, making repeated sample() calls much faster than sampling directly from the graph.

eliminate_to_dag

phasic.Graph.eliminate_to_dag()

Perform symbolic graph elimination to create a reusable DAG structure.

This method performs the O(n³) graph elimination algorithm ONCE and returns a symbolic DAG where edges contain expression trees instead of concrete values. The DAG can then be instantiated with different parameters in O(n) time each.

This is the key optimization for SVGD and other inference methods that require evaluating the same graph structure with many different parameter vectors.

Returns

: SymbolicDAG: Symbolic DAG that can be instantiated with parameters

Raises

: RuntimeError: If the graph is not parameterized or elimination fails

Examples

>>> # Create parameterized graph
>>> g = Graph(1)
>>> v_a = g.create_vertex([0])
>>> v_b = g.create_vertex([1])
>>> v_c = g.create_vertex([2])
>>> v_a.add_edge_parameterized(v_b, 0.0, [1.0, 0.0, 0.0])
>>> v_b.add_edge_parameterized(v_c, 0.0, [0.0, 1.0, 0.0])

>>> # Eliminate to symbolic DAG (once)
>>> dag = g.eliminate_to_dag()
>>> print(dag)  # SymbolicDAG(vertices=3, params=3, acyclic=True)

>>> # Fast instantiation for SVGD (100-1000× faster!)
>>> for theta in particle_swarm:
...     g_concrete = dag.instantiate(theta)
...     log_prob = -g_concrete.expectation()  # Fast!

Performance

Elimination: O(n³) - performed once
Instantiation: O(n) - performed per particle
Expected speedup for SVGD: 100-1000×

expectation

phasic.Graph.expectation(rewards=[], **kwargs)

Compute expected value (first moment) of the phase-type distribution.

Parameters

rewards : ArrayLike = []: Reward vector for reward-transformed expectation. If not provided, computes E[T] where T is time until absorption.
****kwargs** : dict = {}: Additional keyword arguments passed to C++ implementation.

Returns

: float: Expected value E[T] or reward-transformed expectation.

Notes

For parameterized graphs, this method uses trace-based computation which requires O(n) memory instead of O(n²) for the matrix-based approach. Set cache_trace=True when creating the graph to cache the trace for faster repeated evaluations.

This is equivalent to moments(1, rewards).

expected_sojourn_time

phasic.Graph.expected_sojourn_time(*args, **kwargs)

Compute expected sojourn time (residence time) in each state.

Parameters

*args : tuple = (): Positional arguments passed to C++ implementation.
****kwargs** : dict = {}: Keyword arguments passed to C++ implementation.

Returns

: ArrayLike: Expected sojourn time for each vertex.

Notes

This method wraps the C++ implementation of expected_sojourn_time. Returns the expected accumulated reward for the starting vertex propagated through all paths. This is what expectation() uses internally. The difference from expected_residence_time is subtle - expected_residence_time decomposes this into per-vertex contributions.

Available for both continuous and discrete phase-type distributions.

expected_waiting_time

phasic.Graph.expected_waiting_time(*args, **kwargs)

Compute expected waiting time until absorption.

Parameters

*args : tuple = (): Positional arguments passed to C++ implementation.
****kwargs** : dict = {}: Keyword arguments passed to C++ implementation.

Returns

: float: Expected waiting time until absorption.

Notes

This method wraps the C++ implementation of expected_waiting_time. Only available for continuous-time phase-type distributions.

extend

phasic.Graph.extend(callback=None, **kwargs)

Extend the graph by continuing to visit unvisited vertices using a callback.

After manually adding vertices to the graph (e.g., via find_or_create_vertex), this method continues the callback-based graph building process, visiting all newly added unvisited vertices.

Parameters

callback : callable = None: Callback function to use for extending the graph. If None, uses the callback from the original Graph construction. The callback should accept a state array and return: - For non-parameterized: list of (state, weight) tuples - For parameterized: list of (state, base_weight, edge_state) tuples
****kwargs** : Any = {}: Additional keyword arguments to pass to the callback function. If None, uses the kwargs from the original Graph construction.

Raises

: RuntimeError: If no callback is provided and the graph was not constructed with a callback.

Examples

>>> # Build initial graph
>>> graph = Graph(my_callback, nr_samples=5)
>>>
>>> # Manually add new vertex
>>> special_vertex = graph.find_or_create_vertex([100, 200])
>>> graph.starting_vertex().add_edge_parameterized(special_vertex, 0.0, [1.5])
>>>
>>> # Continue with callback to explore new vertices
>>> graph.extend()  # Uses original callback
>>>
>>> # Or use different callback
>>> graph.extend(callback=my_other_callback, param=value)

find_or_create_vertex

phasic.Graph.find_or_create_vertex(state)

Find or create a vertex with the given state.

This method wraps the C++ implementation to track trace invalidation.

from_matrices

phasic.Graph.from_matrices(ipv, sim, states=None)

Construct a Graph from matrix representation.

Parameters

ipv : np.ndarray: Initial probability vector, shape (n_states,)
sim : np.ndarray: Sub-intensity matrix, shape (n_states, n_states)
states : np.ndarray = None: State vectors, shape (n_states, state_dim), dtype=int32 If None, uses default states [0], [1], [2], …

Returns

: Graph: The reconstructed phase-type distribution graph

Examples

>>> ipv = np.array([0.6, 0.4])
>>> sim = np.array([[-2.0, 1.0], [0.0, -3.0]])
>>> g = Graph.from_matrices(ipv, sim)
>>> pdf = g.pdf(1.0)

from_serialized

phasic.Graph.from_serialized(data)

Reconstruct Graph from serialize() output.

This method enables distributed trace recording by allowing graphs to be serialized to JSON, sent across the network via JAX pmap, and reconstructed on worker processes.

Parameters

data : dict[str, Any]: Dictionary returned by Graph.serialize() containing: - states: np.ndarray of shape (n_vertices, state_length) - edges: np.ndarray of shape (n_edges, 3) with [from_idx, to_idx, weight] - start_edges: np.ndarray of shape (n_start_edges, 2) with [to_idx, weight] - param_edges: np.ndarray of shape (n_param_edges, 2+param_length) Format (v0.22.0+): [from_idx, to_idx, coeff1, coeff2, …] - start_param_edges: np.ndarray of shape (n_start_param_edges, 1+param_length) Format (v0.22.0+): [to_idx, coeff1, coeff2, …] - param_length: int - state_length: int - n_vertices: int

Returns

: Graph: Reconstructed graph with identical structure to the original

Raises

: ValueError: If data is missing required fields, has wrong shapes, or malformed arrays
: RuntimeError: If graph reconstruction fails (e.g., edges to non-existent vertices)

Examples

>>> import json
>>> import numpy as np
>>> from phasic import Graph
>>>
>>> # Create and serialize graph
>>> g = Graph(1)
>>> v0 = g.starting_vertex()
>>> v1 = g.find_or_create_vertex([1])
>>> v0.add_edge_parameterized(v1, 0.0, [2.0])
>>> serialized = g.serialize(theta_dim=1)
>>>
>>> # Convert to JSON (for network transmission)
>>> json_dict = {k: v.tolist() if isinstance(v, np.ndarray) else v
...              for k, v in serialized.items()}
>>> json_str = json.dumps(json_dict)
>>>
>>> # Reconstruct on worker (e.g., different machine)
>>> received_dict = json.loads(json_str)
>>> received_dict['states'] = np.array(received_dict['states'], dtype=np.int32)
>>> received_dict['edges'] = np.array(received_dict['edges'], dtype=np.float64)
>>> received_dict['start_edges'] = np.array(received_dict['start_edges'], dtype=np.float64)
>>> received_dict['param_edges'] = np.array(received_dict['param_edges'], dtype=np.float64)
>>> received_dict['start_param_edges'] = np.array(received_dict['start_param_edges'], dtype=np.float64)
>>> g_reconstructed = Graph.from_serialized(received_dict)
>>> assert g_reconstructed.vertices_length() == g.vertices_length()

laplace_transform

phasic.Graph.laplace_transform(theta)

Create a Laplace-transformed graph.

Returns a new graph where each transient state has an additional edge to the absorbing state with weight theta (or theta added to existing absorbing edge weight).

To compute the Laplace transform value L(theta) = E[exp(-theta * T)], call expectation() on the result with a reward vector that is 1 for states that had edges to absorbing in the original graph and 0 otherwise.

Parameters

theta : float: The Laplace transform parameter.

Returns

: Graph: A new graph representing the Laplace-transformed distribution.

Examples

>>> import numpy as np
>>> graph = Graph(coalescent_callback)
>>> # Get reward vector: 1 for states with absorbing edges, 0 otherwise
>>> rewards = graph.absorbing_state_rewards()
>>> laplace_graph = graph.laplace_transform(0.5)
>>> laplace_value = laplace_graph.expectation(rewards=rewards)

Notes

The Laplace transform of a phase-type distribution PH(alpha, S) is:

L(theta) = E[exp(-theta * T)] = alpha * (theta*I - S)^(-1) * s

where s is the exit rate vector indicating states with direct transitions to absorption.

For continuous distributions only. For discrete distributions, use the z-transform instead.

method_of_moments

phasic.Graph.method_of_moments(
    observed_data,
    nr_moments=None,
    rewards=None,
    fixed=None,
    theta_dim=None,
    theta_init=None,
    std_multiplier=2.0,
    discrete=None,
    verbose=True,
)

Find parameter estimates by matching model moments to sample moments.

Solves the nonlinear least-squares problem::

minimize ||moments_fn(theta) - sample_moments||^2
subject to  theta > 0

The returned MoMResult.prior list can be passed directly to :meth:Graph.svgd as the prior argument, providing data-informed priors centred on the method-of-moments estimates.

Parameters

observed_data : ArrayLike: Observed data. 1-D for univariate, 2-D (n_times, n_features) for multivariate models with 2-D reward vectors.
nr_moments : int = None: Number of moments to match per feature dimension. If None (default), automatically chosen to overdetermine the system: max(2 * n_free_params, 4) for univariate data, adjusted per feature for multivariate data. Still auto-increased if a user-specified value gives fewer equations than free parameters.
rewards : np.ndarray = None: Reward vectors. None for standard moments, 1-D for a single reward vector, 2-D (n_features, n_vertices) for multivariate.
fixed : list = None: List of (index, value) tuples pinning specific parameters.
theta_dim : int = None: Number of model parameters. Inferred from the graph when None.
theta_init : np.ndarray = None: Initial guess for the free parameters (excluding fixed ones). If None a coordinate-wise grid search is used.
std_multiplier : float = 2.0: Factor applied to the asymptotic standard error to obtain the prior standard deviation: prior_std = std_multiplier * se.
discrete : bool = None: True for discrete (PMF) models, False for continuous (PDF). If None, inferred from self.is_discrete.
verbose : bool = True: Print progress information.

Returns

: MoMResult: Dataclass with theta, std, prior, success, etc.

Examples

>>> g = Graph(...)  # parameterized graph
>>> data = np.random.exponential(0.5, size=200)
>>> mom = g.method_of_moments(data)
>>> print(mom.theta)            # parameter estimate
>>> svgd = g.svgd(data, prior=mom.prior)  # use as SVGD prior

moments

phasic.Graph.moments(power, rewards=[], discrete=False, **kwargs)

Compute k-th moment of the phase-type distribution.

Parameters

power : int: Moment order (1 for first moment, 2 for second moment, etc.).
rewards : ArrayLike = []: Reward vector for reward-transformed moments. If not provided, uses unit rewards (standard moments).
discrete : bool = False: If True, compute discrete-time moments (DPH distribution). Requires that the graph was discretized via discretize().
****kwargs** : dict = {}: Additional keyword arguments passed to C++ implementation.

Returns

: float: The k-th moment E[T^k] where T is the time until absorption.

Notes

If the graph is parameterized, this method uses trace-based computation for 5-10x speedup on repeated evaluations and O(n) memory usage. Otherwise falls back to direct C++ graph elimination.

For higher moments (k > 2), numerical stability may become an issue for complex distributions.

moments_from_graph

phasic.Graph.moments_from_graph(graph, nr_moments=2, use_ffi=False)

Convert a parameterized Graph to a JAX-compatible function that computes moments.

This method creates a function that computes the first nr_moments moments of the phase-type distribution: [E[T], E[T^2], …, E[T^nr_moments]].

Moments are computed using the existing C++ graph.moments(power) method for efficiency.

Parameters

graph : Graph: Parameterized graph built using the Python API with parameterized edges. Must have edges created with add_edge_parameterized().
nr_moments : int = 2: Number of moments to compute. For example: - 1: Returns [E[T]] (mean only) - 2: Returns [E[T], E[T^2]] (mean and second moment) - 3: Returns [E[T], E[T^2], E[T^3]]
use_ffi : bool = False: If True, uses Foreign Function Interface approach.

Returns

: callable: JAX-compatible function with signature: moments_fn(theta) -> jnp.array(nr_moments,) Returns array of moments: [E[T], E[T^2], …, E[T^k]]

Examples

>>> # Create parameterized coalescent model
>>> def coalescent(state, nr_samples=2):
...     if len(state) == 0:
...         return [(np.array([nr_samples]), 1.0, [1.0])]
...     if state[0] > 1:
...         n = state[0]
...         rate = n * (n - 1) / 2
...         return [(np.array([n-1]), 0.0, [rate])]
...     return []
>>>
>>> graph = Graph(coalescent, nr_samples=3)
>>> moments_fn = Graph.moments_from_graph(graph, nr_moments=2)
>>>
>>> # Compute moments for given theta
>>> theta = jnp.array([0.5])
>>> moments = moments_fn(theta)  # [E[T], E[T^2]]
>>> print(f"Mean: {moments[0]}, Second moment: {moments[1]}")
>>>
>>> # Variance can be computed as: Var[T] = E[T^2] - E[T]^2
>>> variance = moments[1] - moments[0]**2

Notes

Requires graph to have parameterized edges (created with parameterized=True)
Moments are raw moments, not central moments
For variance, compute: Var[T] = E[T^2] - E[T]^2
For standard deviation: std[T] = sqrt(Var[T])

normalize

phasic.Graph.normalize(*args, **kwargs)

Normalize edge weights to make the graph a proper probability distribution.

Parameters

*args : tuple = (): Additional positional arguments passed to C++ implementation.
****kwargs** : dict = {}: Additional keyword arguments passed to C++ implementation.

Returns

: float: Scaling factor applied to normalize the weights.

Notes

This method modifies the graph in-place and invalidates any cached trace.

For continuous distributions: Normalizes transition rates. For discrete distributions: Normalizes transition probabilities to sum to 1.

The normalization ensures that the initial probability vector and transition matrix satisfy the requirements for a valid phase-type distribution.

pdf

phasic.Graph.pdf(time, **kwargs)

Compute probability density/mass function using forward algorithm.

Parameters

time : float or ArrayLike: Time point(s) at which to evaluate the PDF/PMF.
granularity : int: Granularity for uniformization (default: auto-detected as 2*max_rate). Higher values improve accuracy but increase computation time.
****kwargs** : dict = {}: Additional keyword arguments passed to C++ implementation.

Returns

: float or np.ndarray: PDF/PMF value(s) at the specified time point(s).

Notes

This method uses the forward algorithm (Algorithm 4) via uniformization to compute the exact phase-type PDF/PMF, not an approximation.

For continuous distributions: f(t) = α · exp(S·t) · s* For discrete distributions: p(n) = probability of absorption at jump n

plot

phasic.Graph.plot(*args, **kwargs)

Plots the graph using graphviz. See plot::plot_graph.py for more details.

Returns

: Any: description

pmf_and_moments_from_graph

phasic.Graph.pmf_and_moments_from_graph(
    graph,
    nr_moments=2,
    discrete=False,
    use_ffi=False,
    theta_dim=None,
)

Convert a parameterized Graph to a function that computes both PMF/PDF and moments.

This is more efficient than calling pmf_from_graph() and moments_from_graph() separately because it builds the graph once and computes both quantities.

Parameters

graph : Graph: Parameterized graph built using the Python API with parameterized edges.
nr_moments : int = 2: Number of moments to compute
discrete : bool = False: If True, computes discrete PMF. If False, computes continuous PDF.
use_ffi : bool = False: If True, uses Foreign Function Interface approach.
theta_dim : int = None: Number of parameters for parameterized edges. If not provided, will be auto-detected by probing edge states. Providing this explicitly avoids potential issues with auto-detection reading garbage memory.

Returns

: callable: JAX-compatible function with signature: model(theta, times, rewards=None) -> (pmf_values, moments) Where: - theta: Parameter vector - times: Time points or jump counts - rewards: Optional reward vector (one per vertex). If None, computes standard moments. - pmf_values: jnp.array(len(times),) - PMF/PDF values at each time - moments: jnp.array(nr_moments,) - [E[T], E[T^2], …] or [E[R·T], E[R·T^2], …]

Examples

>>> # Create parameterized model
>>> graph = Graph(coalescent, nr_samples=3)
>>> model = Graph.pmf_and_moments_from_graph(graph, nr_moments=2)
>>>
>>> # Compute both PMF and moments
>>> theta = jnp.array([0.5])
>>> times = jnp.array([1.0, 2.0, 3.0])
>>> pmf_vals, moments = model(theta, times)
>>>
>>> print(f"PMF at times: {pmf_vals}")
>>> print(f"Moments: {moments}")  # [E[T], E[T^2]]
>>>
>>> # Compute reward-transformed moments
>>> rewards = jnp.array([1.0, 2.0, 0.5, 1.5])  # One per vertex
>>> pmf_vals, reward_moments = model(theta, times, rewards=rewards)
>>> print(f"Reward moments: {reward_moments}")  # [E[R·T], E[R·T^2]]
>>>
>>> # Use in SVGD with moment regularization
>>> svgd = SVGD(model, observed_pmf, theta_dim=1)
>>> svgd.fit_regularized(observed_times=data, nr_moments=2, regularization=1.0)

Notes

More efficient than separate calls to pmf_from_graph() and moments_from_graph()
Required for using moment-based regularization in SVGD.fit_regularized()
The moments are always computed from the same graph used for PMF/PDF

pmf_and_moments_from_graph_multivariate

phasic.Graph.pmf_and_moments_from_graph_multivariate(
    graph,
    nr_moments=2,
    discrete=False,
    use_ffi=False,
    theta_dim=None,
)

Create a multivariate phase-type model that handles 2D observations and rewards.

This wrapper enables computing joint likelihoods for multivariate phase-type distributions where each feature dimension has its own reward vector defining the marginal distribution.

Parameters

graph : Graph: Parameterized graph built using the Python API with parameterized edges.
nr_moments : int = 2: Number of moments to compute per feature dimension
discrete : bool = False: If True, computes discrete PMF. If False, computes continuous PDF.
use_ffi : bool = False: If True, uses Foreign Function Interface approach.
theta_dim : int = None: Number of parameters for parameterized edges.

Returns

: callable: JAX-compatible function with signature: model(theta, times, rewards=None) -> (pmf_values, moments) Where: - theta: Parameter vector (theta_dim,) - times: Time points - can be 1D (n_times,) or 2D (n_times, n_features) - rewards: Reward vectors - can be: * None: Standard moments * 1D (n_vertices,): Single reward vector (backward compatible) * 2D (n_vertices, n_features): Multivariate - one reward per feature - pmf_values: PMF/PDF values - shape matches input: * 1D rewards → 1D output (n_times,) * 2D rewards → 2D output (n_times, n_features) - moments: Moment values: * 1D rewards → 1D output (nr_moments,) * 2D rewards → 2D output (n_features, nr_moments)

Examples

>>> # Create parameterized model
>>> graph = Graph(coalescent, nr_samples=3)
>>> model = Graph.pmf_and_moments_from_graph_multivariate(graph, nr_moments=2)
>>>
>>> # 1D case (backward compatible)
>>> theta = jnp.array([0.5])
>>> times = jnp.array([1.0, 2.0, 3.0])
>>> rewards_1d = jnp.array([1.0, 2.0, 0.5, 1.5])  # (n_vertices,)
>>> pmf_vals, moments = model(theta, times, rewards_1d)
>>> print(pmf_vals.shape)  # (3,)
>>> print(moments.shape)   # (2,)
>>>
>>> # 2D case (multivariate)
>>> rewards_2d = jnp.array([[1.0, 2.0, 0.5, 1.5],   # Feature 1 reward vector
...                          [0.5, 1.0, 2.0, 0.8]])  # Feature 2 reward vector
...                                                  # (n_features, n_vertices)
>>> times_2d = jnp.array([[1.0, 1.5],
...                        [2.0, 2.5],
...                        [3.0, 3.5]])  # (n_times, n_features)
>>> pmf_vals, moments = model(theta, times_2d, rewards_2d)
>>> print(pmf_vals.shape)  # (3, 2)
>>> print(moments.shape)   # (2, 2)
>>>
>>> # Use in SVGD with 2D observations
>>> observed_data = jnp.array([[1.5, 2.1], [0.8, 1.2], [2.3, 3.1]])
>>> svgd = SVGD(model, observed_data, theta_dim=1, rewards=rewards_2d)
>>> results = svgd.optimize()

Notes

For 2D rewards, each feature dimension is computed independently using the corresponding row of the rewards matrix (rewards[j, :] for feature j)
Log-likelihood is computed as sum over all observation elements
Backward compatible: 1D rewards behave exactly as pmf_and_moments_from_graph()

pmf_from_cpp

phasic.Graph.pmf_from_cpp(cpp_file, discrete=False)

Load a phase-type model from a user’s C++ file and return a JAX-compatible function.

The C++ file should include ‘user_model.h’ and implement:

phasic::Graph build_model(const double* theta, int n_params) { // Build and return Graph instance }

For efficient repeated evaluations with the same parameters without JAX overhead, use load_cpp_builder() instead to get a builder function that creates Graph objects.

Parameters

cpp_file : str or pathlib.Path: Path to the user’s C++ file
discrete : bool = False: If True, uses discrete phase-type distribution (DPH) computation. If False, uses continuous phase-type distribution (PDF).

Raises

: ImportError: If JAX is not installed. Install with: pip install jax jaxlib
: FileNotFoundError: If the specified C++ file does not exist

Returns

: callable: JAX-compatible function (theta, times) -> pmf_values that supports JIT, grad, vmap, etc.

Examples

JAX-compatible approach (default - for SVGD, gradients, optimization)

>>> model = Graph.pmf_from_cpp("my_model.cpp")
>>> theta = jnp.array([1.0, 2.0])
>>> times = jnp.linspace(0, 10, 100)
>>> pmf = model(theta, times)
>>> gradient = jax.grad(lambda p: jnp.sum(model(p, times)))(theta)

Discrete phase-type distribution

>>> model = Graph.pmf_from_cpp("my_model.cpp", discrete=True)
>>> theta = jnp.array([1.0, 2.0])
>>> jumps = jnp.array([1, 2, 3, 4, 5])
>>> dph_pmf = model(theta, jumps)

For direct C++ access without JAX (faster for repeated evaluations):

>>> builder = load_cpp_builder("my_model.cpp")
>>> graph = builder(np.array([1.0, 2.0]))  # Build graph once
>>> pdf1 = graph.pdf(1.0)  # Use many times
>>> pdf2 = graph.pdf(2.0)  # No rebuild needed

pmf_from_graph

phasic.Graph.pmf_from_graph(
    graph,
    discrete=False,
    use_cache=True,
    theta_dim=None,
)

Convert a Python-built Graph to a JAX-compatible function with full gradient support.

This method automatically detects if the graph has parameterized edges (edges with state vectors) and generates optimized C++ code to enable full JAX transformations including gradients, vmap, and jit compilation.

For direct C++ access without JAX wrapping, use the Graph object’s methods directly: graph.pdf(time), graph.dph_pmf(jump), graph.moments(power), etc.

Raises

: ImportError: If JAX is not installed. Install with: pip install jax jaxlib

Parameters

graph : Graph: Graph built using the Python API. Can have regular edges or parameterized edges.
discrete : bool = False: If True, uses discrete phase-type distribution (DPH) computation. If False, uses continuous phase-type distribution (PDF).
use_cache : bool = True: If True, uses symbolic DAG cache to avoid re-computing expensive symbolic elimination for graphs with the same structure. Default: True Set to False to disable caching (useful for testing).

Returns

: callable: If graph has parameterized edges: JAX-compatible function (theta, times) -> pmf_values Supports JIT, grad, vmap, etc. If graph has no parameterized edges: JAX-compatible function (times) -> pmf_values Supports JIT (backward compatible signature)

Examples

Non-parameterized graph (regular edges only)

>>> g = Graph(1)
>>> start = g.starting_vertex()
>>> v0 = g.find_or_create_vertex([0])
>>> v1 = g.find_or_create_vertex([1])
>>> start.add_edge(v0, 1.0)
>>> v0.add_edge(v1, 2.0)  # fixed weight
>>>
>>> model = Graph.pmf_from_graph(g)
>>> times = jnp.linspace(0, 5, 50)
>>> pdf = model(times)  # No theta needed

Parameterized graph (with edge states for gradient support)

>>> g = Graph(1)
>>> start = g.starting_vertex()
>>> v0 = g.find_or_create_vertex([0])
>>> v1 = g.find_or_create_vertex([1])
>>> start.add_edge(v0, 1.0)
>>> v0.add_edge_parameterized(v1, 0.0, [2.0, 0.5])  # weight = 2.0*theta[0] + 0.5*theta[1]
>>>
>>> model = Graph.pmf_from_graph(g)
>>> theta = jnp.array([1.0, 3.0])
>>> pdf = model(theta, times)  # weight becomes 2.0*1.0 + 0.5*3.0 = 3.5
>>>
>>> # Full JAX support for parameterized graphs
>>> grad_fn = jax.grad(lambda t: jnp.sum(model(t, times)))
>>> gradient = grad_fn(theta)  # Gradients work!

For direct C++ access (no JAX overhead), use graph methods:

>>> pdf_value = g.pdf(1.5)  # Direct C++ call
>>> pmf_value = g.dph_pmf(3)  # Direct C++ call

With symbolic DAG caching (default)

>>> model = Graph.pmf_from_graph(g, use_cache=True)  # First call: computes and caches
>>> model2 = Graph.pmf_from_graph(g, use_cache=True)  # Subsequent: instant from cache!

pmf_from_graph_joint_index

phasic.Graph.pmf_from_graph_joint_index(graph, theta_dim=None, fixed_mask=None)

Create a JAX-compatible model for joint index distributions.

In joint index mode, likelihood is computed from exact expected sojourn times rather than PDF/PMF values. The observed_data should contain vertex indices (integers) instead of time values.

This is used for joint index distributions in population genetics models where the observed data represents which states (vertices) were visited.

Uses the fast expected_sojourn_time() method which computes exact sojourn times for all states in a single pass through the elimination trace.

Parameters

graph : Graph: Parameterized graph built using the Python API with parameterized edges.
param_length : int: Number of parameters for parameterized edges. If not provided, will be auto-detected from the graph.
fixed_mask : jnp.ndarray = None: Binary mask indicating which parameters are fixed (1=fixed, 0=learnable). If provided, gradients for fixed dimensions will be zero in the custom VJP. This is used to skip finite difference computation for fixed parameters.

Returns

: callable: JAX-compatible function with signature: model(theta, vertex_indices, rewards=None) -> (sojourn_times, dummy_moments) Where: - theta: Parameter vector - vertex_indices: Array of vertex indices (integers) - rewards: Ignored (must be None for joint_index mode) - sojourn_times: Expected sojourn times for the specified vertices - dummy_moments: Zeros array (moments not supported in joint_index mode)

Notes

Uses expected_sojourn_time() for fast exact computation
Much faster than iterating accumulated_visiting_time() until convergence
Moment regularization is not supported (regularization must be 0)
Reward transformation is not supported (rewards must be None)

pmf_from_graph_parameterized

phasic.Graph.pmf_from_graph_parameterized(graph_builder, discrete=False)

Convert a parameterized Python graph builder to a JAX-compatible function.

This allows users to define parameterized models where the graph structure or edge weights depend on parameters.

Parameters

graph_builder : callable: Function (theta) -> Graph that builds a graph with given parameters
discrete : bool = False: If True, uses discrete phase-type distribution (DPH) computation. If False, uses continuous phase-type distribution (PDF).

Returns

: callable: JAX-compatible function (theta, times) -> pdf_values that supports JIT, grad, vmap, etc.

Examples

>>> def build_exponential(rate):
...     g = Graph(1)
...     start = g.starting_vertex()
...     v0 = g.find_or_create_vertex([0])
...     v1 = g.find_or_create_vertex([1])
...     start.add_edge(v0, 1.0)
...     v0.add_edge(v1, float(rate))
...     return g
>>>
>>> model = Graph.pmf_from_graph_parameterized(build_exponential)
>>> theta = jnp.array([1.5])
>>> times = jnp.linspace(0, 5, 50)
>>> pdf = model(theta, times)

probability_matching

phasic.Graph.probability_matching(
    observed_data,
    fixed=None,
    theta_dim=None,
    theta_init=None,
    std_multiplier=2.0,
    verbose=True,
)

Find parameter estimates by matching model probabilities to empirical probabilities.

For joint probability graphs (created via :meth:joint_prob_graph), observations are feature-count tuples and the model outputs a probability table. This method minimises the squared difference between model and empirical probabilities.

The returned ProbMatchResult.prior list can be passed directly to :meth:Graph.svgd as the prior argument.

Parameters

observed_data : list of tuples: Observed feature-count tuples, e.g. [(0, 1, 0), (1, 0, 0), ...]. Each tuple must match a row in the joint probability table.
fixed : list = None: List of (index, value) tuples pinning specific parameters.
theta_dim : int = None: Number of model parameters. Inferred from the graph when None.
theta_init : np.ndarray = None: Initial guess for the free parameters (excluding fixed ones). If None a coordinate-wise grid search is used.
std_multiplier : float = 2.0: Factor applied to the asymptotic standard error to obtain the prior standard deviation: prior_std = std_multiplier * se.
verbose : bool = True: Print progress information.

Returns

: ProbMatchResult: Dataclass with theta, std, prior, empirical_probs, model_probs, unique_indices, etc.

Raises

: ValueError: If the graph is not a joint probability graph.

Examples

>>> jg = graph.joint_prob_graph(indexer, tot_reward_limit=2)
>>> jg.update_weights([1.0, 0.5])
>>> obs = [tuple(row) for row in jg.joint_prob_table().iloc[:, :-1].values]
>>> result = jg.probability_matching(obs, fixed=[(1, 0.5)])
>>> svgd = jg.svgd(obs, prior=result.prior, fixed=[(1, 0.5)])

reward_transform

phasic.Graph.reward_transform(rewards)

Apply reward transformation to create a new graph with modified rewards.

Parameters

rewards : np.ndarray: Reward vector of length n_vertices. Each element specifies the reward associated with visiting the corresponding vertex.

Returns

: Graph: New graph with reward-transformed distribution.

Notes

Reward transformation is used to compute moments and expectations for different reward structures. The transformation modifies the graph to compute E[∑ rewards[i] * time_in_state_i] instead of E[T].

For continuous distributions: Uses continuous reward transformation. For discrete distributions: Uses discrete reward transformation.

The returned graph is a new Graph object (not modified in-place).

reward_transform_discrete

phasic.Graph.reward_transform_discrete(rewards)

Apply reward transformation for discrete-time distributions.

Parameters

rewards : np.ndarray: Reward vector of length n_vertices.

Returns

: Graph: New graph with discrete reward transformation applied.

Notes

This method is specific to discrete-time phase-type (DPH) distributions. For automatic dispatch, use reward_transform() instead.

sample

phasic.Graph.sample(n, **kwargs)

Generate random samples from the phase-type distribution.

Parameters

n : int: Number of samples to generate.
****kwargs** : dict = {}: Additional keyword arguments passed to C++ implementation.

Returns

: np.ndarray: Array of n samples from the distribution.

Notes

Sampling is done by simulating the underlying Markov chain until absorption. For more efficient repeated sampling, first create a distribution context using distribution_context().

serialize

phasic.Graph.serialize(theta_dim=None)

Serialize graph to array representation for efficient computation.

Parameters

theta_dim : int = None: Number of parameters for parameterized edges. If not provided, will be auto-detected by probing edge states. Providing this explicitly avoids potential issues with auto-detection reading garbage memory.

Returns

: dict: Dictionary containing: - ‘states’: Array of vertex states (n_vertices, state_dim) - ‘edges’: Array of regular edges [from_idx, to_idx, weight] (n_edges, 3) - ‘start_edges’: Array of starting vertex regular edges [to_idx, weight] (n_start_edges, 2) - ‘param_edges’: Array of parameterized edges [from_idx, to_idx, x1, x2, …] (n_param_edges, theta_dim+2) - ‘start_param_edges’: Array of starting vertex parameterized edges [to_idx, x1, x2, …] (n_start_param_edges, theta_dim+1) NOTE: start_param_edges should be empty (starting edges are not parameterized) - ‘param_length’: Length of parameter vector (0 if no parameterized edges) - ‘state_length’: Integer state dimension - ‘n_vertices’: Number of vertices

state_probability

phasic.Graph.state_probability(time, **kwargs)

Compute probability of being in each state at a given time.

Parameters

time : float or int: Time point (continuous) or jump number (discrete) at which to evaluate state probabilities.
****kwargs** : dict = {}: Additional keyword arguments passed to C++ implementation.

Returns

: np.ndarray: Probability of being in each state at the specified time.

Notes

For continuous distributions: probability vector at time t For discrete distributions: probability vector after n jumps Computed via matrix exponentiation or uniformization.

stop_probability

phasic.Graph.stop_probability(time, **kwargs)

Compute probability of being in each state at a given time.

Parameters

time : float or int: Time point (continuous) or jump number (discrete) at which to evaluate state probabilities.
****kwargs** : dict = {}: Additional keyword arguments passed to C++ implementation.

Returns

: np.ndarray: Probability of being in each state at the specified time.

Notes

For continuous distributions: probability vector at time t For discrete distributions: probability vector after n jumps Computed via matrix exponentiation or uniformization.

svgd

phasic.Graph.svgd(
    observed_data,
    discrete=None,
    prior=None,
    n_particles=None,
    n_iterations=100,
    optimizer=None,
    learning_rate=None,
    bandwidth='median_per_dim',
    theta_init=None,
    theta_dim=None,
    return_history=True,
    seed=None,
    verbose=False,
    progress=True,
    jit=None,
    parallel=None,
    n_devices=None,
    precompile=True,
    compilation_config=None,
    regularization=0.0,
    nr_moments=2,
    positive_params=True,
    param_transform=None,
    joint_index=False,
    rewards=None,
    fixed=None,
    preconditioner='auto',
)

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

SVGD finds the posterior distribution p(theta | data) by optimizing a set of particles to approximate the posterior. This method works with parameterized models created by pmf_from_graph() or pmf_from_cpp() where the model signature is model(theta, times).

Parameters

observed_data : ArrayLike: Observed data points. For continuous models (PDF), these are time points where the density was observed. For discrete models (PMF), these are jump counts.
discrete : bool = None: If True, computes discrete PMF. If False, computes continuous PDF. If undefined it is inferred from the graph.is_discrete attribute.
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(ci=[0,1]), None, GaussPrior(ci=[0,1])] - DataPrior instance: Data-informed prior estimated from the observed data. If None (default), constructs DataPrior(graph, observed_data, sd=5) which uses method-of-moments (or probability matching for joint probability graphs) to estimate prior means from the data, with a wide spread (5x the asymptotic standard error). Falls back to standard normal if DataPrior construction fails. 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 = 50: Number of SVGD particles. More particles = better posterior approximation but slower.
n_iterations : int = 1000: Number of SVGD optimization steps
optimizer : object = None: Learning rate optimizer instance from phasic.optimizers. Default is Adam when learning_rate=None and regularization=0. Options include Adamelia, Adam, SGDMomentum, RMSprop, Adagrad. When an optimizer is used, the learning_rate parameter is ignored (the optimizer has its own learning rate).
learning_rate : float or None = None: SVGD step size. If None (default), uses Adame optimizer with adaptive learning rates. If a float is provided, uses fixed learning rate approach. Larger values = faster convergence but may be unstable.
bandwidth : str, float, or np.ndarray = 'median_per_dim': Kernel bandwidth selection method: - ‘median_per_dim’: Per-dimension median heuristic (default). Uses a separate bandwidth per parameter dimension for an anisotropic kernel. - ‘median’: Scalar median heuristic (isotropic kernel) - float: Fixed scalar bandwidth value - np.ndarray: Fixed per-dimension bandwidth vector
theta_init : ArrayLike = None: Initial particle positions (n_particles, theta_dim). If None, initializes randomly from standard normal.
theta_dim : int = None: Dimension of theta parameter vector. Can be: - Set at graph construction: Graph(callback, theta_dim=2) - Overridden here for SVGD inference (if graph was modified/augmented) If None, inferred from the graph’s parameterized edge structure via param_length(). Only required if theta_init is None and the graph has no parameterized edges. The value specified here overrides any theta_dim set during graph construction, which is useful if you’ve modified the graph structure (e.g., via extend()).
return_history : bool = True: If True, return particle positions throughout optimization
seed : int = None: Random seed for reproducibility
verbose : bool = True: Print progress information
jit : bool or None = None: Enable JIT compilation. If None, uses value from phasic.get_config().jit. 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)
n_devices : int or None = None: Number of devices to use for pmap. Only used when parallel=‘pmap’. If None, uses all available devices.
precompile : bool = True: (Deprecated: use jit parameter instead) Precompile model and gradient functions for faster execution. First run will take longer but subsequent iterations will be much faster.
compilation_config : CompilationConfig, dict, str, or pathlib.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)
positive_params : bool = True: If True, applies softplus transformation to ensure all parameters are positive. Recommended for phase-type models where parameters represent rates. SVGD operates in unconstrained space, but model receives positive parameters.
param_transform : callable = None: Custom parameter transformation function: transform(theta_unconstrained) -> theta_constrained. If provided, SVGD optimizes in unconstrained space and applies this transformation before calling the model. Cannot be used together with positive_params. Example: lambda theta: jnp.concatenate([jnp.exp(theta[:1]), jax.nn.softplus(theta[1:])])
joint_index : bool = False: If True, use joint index mode where observed_data contains vertex indices (integers) instead of time values. In this mode, likelihood is computed from converged accumulated_visits() values rather than PDF/PMF values. This is used for joint index distributions in population genetics models. When joint_index=True: - observed_data should contain vertex indices (integers) - Forces discrete=True behavior - Moment regularization is not supported (raises NotImplementedError if regularization > 0)
rewards : ArrayLike = None: Reward vectors for computing reward-transformed likelihoods. Can be: - None: Standard phase-type likelihood (default) - 1D array (n_vertices,): Single reward vector for univariate models - 2D array (n_vertices, n_features): Multivariate rewards - one reward vector per feature dimension. Requires use of pmf_and_moments_from_graph_multivariate() model. For multivariate models, observed_data should also be 2D (n_times, n_features).
preconditioner : str, preconditioner instance, or None = 'auto': Preconditioning method for multi-scale parameters: - ‘auto’ or ‘jacobian’: Moment Jacobian preconditioning (default, recommended). Uses column norms of the moment Jacobian matrix for scaling. Simpler and more robust than Fisher preconditioning. - ‘fisher’: Fisher diagonal preconditioning. Uses empirical Fisher information matrix diagonal. Can be unstable when PMF values are small. - None or ‘none’: No preconditioning (original behavior) - MomentJacobianPreconditioner or FisherPreconditioner instance: Custom preconditioner

Returns

: dict: Inference results containing: - ‘particles’: Final posterior samples (n_particles, theta_dim) - ‘theta_mean’: Posterior mean estimate - ‘theta_std’: Posterior standard deviation - ‘history’: Particle evolution over iterations (if return_history=True)

Raises

: ImportError: If JAX is not installed
: ValueError: If model is not parameterized or theta_dim cannot be inferred

Examples

>>> import jax.numpy as jnp
>>> from phasic import Graph
>>>
>>> # Build parameterized coalescent model
>>> def coalescent_callback(state, nr_samples=3):
...     if len(state) == 0:
...         return [(np.array([nr_samples]), 1.0, [1.0])]
...     if state[0] > 1:
...         n = state[0]
...         rate = n * (n - 1) / 2
...         return [(np.array([n - 1]), 0.0, [rate])]
...     return []
>>>
>>> g = Graph.from_callback_parameterized(coalescent_callback, nr_samples=4)
>>> model = Graph.pmf_from_graph(g, discrete=False)
>>>
>>> # Generate synthetic observed data
>>> true_theta = jnp.array([2.0])
>>> times = jnp.linspace(0.1, 3.0, 15)
>>> observed_pdf = model(true_theta, times)
>>>
>>> # Run SVGD inference
>>> results = Graph.svgd(
...     observed_data=observed_pdf,
...     theta_dim=1,
...     n_particles=30,
...     n_iterations=500,
...     learning_rate=0.01
... )
>>>
>>> print(f"True theta: {true_theta}")
>>> print(f"Posterior mean: {results['theta_mean']}")
>>> print(f"Posterior std: {results['theta_std']}")

Notes

SVGD requires a parameterized model. Non-parameterized models (signature: model(times)) cannot be used for inference as there are no parameters to estimate.
The likelihood is computed as sum(log(model(theta, observed_data)))
For better results, ensure observed_data has sufficient information about the parameters
Learning rate and number of iterations may need tuning for different problems

update_weights

phasic.Graph.update_weights(theta, callback=None, log=False)

Update parameterized edge weights with given parameters.

This method wraps the C++ implementation to cache theta for use with trace-based computation.

Parameters

theta : ArrayLike: Parameter vector to set edge weights.
callback : callable = None: Custom callback for weight computation.
log : bool = False: If True, use log-space computation.

Notes

This does NOT invalidate the cached trace since it only changes parameter values, not graph structure.

variance

phasic.Graph.variance(rewards=[], **kwargs)

Compute variance of the phase-type distribution.

Parameters

rewards : ArrayLike = []: Reward vector for reward-transformed variance. If not provided, computes Var(T) where T is time until absorption.
****kwargs** : dict = {}: Additional keyword arguments passed to C++ implementation.

Returns

: float: Variance Var(T) or reward-transformed variance.

Notes

Computed as Var(T) = E[T^2] - E[T]^2 using moments.

Graph

Attributes

Methods

absorbing_state_rewards

Returns

Examples

See Also

accumulated_occupancy

Parameters

Returns

Notes

accumulated_visiting_time

Parameters

Returns

Notes

accumulated_visits

Parameters

Returns

Raises

Notes

as_matrices

Returns

Examples

cdf

Parameters

Returns

Notes

clone

Returns

compute_trace

Parameters

Returns

Notes

Examples

copy

Returns

covariance

Parameters

Returns

Notes

discretize

Parameters

Returns

distribution_context

Parameters

Returns

Notes

eliminate_to_dag

Returns

Raises

Examples

Performance

See Also

expectation

Parameters

Returns

Notes

expected_sojourn_time

Parameters

Returns

Notes

expected_waiting_time

Parameters

Returns

Notes

extend

Parameters

Raises

Examples

find_or_create_vertex

from_matrices

Parameters

Returns

Examples

from_serialized

Parameters

Returns

Raises

Examples

laplace_transform