Source code for mixle.stats.sequences.markov_transform
"""Create, estimate, and sample from a Markov transform model for pairs of count-sets producing a third.
Defines the MarkovTransformDistribution, MarkovTransformSampler, MarkovTransformAccumulatorFactory,
MarkovTransformAccumulator, MarkovTransformEstimator, and the MarkovTransformDataEncoder classes for use with
mixle.
Data type: Tuple[List[Tuple[int, float]], List[Tuple[int, float]], List[Tuple[int, float]]]: An observation
(S1, S2, S3) consists of three bags of integer values in {0,...,W-1}, each given as a list of (value, count)
pairs. The members of S1 and S2 are drawn iid from an initial probability vector p, and each member w of S3 is
generated by picking a pair (u, v) uniformly from S1 x S2 and drawing w from the conditional distribution
P(w | u, v), stored as a sparse (W*W by W) matrix with row index u*W + v. With regularization weight alpha, the
log-likelihood is
log(P(S1, S2, S3)) = sum_{u in S1} c_u*log(p_u) + sum_{v in S2} c_v*log(p_v)
+ sum_{w in S3} c_w*log( sum_{u,v} ((1-alpha)*P(w|u,v) + alpha/W)*c_u*c_v/(n1*n2) ),
where c_u is the count attached to value u and n1, n2 are the total counts of S1 and S2. An optional length
distribution len_dist models the total counts [n1, n2, n3].
"""
import itertools
import numpy as np
from scipy.sparse import csc_matrix
from mixle.engines.arithmetic import *
from mixle.engines.arithmetic import maxrandint
from mixle.stats.compute.pdist import (
DataSequenceEncoder,
DistributionSampler,
ParameterEstimator,
SequenceEncodableProbabilityDistribution,
SequenceEncodableStatisticAccumulator,
StatisticAccumulatorFactory,
)
from mixle.stats.sequences._keyed_accumulator import InitTransKeyedAccumulator
from mixle.utils.aliasing import MISSING, coalesce_alias
from mixle.utils.optsutil import count_by_value
[docs]
class MarkovTransformDistribution(SequenceEncodableProbabilityDistribution):
"""MarkovTransformDistribution object modeling two count-sets transforming into a third count-set."""
def __init__(self, init_prob_vec, cond_prob_mat, alpha=0.0, len_dist=None):
"""MarkovTransformDistribution object.
Args:
init_prob_vec: Probability vector of length W for the values in S1 and S2.
cond_prob_mat: Sparse (W*W by W) matrix with row u*W + v holding P(. | u, v).
alpha (float): Regularization weight in [0, 1] mixed with a uniform background.
len_dist (Optional[SequenceEncodableProbabilityDistribution]): Distribution for the total counts
[n1, n2, n3]. Required for sampling.
Attributes:
init_prob_vec (np.ndarray): Probability vector for the values in S1 and S2.
cond_prob_mat (csc_matrix): Sparse (W*W by W) conditional probability matrix.
alpha (float): Regularization weight in [0, 1].
len_dist (Optional[SequenceEncodableProbabilityDistribution]): Distribution for the total counts.
num_vals (int): Number of possible values W.
"""
self.init_prob_vec = np.asarray(init_prob_vec, dtype=np.float64)
self.cond_prob_mat = csc_matrix(cond_prob_mat, dtype=np.float64)
self.len_dist = len_dist
self.num_vals = len(init_prob_vec)
self.alpha = alpha
def __str__(self):
"""Returns string representation of MarkovTransformDistribution object."""
s1 = ",".join(map(str, self.init_prob_vec))
temp = self.cond_prob_mat.nonzero()
tt = np.asarray(self.cond_prob_mat[temp[0], temp[1]]).flatten()
s20 = ",".join(map(str, tt))
s21 = ",".join(map(str, temp[0]))
s22 = ",".join(map(str, temp[1]))
s2 = "([%s], ([%s],[%s]))" % (s20, s21, s22)
s3 = str(self.alpha)
s4 = str(self.len_dist)
return "MarkovTransformDistribution([%s], %s, alpha=%s, len_dist=%s)" % (s1, s2, s3, s4)
[docs]
def density(self, x):
"""Density of the Markov transform model at observation x.
See log_density() for details.
Args:
x: Observation tuple (S1, S2, S3), each a list of (value, count) pairs.
Returns:
Density at observation x.
"""
return exp(self.log_density(x))
[docs]
def log_density(self, x):
"""Log-density of the Markov transform model at observation x.
Computes log(P(S1)) + log(P(S2)) + log(P(S3 | S1, S2)), plus the log-density of the total counts
[n1, n2, n3] under len_dist when provided. See the module docstring for the likelihood form.
Args:
x: Observation tuple (S1, S2, S3), each a list of (value, count) pairs.
Returns:
Log-density at observation x.
"""
nw = self.num_vals
a = self.alpha / nw
b = 1 - self.alpha
xx, yy, zz = x
ll1 = 0.0
ll2 = 0.0
ll3 = 0.0
n1 = 0
n2 = 0
n3 = 0
for u, c in xx:
ll1 += np.log(self.init_prob_vec[u]) * c
n1 += c
for u, c in yy:
ll2 += np.log(self.init_prob_vec[u]) * c
n2 += c
nn = n1 * n2
for w, cw in zz:
loc_ll = 0.0
for u, cu in xx:
for v, cv in yy:
loc_ll += (b * self.cond_prob_mat[u * nw + v, w] + a) * cu * cv / nn
ll3 += np.log(loc_ll) * cw
n3 += cw
rv = ll1 + ll2 + ll3
if self.len_dist is not None:
rv += self.len_dist.log_density([n1, n2, n3])
return rv
[docs]
def seq_log_density(self, x):
"""Vectorized evaluation of log-density at sequence encoded input x.
Args:
x: Encoded sequence (from seq_encode or MarkovTransformDataEncoder.seq_encode).
Returns:
Numpy array of log-densities, one per encoded observation.
"""
nw = self.num_vals
a = self.alpha / nw
b = 1 - self.alpha
rv = np.zeros(len(x[0]), dtype=np.float64)
for i, entry in enumerate(x[0]):
xx, cx, yy, cy, zz, cz = entry
ridx = np.reshape(xx * nw, (-1, 1)) + np.reshape(yy, (1, -1))
ridx = ridx.flatten()
cc = np.reshape(cx, (-1, 1)) * np.reshape(cy, (1, -1))
cc = cc.flatten()
cc /= cc.sum()
loc_cprob = self.cond_prob_mat[:, zz]
loc_cprob = ((loc_cprob[ridx, :].toarray().T) * b) + a
ll3 = np.dot(np.log(np.dot(loc_cprob, cc)), cz)
ll1 = np.dot(np.log(self.init_prob_vec[xx]), cx)
ll2 = np.dot(np.log(self.init_prob_vec[yy]), cy)
rv[i] = ll1 + ll2 + ll3
if self.len_dist is not None:
lln = self.len_dist.seq_log_density(x[1])
rv += lln
return rv
[docs]
def compute_capabilities(self):
"""Return backend capability metadata for this concrete Markov-transform instance."""
from mixle.stats.compute.capabilities import DistributionCapabilities, intersect_engine_ready
ready = intersect_engine_ready((self.len_dist,)) if self.len_dist is not None else ("numpy", "torch")
return DistributionCapabilities(engine_ready=ready, kernel_status="generic_object")
[docs]
def backend_seq_log_density(self, x, engine):
"""Engine-neutral Markov-transform scoring.
The sparse conditional-probability gather stays on the host (scipy sparse), but the
per-observation dense reductions (log-likelihood of the transform plus the two marginal
terms) run on the active engine (numpy or torch).
"""
from mixle.stats.compute.backend import backend_seq_log_density
nw = self.num_vals
a = self.alpha / nw
b = 1 - self.alpha
init_log = engine.log(engine.asarray(self.init_prob_vec))
vals = []
for entry in x[0]:
xx, cx, yy, cy, zz, cz = entry
ridx = (np.reshape(xx * nw, (-1, 1)) + np.reshape(yy, (1, -1))).flatten()
cc = (np.reshape(cx, (-1, 1)) * np.reshape(cy, (1, -1))).flatten()
cc = cc / cc.sum()
loc_dense = (self.cond_prob_mat[ridx, :][:, zz]).toarray().T # (len(zz), len(ridx))
loc = engine.asarray(loc_dense) * b + a
inner = engine.matmul(loc, engine.asarray(cc)) # (len(zz),)
ll3 = engine.sum(engine.log(inner) * engine.asarray(np.asarray(cz, dtype=np.float64)))
ll1 = engine.sum(
init_log[engine.asarray(np.asarray(xx, dtype=np.int64))]
* engine.asarray(np.asarray(cx, dtype=np.float64))
)
ll2 = engine.sum(
init_log[engine.asarray(np.asarray(yy, dtype=np.int64))]
* engine.asarray(np.asarray(cy, dtype=np.float64))
)
vals.append(float(engine.to_numpy(ll1 + ll2 + ll3)))
rv = engine.asarray(np.asarray(vals, dtype=np.float64))
if self.len_dist is not None:
rv = rv + backend_seq_log_density(self.len_dist, x[1], engine)
return rv
[docs]
def seq_encode(self, x):
"""Encode a sequence of observations for vectorized calls (legacy method).
Note: this legacy method encodes the lengths with self.len_dist.seq_encode(); prefer
dist_to_encoder().seq_encode(x), which uses the length distribution's DataSequenceEncoder.
Args:
x: Sequence of observation tuples (S1, S2, S3), each a list of (value, count) pairs.
Returns:
Tuple (rv, nn, vv) where rv holds per-observation (values, counts) arrays for S1, S2, S3, nn is the
encoded length data (None if len_dist is None), and vv is the array of distinct (u, v, w) triples.
"""
rv = []
nn = []
vset = set()
for xx in x:
rv0 = []
nn0 = []
for cvec in xx:
rv0.append(np.asarray([v for v, c in cvec], dtype=int))
rv0.append(np.asarray([c for v, c in cvec], dtype=float))
nn0.append(np.sum(rv0[-1]))
vset.update(itertools.product(rv0[0], rv0[2], rv0[4]))
rv.append(tuple(rv0))
nn.append(tuple(nn0))
if self.len_dist is not None:
nn = self.len_dist.seq_encode(nn)
else:
nn = None
vv = np.zeros((len(vset), 3), dtype=int)
for i, vvv in enumerate(vset):
vv[i, :] = vvv[:]
return rv, nn, vv
[docs]
def sampler(self, seed=None):
"""Create a MarkovTransformSampler object from this instance.
Requires len_dist to be set (it samples the total counts [n1, n2, n3]).
Args:
seed (Optional[int]): Used to set seed in random sampler.
Returns:
MarkovTransformSampler object.
"""
return MarkovTransformSampler(self, seed)
[docs]
def estimator(self, pseudo_count=None):
"""Create a MarkovTransformEstimator object from this instance.
Args:
pseudo_count (Optional[float]): Used to inflate sufficient statistics.
Returns:
MarkovTransformEstimator object.
"""
len_estimator = None if self.len_dist is None else self.len_dist.estimator(pseudo_count)
return MarkovTransformEstimator(
self.num_vals, alpha=self.alpha, len_estimator=len_estimator, pseudo_count=pseudo_count
)
[docs]
def dist_to_encoder(self):
"""Returns a MarkovTransformDataEncoder object for encoding sequences of data."""
len_encoder = None if self.len_dist is None else self.len_dist.dist_to_encoder()
return MarkovTransformDataEncoder(len_encoder=len_encoder)
[docs]
class MarkovTransformSampler(DistributionSampler):
"""MarkovTransformSampler object for sampling observations from a MarkovTransformDistribution."""
def __init__(self, dist: MarkovTransformDistribution, seed: int | None = None):
"""MarkovTransformSampler object.
Args:
dist (MarkovTransformDistribution): Distribution to sample from. Must have len_dist set.
seed (Optional[int]): Used to set seed in random sampler.
Attributes:
dist (MarkovTransformDistribution): Distribution to sample from.
rng (RandomState): RandomState with seed set if passed in args.
size_sampler (DistributionSampler): Sampler for the total counts [n1, n2, n3].
"""
self.rng = np.random.RandomState(seed)
self.dist = dist
# self.init_sampler = np.random.RandomState(self.rng.tomaxint())
# self.next_sampler = np.random.RandomState(self.rng.tomaxint())
# self.tran_sampler = np.random.RandomState(self.rng.tomaxint())
# self.flat_sampler = np.random.RandomState(self.rng.tomaxint())
self.size_sampler = self.dist.len_dist.sampler(seed=self.rng.randint(0, maxrandint))
[docs]
def sample(self, size: int | None = None):
"""Draw 'size' iid observations from the Markov transform model.
Each observation is a tuple (S1, S2, S3) of lists of (value, count) pairs. If size is None a single
observation is returned, else a list of 'size' observations is returned.
Args:
size (Optional[int]): Number of observations to draw. Treated as a single draw if None.
Returns:
A single observation tuple, or a list of observation tuples when size is not None.
"""
if size is None:
slens = self.size_sampler.sample()
rng = np.random.RandomState(self.rng.randint(0, maxrandint))
v1 = list(rng.choice(len(self.dist.init_prob_vec), p=self.dist.init_prob_vec, replace=True, size=slens[0]))
v2 = list(rng.choice(len(self.dist.init_prob_vec), p=self.dist.init_prob_vec, replace=True, size=slens[1]))
v3 = []
z1 = list(rng.choice(len(v1), replace=True, size=slens[2]))
z2 = list(rng.choice(len(v2), replace=True, size=slens[2]))
nw = self.dist.num_vals
for zz1, zz2 in zip(z1, z2):
if rng.rand() > self.dist.alpha:
p = self.dist.cond_prob_mat[v1[zz1] * nw + v2[zz2], :].toarray().flatten()
v3.append(rng.choice(nw, p=p))
else:
v3.append(rng.choice(nw))
return list(count_by_value(v1).items()), list(count_by_value(v2).items()), list(count_by_value(v3).items())
else:
return [self.sample() for i in range(size)]
[docs]
class MarkovTransformAccumulator(InitTransKeyedAccumulator, SequenceEncodableStatisticAccumulator):
"""MarkovTransformAccumulator object for accumulating sufficient statistics of the Markov transform model."""
def __init__(self, num_vals, size_acc=None, keys=(None, None)):
"""MarkovTransformAccumulator object.
Args:
num_vals (int): Number of possible values W.
size_acc (Optional[SequenceEncodableStatisticAccumulator]): Accumulator for the total counts.
keys (Tuple[Optional[str], Optional[str]]): Keys for initial and transition statistics.
Attributes:
init_count (np.ndarray): Weighted counts for the initial probability vector.
trans_count (csc_matrix): Weighted (W*W by W) counts for the conditional probability matrix.
size_accumulator (Optional[SequenceEncodableStatisticAccumulator]): Accumulator for total counts.
num_vals (int): Number of possible values W.
init_key (Optional[str]): Key for the initial-count statistics.
trans_key (Optional[str]): Key for the transition-count statistics.
"""
self.init_count = np.zeros(num_vals)
self.trans_count = csc_matrix((num_vals * num_vals, num_vals))
self.size_accumulator = size_acc
self.num_vals = num_vals
self.init_key = keys[0]
self.trans_key = keys[1]
# Data log-likelihood accumulated as a byproduct of the E-step (the per-observation
# log_density), only when _track_ll is enabled. Used by the fused-EM fast path in
# optimize(reuse_estep_ll=True); not part of value(). Off by default so the standard path
# pays nothing.
self._track_ll = False
self._seq_ll = 0.0
[docs]
def update(self, x, weight, estimate):
"""Update sufficient statistics with a single weighted observation.
Args:
x: Observation tuple (S1, S2, S3), each a list of (value, count) pairs.
weight (float): Weight of the observation.
estimate (MarkovTransformDistribution): Previous estimate used to assign transition responsibility.
Returns:
None.
"""
nw = self.num_vals
xx, yy, zz = x
vx = np.asarray([u[0] for u in xx], dtype=int)
cx = np.asarray([u[1] for u in xx], dtype=float)
vy = np.asarray([u[0] for u in yy], dtype=int)
cy = np.asarray([u[1] for u in yy], dtype=float)
vz = np.asarray([u[0] for u in zz], dtype=int)
cz = np.asarray([u[1] for u in zz], dtype=float)
ridx = np.reshape(vx * nw, (-1, 1)) + np.reshape(vy, (1, -1))
ridx = ridx.flatten()[:, None]
cc = np.reshape(cx, (-1, 1)) * np.reshape(cy, (1, -1))
cc = cc.flatten()[:, None]
cs = cc.sum()
a = estimate.alpha / nw
b = 1 - estimate.alpha
temp = estimate.cond_prob_mat[ridx, vz].toarray()
loc_cprob = temp * cc
w = loc_cprob.sum(axis=0)
loc_cprob *= (cz * b / (b * w + a * cs)) * weight
self.trans_count[ridx, vz] += loc_cprob
self.init_count[vx] += cx * weight
self.init_count[vy] += cy * weight
if self.size_accumulator is not None:
self.size_accumulator.update((cx.sum(), cy.sum(), cz.sum()), weight, estimate.len_dist)
[docs]
def initialize(self, x, weight, rng):
"""Initialize sufficient statistics with a single weighted observation (no previous estimate).
Args:
x: Observation tuple (S1, S2, S3), each a list of (value, count) pairs.
weight (float): Weight of the observation.
rng (RandomState): Used to initialize the size accumulator if present.
Returns:
None.
"""
nw = self.num_vals
xx, yy, zz = x
vx = np.asarray([u[0] for u in xx], dtype=int)
cx = np.asarray([u[1] for u in xx], dtype=float)
vy = np.asarray([u[0] for u in yy], dtype=int)
cy = np.asarray([u[1] for u in yy], dtype=float)
vz = np.asarray([u[0] for u in zz], dtype=int)
cz = np.asarray([u[1] for u in zz], dtype=float)
ridx = np.reshape(vx * nw, (-1, 1)) + np.reshape(vy, (1, -1))
ridx = ridx.flatten()[:, None]
cc = np.reshape(cx, (-1, 1)) * np.reshape(cy, (1, -1))
cc = cc.flatten()
loc_cprob = np.outer(cc / cc.sum(), weight * cz)
# umat = lil_matrix((nw * nw, nw))
# umat[ridx, vz] = loc_cprob
self.trans_count[ridx, vz] += loc_cprob
# self.trans_count += umat
self.init_count[vx] += cx * weight
self.init_count[vy] += cy * weight
if self.size_accumulator is not None:
self.size_accumulator.initialize((cx.sum(), cy.sum(), cz.sum()), weight, rng)
[docs]
def seq_initialize(self, x, weights, rng):
"""Initialize sufficient statistics with a sequence of weighted encoded observations.
Applies the same updates as initialize() to each encoded observation.
Args:
x: Encoded sequence (from MarkovTransformDataEncoder.seq_encode).
weights (np.ndarray): Weights, one per encoded observation.
rng (RandomState): Used to initialize the size accumulator if present.
Returns:
None.
"""
nw = self.num_vals
for entry, ww in zip(x[0], weights):
xx, cx, yy, cy, zz, cz = entry
ridx = np.reshape(xx * nw, (-1, 1)) + np.reshape(yy, (1, -1))
ridx = ridx.flatten()[:, None]
cc = np.reshape(cx, (-1, 1)) * np.reshape(cy, (1, -1))
cc = cc.flatten()
loc_cprob = np.outer(cc / cc.sum(), ww * cz)
self.trans_count[ridx, zz] += loc_cprob
self.init_count[xx] += cx * ww
self.init_count[yy] += cy * ww
if self.size_accumulator is not None:
self.size_accumulator.seq_initialize(x[1], weights, rng)
[docs]
def seq_update(self, x, weights, estimate):
"""Update sufficient statistics with a sequence of weighted encoded observations.
Args:
x: Encoded sequence (from MarkovTransformDataEncoder.seq_encode).
weights (np.ndarray): Weights, one per encoded observation.
estimate (MarkovTransformDistribution): Previous estimate used to assign transition responsibility.
Returns:
None.
"""
nw = self.num_vals
nzv = x[2]
a = estimate.alpha / nw
b = 1 - estimate.alpha
umat = csc_matrix((np.zeros(nzv.shape[0]), (nzv[:, 0] * nw + nzv[:, 1], nzv[:, 2])), shape=(nw * nw, nw))
track = self._track_ll
log_init = np.log(estimate.init_prob_vec) if track else None
obs_ll = np.zeros(len(x[0]), dtype=np.float64) if track else None
for i, (entry, ww) in enumerate(zip(x[0], weights)):
xx, cx, yy, cy, zz, cz = entry
ridx = np.reshape(xx * nw, (-1, 1)) + np.reshape(yy, (1, -1))
ridx = ridx.flatten()[:, None]
cc = np.reshape(cx, (-1, 1)) * np.reshape(cy, (1, -1))
cc = cc.flatten()[:, None]
cs = cc.sum()
temp = estimate.cond_prob_mat[ridx, zz].toarray()
loc_cprob = temp * cc
w = loc_cprob.sum(axis=0)
if track:
# Per-observation log-density (== MarkovTransformDistribution.log_density), reusing
# ``w`` (== sum_{u,v} P(z|u,v)*c_u*c_v) and ``cs`` before the responsibility scaling
# below. seq_log_density normalizes cc by cs, so inner = (b*w + a*cs)/cs.
with np.errstate(divide="ignore"):
ll3 = float(np.dot(np.log((b * w + a * cs) / cs), cz))
ll1 = float(np.dot(log_init[xx], cx))
ll2 = float(np.dot(log_init[yy], cy))
obs_ll[i] = ll1 + ll2 + ll3
loc_cprob *= (cz * b / (b * w + a * cs)) * ww
umat[ridx, zz] += loc_cprob
self.init_count[xx] += cx * ww
self.init_count[yy] += cy * ww
if self.size_accumulator is not None:
self.size_accumulator.seq_update(x[1], weights, estimate.len_dist)
if track:
if self.size_accumulator is not None:
obs_ll += estimate.len_dist.seq_log_density(x[1])
self._seq_ll += float(np.dot(np.asarray(weights, dtype=np.float64), obs_ll))
self.trans_count += umat
[docs]
def seq_update_engine(self, x, weights, estimate, engine):
"""Engine-aware E-step. The per-observation transition responsibilities are computed on the
active engine (numpy or torch); the sparse conditional gather and the sparse count scatter
stay on the host, since the sufficient statistic is a sparse matrix. Mirrors seq_update.
"""
nw = self.num_vals
nzv = x[2]
a = estimate.alpha / nw
b = 1 - estimate.alpha
weights_np = np.asarray(engine.to_numpy(weights) if hasattr(engine, "to_numpy") else weights, dtype=np.float64)
umat = csc_matrix((np.zeros(nzv.shape[0]), (nzv[:, 0] * nw + nzv[:, 1], nzv[:, 2])), shape=(nw * nw, nw))
for i, (entry, ww) in enumerate(zip(x[0], weights_np)):
xx, cx, yy, cy, zz, cz = entry
ridx = (np.reshape(xx * nw, (-1, 1)) + np.reshape(yy, (1, -1))).flatten()[:, None]
cc = (np.reshape(cx, (-1, 1)) * np.reshape(cy, (1, -1))).flatten()[:, None]
cs = float(cc.sum())
temp = estimate.cond_prob_mat[ridx, zz].toarray() # (len(ridx), len(zz))
loc_cprob = engine.asarray(temp) * engine.asarray(cc)
w = engine.sum(loc_cprob, axis=0) # (len(zz),)
scale = engine.asarray(np.asarray(cz, dtype=np.float64)) * b / (b * w + a * cs) * float(ww)
loc_cprob = loc_cprob * scale[None, :]
umat[ridx, zz] += np.asarray(engine.to_numpy(loc_cprob))
self.init_count[xx] += cx * ww
self.init_count[yy] += cy * ww
if self.size_accumulator is not None:
self.size_accumulator.seq_update(x[1], weights_np, estimate.len_dist)
self.trans_count += umat
[docs]
def combine(self, suff_stat):
"""Merge the sufficient statistics of another accumulator into this one.
Args:
suff_stat: Tuple (init_count, trans_count, size_value) from another accumulator's value().
Returns:
This MarkovTransformAccumulator object.
"""
init_count, trans_count, size_acc = suff_stat
if self.size_accumulator is not None:
self.size_accumulator.combine(size_acc)
self.init_count += init_count
self.trans_count += trans_count
return self
[docs]
def value(self):
"""Returns the sufficient statistic tuple (init_count, trans_count, size_value)."""
if self.size_accumulator is not None:
return self.init_count, self.trans_count, self.size_accumulator.value()
else:
return self.init_count, self.trans_count, None
[docs]
def from_value(self, x):
"""Set the sufficient statistics from a value() tuple.
Args:
x: Tuple (init_count, trans_count, size_value).
Returns:
This MarkovTransformAccumulator object.
"""
init_count, trans_count, size_acc = x
self.init_count = init_count
self.trans_count = trans_count
if self.size_accumulator is not None:
self.size_accumulator.from_value(size_acc)
return self
# key_merge / key_replace: provided by InitTransKeyedAccumulator (shared two-key plumbing).
[docs]
def acc_to_encoder(self):
"""Returns a MarkovTransformDataEncoder object for encoding sequences of data."""
len_encoder = None if self.size_accumulator is None else self.size_accumulator.acc_to_encoder()
return MarkovTransformDataEncoder(len_encoder=len_encoder)
[docs]
class MarkovTransformAccumulatorFactory(StatisticAccumulatorFactory):
"""MarkovTransformAccumulatorFactory object for creating MarkovTransformAccumulator objects."""
def __init__(self, num_vals, len_factory, keys):
"""MarkovTransformAccumulatorFactory object.
Args:
num_vals (int): Number of possible values W.
len_factory (Optional[StatisticAccumulatorFactory]): Factory for the total-count accumulator.
keys (Tuple[Optional[str], Optional[str]]): Keys for initial and transition statistics.
Attributes:
num_vals (int): Number of possible values W.
len_factory (Optional[StatisticAccumulatorFactory]): Factory for the total-count accumulator.
keys (Tuple[Optional[str], Optional[str]]): Keys for initial and transition statistics.
"""
self.len_factory = len_factory
self.keys = keys
self.num_vals = num_vals
[docs]
def make(self):
"""Returns a new MarkovTransformAccumulator object."""
if self.len_factory is None:
return MarkovTransformAccumulator(self.num_vals, size_acc=None, keys=self.keys)
else:
return MarkovTransformAccumulator(self.num_vals, size_acc=self.len_factory.make(), keys=self.keys)
[docs]
class MarkovTransformEstimator(ParameterEstimator):
"""MarkovTransformEstimator object for estimating MarkovTransformDistribution objects from statistics."""
def __init__(
self,
num_vals=MISSING,
alpha=0.0,
len_estimator=None,
suff_stat=None,
pseudo_count=None,
keys=(None, None),
num_values=MISSING,
):
"""MarkovTransformEstimator object.
Args:
num_vals (int): Number of possible values W.
alpha (float): Regularization weight in [0, 1] for the estimated distribution.
len_estimator (Optional[ParameterEstimator]): Estimator for the total counts [n1, n2, n3].
suff_stat (Optional[Any]): Kept for consistency with the estimate function.
pseudo_count (Optional[float]): Kept for consistency (unused in estimation).
keys (Tuple[Optional[str], Optional[str]]): Keys for initial and transition statistics.
Attributes:
num_vals (int): Number of possible values W.
alpha (float): Regularization weight in [0, 1] for the estimated distribution.
len_estimator (Optional[ParameterEstimator]): Estimator for the total counts [n1, n2, n3].
suff_stat (Optional[Any]): Kept for consistency with the estimate function.
pseudo_count (Optional[float]): Kept for consistency (unused in estimation).
keys (Tuple[Optional[str], Optional[str]]): Keys for initial and transition statistics.
"""
self.keys = keys
self.len_estimator = len_estimator
self.pseudo_count = pseudo_count
self.suff_stat = suff_stat
self.num_vals = coalesce_alias("num_vals", num_vals, "num_values", num_values, default=MISSING)
self.alpha = alpha
[docs]
def accumulator_factory(self):
"""Returns a MarkovTransformAccumulatorFactory object for this estimator."""
len_factory = None if self.len_estimator is None else self.len_estimator.accumulator_factory()
return MarkovTransformAccumulatorFactory(self.num_vals, len_factory, self.keys)
[docs]
def accumulatorFactory(self):
"""Deprecated alias for accumulator_factory()."""
return self.accumulator_factory()
[docs]
def estimate(self, nobs, suff_stat):
"""Estimate a MarkovTransformDistribution from aggregated sufficient statistics.
Arg suff_stat is a tuple of length 3 containing:
suff_stat[0] (np.ndarray): Weighted counts for the initial probability vector.
suff_stat[1] (csc_matrix): Weighted (W*W by W) counts for the conditional probability matrix.
suff_stat[2]: Sufficient statistics for the total-count distribution (None if not tracked).
Args:
nobs (Optional[float]): Weighted number of observations.
suff_stat: See above for details.
Returns:
MarkovTransformDistribution object.
"""
init_count, trans_count, size_stats = suff_stat
if self.len_estimator is not None:
len_dist = self.len_estimator.estimate(nobs, size_stats)
else:
len_dist = None
trans_count = trans_count.tocsc()
row_sum = trans_count * csc_matrix(np.ones((trans_count.shape[1], 1)))
init_prob = init_count / np.sum(init_count)
trans_prob = trans_count.multiply(row_sum.power(-1))
return MarkovTransformDistribution(init_prob, trans_prob, self.alpha, len_dist)
[docs]
class MarkovTransformDataEncoder(DataSequenceEncoder):
"""MarkovTransformDataEncoder object for encoding sequences of Markov transform observations."""
def __init__(self, len_encoder=None):
"""MarkovTransformDataEncoder object.
Args:
len_encoder (Optional[DataSequenceEncoder]): Encoder for the total counts [n1, n2, n3].
Attributes:
len_encoder (Optional[DataSequenceEncoder]): Encoder for the total counts [n1, n2, n3].
"""
self.len_encoder = len_encoder
def __str__(self):
"""Returns string representation of MarkovTransformDataEncoder object."""
return "MarkovTransformDataEncoder(len_encoder=%s)" % (str(self.len_encoder))
def __eq__(self, other):
"""Encoders are interchangeable iff other is a MarkovTransformDataEncoder with an equal len_encoder.
Args:
other (object): Object to compare against.
Returns:
True if other is an equivalent MarkovTransformDataEncoder instance.
"""
if isinstance(other, MarkovTransformDataEncoder):
return other.len_encoder == self.len_encoder
else:
return False
[docs]
def seq_encode(self, x):
"""Encode a sequence of observations for vectorized calls.
Args:
x: Sequence of observation tuples (S1, S2, S3), each a list of (value, count) pairs.
Returns:
Tuple (rv, nn, vv) where rv holds per-observation (values, counts) arrays for S1, S2, S3, nn is the
encoded length data (None if len_encoder is None), and vv is the array of distinct (u, v, w) triples.
"""
rv = []
nn = []
vset = set()
for xx in x:
rv0 = []
nn0 = []
for cvec in xx:
rv0.append(np.asarray([v for v, c in cvec], dtype=int))
rv0.append(np.asarray([c for v, c in cvec], dtype=float))
nn0.append(np.sum(rv0[-1]))
vset.update(itertools.product(rv0[0], rv0[2], rv0[4]))
rv.append(tuple(rv0))
nn.append(tuple(nn0))
if self.len_encoder is not None:
nn = self.len_encoder.seq_encode(nn)
else:
nn = None
vv = np.zeros((len(vset), 3), dtype=int)
for i, vvv in enumerate(vset):
vv[i, :] = vvv[:]
return rv, nn, vv