Dynamic Correlations with State-Space Smoothing

Import Libraries

Import the modules required for this session. Run all cells in this section before starting the exercises.

Note: This notebook uses a local libary that sits in the /ssll folder. Before we can import from it we must call sys.path.insert(0, 'ssll') to add that folder to the path where Python can find it.

import numpy as np
import matplotlib.pyplot as plt

# Import modules from local library /ssll
import sys
sys.path.insert(0, 'ssll')
import ssll
import synthesis
import transforms
import energies

When working with custom libraries, it is a good practice to test that it is working as expected. The cell below runs the code with some example paramerters and uses assert to verify the output is as expected. The cell should print: "Setup OK: shapes (5, 6) (5, 8) (5, 10, 3)"

N_test, O_test, T_test, R_test = 3, 2, 10, 20
transforms.initialise(N_test, O_test)
theta_test = synthesis.generate_thetas(N_test, O_test, T_test)
p_test = np.array([transforms.compute_p(theta_test[t]) for t in range(T_test)])
spikes_test = synthesis.generate_spikes(p_test, R_test)
assert spikes_test.shape == (T_test, R_test, N_test), f"Expected (10,20,3), got {spikes_test.shape}"
print(f"Setup OK: shapes {theta_test.shape} {p_test.shape} {spikes_test.shape}")

Setup OK: shapes (10, 6) (10, 8) (10, 20, 3)

Utility Functions

The utils class below provides two plotting helpers used throughout the exercises. Run this cell once as part of Setup — you do not need to study the implementation.

class utils:
    @staticmethod
    def plot_trajectory(theta_true, theta_est, ylabel, label='Smoothed estimate', ax=None):
        """Ground-truth vs a single estimated parameter trajectory."""
        if ax is None:
            fig, ax = plt.subplots(figsize=(10, 3))
        ax.plot(theta_true, 'k--', label='Ground truth')
        ax.plot(theta_est, label=label)
        ax.set_xlabel('Time bin')
        ax.set_ylabel(ylabel)
        ax.legend()
        plt.tight_layout()
        return ax

    @staticmethod
    def plot_trajectory_with_band(theta_true, theta_est, std, ylabel):
        """Ground-truth vs estimated trajectory with a ±2σ credible band."""
        fig, ax = plt.subplots(figsize=(10, 3))
        t_ax = np.arange(len(theta_est))
        ax.plot(theta_true, 'k--', label='Ground truth')
        ax.plot(theta_est, label=r'Smoothed $\pm 2\sigma$')
        ax.fill_between(t_ax, theta_est - 2*std, theta_est + 2*std, alpha=0.3)
        ax.set_xlabel('Time bin')
        ax.set_ylabel(ylabel)
        ax.legend()
        plt.tight_layout()

In many neural circuits, the statistical relationships between neurons change over time, i.e., activity is non-stationary. Neurons may display synchronization during a stimulus or internal brain state, and decouple at other times. To study such dynamic changes it is necessary to track how the pairwise interaction parameters $\theta$ evolve over time.

The function synthesis.generate_thetas draws a smooth trajectory of natural parameters $\theta_t$ from a Gaussian process (GP) prior, producing continuous curves for the biases and couplings. Each time step has its own $\theta$ vector, so the probability distribution $P(x|\theta_t)$ over spike patterns changes continuously. The parameter sigma1 controls the smoothness of the bias trajectories, and sigma2 controls the smoothness of the coupling trajectories. Both default to 50 time bins, meaning the parameters fluctuate on a timescale of roughly 50 bins. Smaller values produce faster, more jagged dynamics; larger values produce slower, smoother ones.

Generating spikes from a time-varying model produces a three-dimensional array of shape (T, R, N): T time bins, each with R independent trials drawn from the pattern distribution at that time bin. The trial-averaged firing rate at each time bin, obtained by spikes.mean(axis=1), reflects the instantaneous marginal probabilities and reveals how population activity changes across the recording.

Code Reference

Exercises

Example: Generate N=3, O=2, T=200, R=100. Plot biases and couplings represented by the natural parameters $\theta_t$ over time.

N, O, T, R = 3, 2, 200, 100

transforms.initialise(N, O)
theta = synthesis.generate_thetas(N, O, T, sigma1=50, sigma2=50)  # default values
p = np.array([transforms.compute_p(theta[t]) for t in range(T)])
spikes = synthesis.generate_spikes(p, R)

fig, axes = plt.subplots(2, 1, figsize=(10, 5), sharex=True)
for i in range(N):
    axes[0].plot(theta[:, i], label=f'θ{i+1}')
axes[0].set_ylabel('Bias θ')
axes[0].legend(loc='upper right')
for k in range(N * (N-1) // 2):
    axes[1].plot(theta[:, N + k], label=f'J{k+1}')
axes[1].set_ylabel('Coupling θ')
axes[1].set_xlabel('Time bin')
axes[1].legend(loc='upper right');

Exercise: Print spikes.shape of the example and verify it matches (T, R, N).

print('spikes.shape:', spikes.shape)

spikes.shape: (200, 100, 3)

Exercise: Plot the trial-averaged firing rates spikes.mean(axis=1) alongside the bias parameters theta[:, :N] (e.g., in a two-panel figure with shared x-axis). Does each neuron’s rate track its bias parameter over time?

Solution

rates = spikes.mean(axis=1)
fig, axes = plt.subplots(2, 1, figsize=(10, 5), sharex=True)
for i in range(N):
    axes[0].plot(theta[:, i], label=f'θ{i+1}')
axes[0].set_ylabel('Bias θ')
axes[0].legend(loc='upper right')
for i in range(N):
    axes[1].plot(rates[:, i], label=f'Neuron {i+1}')
axes[1].set_ylabel('Empirical rate')
axes[1].set_xlabel('Time bin')
axes[1].legend(loc='upper right')

Exercise: Regenerate with sigma1=150 (slow variation) and sigma1=20 (fast variation). Plot both bias traces (first natural parameter $\theta$) for neuron 0 on the same axis.

Solution

theta_slow = synthesis.generate_thetas(N, O, T, sigma1=150)
theta_fast = synthesis.generate_thetas(N, O, T, sigma1=20)
fig, ax = plt.subplots(figsize=(10, 3))
ax.plot(theta_slow[:, 0], label='sigma1=150')
ax.plot(theta_fast[:, 0], label='sigma1=20')
ax.set_xlabel('Time bin')
ax.set_ylabel('\u03b8\u2081')
ax.legend();

Exercise: Generate data with a stronger mean coupling (mu2=1.5); store the result as theta_strong_coupling. Plot the coupling traces and compare to the default mu2=0 case.

Solution

theta_strong_coupling = synthesis.generate_thetas(N, O, T, mu2=1.5)
fig, ax = plt.subplots(figsize=(10, 3))
ax.plot(theta[:, N:], label='Default (mu2=0)', alpha=0.6)
ax.plot(theta_strong_coupling[:, N:], label='mu2=1.5', alpha=0.6)
ax.set_xlabel('Time bin')
ax.set_ylabel('Coupling θ')
ax.legend();

Exercise: Generate stationary thetas using synthesis.generate_stationary_thetas(N, O, T). Verify that all parameters are constant over time by printing theta_stat[0] (first time bin) and theta_stat[-1] (last time bin) — they should be identical.

theta_stat = synthesis.generate_stationary_thetas(N, O, T)
print('theta at t=0:   ', theta_stat[0].round(3))
print('theta at t=T-1: ', theta_stat[-1].round(3))

theta at t=0:    [-3. -3. -3.  0.  0.  0.]
theta at t=T-1:  [-3. -3. -3.  0.  0.  0.]

The state-space model treats the natural parameters $\theta_t$ as a hidden state that evolves according to $\theta_t = \theta_{t-1} + \varepsilon_t$, where $\varepsilon_t \sim \mathcal{N}(0, Q)$ is a small random perturbation. This random-walk prior favors smooth trajectories while allowing the parameters to change over time. At each time step, the R observed spike patterns provide a noisy window into the current $\theta_t$.

The SSLL algorithm combines the prior with the data using a Kalman filter (forward pass) and a Rauch-Tung-Striebel smoother (backward pass) to produce the posterior distribution of $\theta_t$ given all $T \times R$ observations. The result emd.theta_s is the smoothed posterior mean and emd.sigma_s is the corresponding posterior covariance. The state noise covariance Q is learned from the data during the M-step.

The learned Q, stored in emd.Q, reflects how much the algorithm believes the parameters change from one time bin to the next. A larger diagonal of Q implies the model has inferred faster dynamics.

Code Reference

Exercises

Example: Fit the model and compare the smoothed bias of neuron 0 to the ground truth.

N, O, T, R = 3, 2, 200, 100

transforms.initialise(N, O)
theta = synthesis.generate_thetas(N, O, T)
p = np.array([transforms.compute_p(theta[t]) for t in range(T)])
spikes = synthesis.generate_spikes(p, R)
emd = ssll.run(spikes, order=O, EM_Info=False)

utils.plot_trajectory(theta[:, 0], emd.theta_s[:, 0], r'$\theta_1$');

Exercise: Print emd.theta_s.shape. Print the diagonal of emd.Q as a relative indicator of how fast parameters are evolving.

print('theta_s shape:', emd.theta_s.shape)
print('Learned Q (diagonal):', np.diag(emd.Q).round(5))

theta_s shape: (200, 6)
Learned Q (diagonal): [0.0084 0.0084 0.0084 0.0084 0.0084 0.0084]

Exercise: Compute the mean squared error of the smoothed estimate, separately for biases (first N columns of emd.theta_s) and couplings (remaining columns).

mse_bias = np.mean((emd.theta_s[:, :N] - theta[:, :N])**2)
mse_coupling = np.mean((emd.theta_s[:, N:] - theta[:, N:])**2)
print(f'MSE biases:    {mse_bias:.4f}')
print(f'MSE couplings: {mse_coupling:.4f}')

MSE biases:    0.0086
MSE couplings: 0.0173

Exercise: Plot the smoothed estimate of the first coupling emd.theta_s[:, N] against the ground truth theta[:, N].

Solution

utils.plot_trajectory(theta[:, N], emd.theta_s[:, N], r'$\theta_{12}$');

Exercise: Plot emd.mll_list and print the number of EM iterations.

Solution

fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(emd.mll_list)
ax.set_xlabel('EM iteration')
ax.set_ylabel('Log marginal likelihood')
print(f'Converged in {len(emd.mll_list)} iterations')

Converged in 4 iterations

Exercise: Increase R to R=300 and refit. How does the MSE change compared to R=100?

p_r300 = np.array([transforms.compute_p(theta[t]) for t in range(T)])
spikes_r300 = synthesis.generate_spikes(p_r300, 300)
emd_r300 = ssll.run(spikes_r300, order=O, EM_Info=False)
mse_r300 = np.mean((emd_r300.theta_s - theta)**2)
mse_r100 = np.mean((emd.theta_s - theta)**2)
print(f'MSE R=100: {mse_r100:.4f}')
print(f'MSE R=300: {mse_r300:.4f}')

MSE R=100: 0.0130
MSE R=300: 0.0096

The forward filter processes the data from left to right in time, estimating $\theta_t$ using only the observations up to time $t$. This is a causal estimate that could be used in real-time processing. The backward smoother then refines these estimates by propagating information from later time bins back to earlier ones, producing the non-causal smoothed estimate emd.theta_s. Because the smoother has access to future data, it is typically more accurate than the filtered estimate emd.theta_f, and the posterior covariance emd.sigma_s is correspondingly smaller.

The posterior standard deviation at each time bin quantifies the algorithm’s confidence in its estimate. At the edges of the recording, the smoother has less data to work with (before the start or after the end of the smoothing window), so posterior uncertainty is typically higher there. Comparing emd.sigma_s and emd.sigma_f reveals exactly how much the future observations reduce uncertainty about the current state.

Code Reference

Exercises

Run the cell below only if you are starting here without having completed Section 2; it regenerates N, O, T, R, theta, p, spikes, and the fitted emd.

# Realign: run this cell if you are starting Section 3 without having run Section 2
# Sets up: N, O, T, R, theta, p, spikes, emd
N, O, T, R = 3, 2, 200, 100

transforms.initialise(N, O)
theta = synthesis.generate_thetas(N, O, T)
p = np.array([transforms.compute_p(theta[t]) for t in range(T)])
spikes = synthesis.generate_spikes(p, R)
emd = ssll.run(spikes, order=O, EM_Info=False)

print(f'Realign OK — N={N}, T={T}, R={R}, emd fitted')

Realign OK — N=3, T=200, R=100, emd fitted

Example: Plot theta_f[:, 0] and theta_s[:, 0] together with the ground truth $\theta$ to compare filtered and smoothed estimates.

fig, ax = plt.subplots(figsize=(10, 3))
ax.plot(theta[:, 0], 'k--', label='Ground truth')
ax.plot(emd.theta_f[:, 0], 'g', alpha=0.7, label='Filtered')
ax.plot(emd.theta_s[:, 0], 'b', label='Smoothed')
ax.set_xlabel('Time bin')
ax.set_ylabel('\u03b8\u2081')
ax.legend();

Exercise: Compute the posterior standard deviation np.sqrt(emd.sigma_s[:, 0, 0]), store it as std_s, and add a shaded $\pm 2\sigma$ band around the smoothed estimate.

Solution

std_s = np.sqrt(emd.sigma_s[:, 0, 0])
utils.plot_trajectory_with_band(theta[:, 0], emd.theta_s[:, 0], std_s, r'$\theta_1$')

Exercise: Using std_s from the previous exercise, compute the fraction of time bins where the true $\theta_1(t)$ falls within the $\pm 2\sigma$ credible interval.

inside = np.abs(theta[:, 0] - emd.theta_s[:, 0]) < 2 * std_s
print(f'Coverage: {inside.mean()*100:.1f}% within \u00b12\u03c3')

Coverage: 100.0% within ±2σ

Exercise: Compute std_f = np.sqrt(emd.sigma_f[:, 0, 0]). Compare the mean posterior standard deviation of the filtered vs smoothed estimate.

std_f = np.sqrt(emd.sigma_f[:, 0, 0])
print(f'Mean std (smoothed): {std_s.mean():.4f}')
print(f'Mean std (filtered): {std_f.mean():.4f}')

Mean std (smoothed): 0.1236
Mean std (filtered): 0.1635

Exercise: At which time bin is the smoothed posterior standard deviation largest? Why might the edges of the recording have higher uncertainty?

# *Answer*: the smoother integrates information from both past and future observations. Near $t=0$ there is no past data and near $t=T-1$ there is no future data, so the posterior is less constrained at the edges than in the middle of the recording.
t_max = np.argmax(std_s)
print(f'Highest uncertainty at t={t_max}, std={std_s[t_max]:.4f}')

Highest uncertainty at t=199, std=0.1563

Exercise: Halve the number of trials (R=50). Regenerate the spikes and refit; store the new posterior std as std_s50. Compare std_s50.mean() to std_s.mean() from above. How do the posterior standard deviations change?

spikes_r50 = synthesis.generate_spikes(p, 50)
emd_r50 = ssll.run(spikes_r50, order=O, EM_Info=False)
std_s50 = np.sqrt(emd_r50.sigma_s[:, 0, 0])
print(f'Mean std with R=100: {std_s.mean():.4f}')
print(f'Mean std with R=50:  {std_s50.mean():.4f}')

Mean std with R=100: 0.1236
Mean std with R=50:  0.1500

The parameter state_cov sets the initial state noise covariance $Q$ and controls what the EM algorithm optimizes. A small value forces the estimated parameters to vary slowly and smoothly over time. A larger value allows the model to track rapid changes but results in noisier estimates. You can set different noise levels for biases and couplings by passing a two-element list: state_cov=[lambda_1, lambda_2] assigns lambda_1 to the $N$ bias parameters and lambda_2 to the $N(N-1)/2$ coupling parameters. Setting lambda_2=0 forces the couplings to remain constant while biases are allowed to vary.

This flexibility is useful when you have prior knowledge about the dynamics of different parameter types. For example, firing rates often drift on a slow timescale driven by neuromodulation, while pairwise correlations may be more stable during a single session. Constraining the coupling dynamics can improve inference quality when the available data is limited.

Code Reference

Exercises

Run the cell below only if you are starting here without having completed Sections 1 and 2; it regenerates N, O, T, R, theta, p, spikes, and the fitted emd.

# Realign: run this cell if you are starting Section 4 without having run Sections 1 and 2
# Sets up: N, O, T, R, theta, p, spikes, emd
N, O, T, R = 3, 2, 200, 100

transforms.initialise(N, O)
theta = synthesis.generate_thetas(N, O, T)
p = np.array([transforms.compute_p(theta[t]) for t in range(T)])
spikes = synthesis.generate_spikes(p, R)
emd = ssll.run(spikes, order=O, EM_Info=False)

print(f'Realign OK — N={N}, T={T}, R={R}, emd fitted')

Realign OK — N=3, T=200, R=100, emd fitted

Example: Compare smoothed bias estimates under two different state_cov values (e.g., 0.001 and 0.1).

emd_slow = ssll.run(spikes, order=O, state_cov=0.001, EM_Info=False)
emd_fast = ssll.run(spikes, order=O, state_cov=0.1, EM_Info=False)

fig, ax = plt.subplots(figsize=(10, 3))
ax.plot(theta[:, 0], 'k--', label='Ground truth')
ax.plot(emd_slow.theta_s[:, 0], label='state_cov=0.001')
ax.plot(emd_fast.theta_s[:, 0], label='state_cov=0.1')
ax.set_xlabel('Time bin')
ax.set_ylabel('\u03b8\u2081')
ax.legend();

Exercise: Run SSLL with state_cov=[0.01, 0.0] (dynamic biases, static couplings); store the result as emd_mixed. Plot emd_mixed.theta_s[:, N] over time against the ground truth theta[:,N] of the first coupling using the plot_trajectory utility function. Is the coupling approximately constant?

Solution

emd_mixed = ssll.run(spikes, order=O, state_cov=[0.01, 0.0], EM_Info=False)
utils.plot_trajectory(theta[:, N], emd_mixed.theta_s[:, N], r'$\theta_{12}$', label='Mixed model coupling');

Exercise: Run with state_cov=0 (fully stationary model); store the result as emd_stat. Plot the bias estimate for neuron 0 emd_stat.theta_s[:, 0] over time and verify it is flat.

Solution

emd_stat = ssll.run(spikes, order=O, state_cov=0, EM_Info=False)
utils.plot_trajectory(theta[:, 0], emd_stat.theta_s[:, 0], r'$\theta_1$', label='Stationary model');

Exercise: Compare emd.mll with emd_stat.mll. Which model achieves a higher log marginal likelihood on the time-varying data?

print(f'Dynamic model    mll: {emd.mll:.4f}')
print(f'Stationary model mll: {emd_stat.mll:.4f}')

Dynamic model    mll: -18779.3437
Stationary model mll: -18784.5708

Exercise: Generate truly stationary data with synthesis.generate_stationary_thetas. Fit both the dynamic and stationary models. Compare their mll values. Which model wins when the data is truly stationary?

theta_stat_gt = synthesis.generate_stationary_thetas(N, O, T)
p_stat_gt = np.array([transforms.compute_p(theta_stat_gt[t]) for t in range(T)])
spikes_stat_gt = synthesis.generate_spikes(p_stat_gt, R)

emd_dyn_stat = ssll.run(spikes_stat_gt, order=O, EM_Info=False)
emd_stat_stat = ssll.run(spikes_stat_gt, order=O, state_cov=0, EM_Info=False)
print(f'Dynamic on stationary data:    mll = {emd_dyn_stat.mll:.4f}')
print(f'Stationary on stationary data: mll = {emd_stat_stat.mll:.4f}')

Dynamic on stationary data:    mll = -11667.1019
Stationary on stationary data: mll = -11648.8223

The thermodynamic analysis from Session 1 can be applied for each time-bin to track how the role of pairwise interactions changes over the recording. When coupling parameters are large at time $t$, the entropy ratio $S_\text{ratio}(t)$ is high: a large fraction of the population variability is explained by pairwise interactions at that moment. When couplings are near zero, $S_\text{ratio}(t)$ is near zero and the population behaves essentially independently.

The KL divergence $D_\text{KL}(t)$ measures how much information would be lost by switching from the pairwise model to the independent model at time $t$. Together, the time courses of $S_\text{ratio}$ and $D_\text{KL}$ indicate when in the recording population interactions are most important, potentially revealing dynamics related to stimulus processing, behavioral state, or network phenomena.

Both $S_\text{ratio}$ and $D_\text{KL}$ are computed from the smoothed posterior mean emd.theta_s, so they inherit the uncertainty of the state estimation. A time bin with high posterior variance will produce less reliable thermodynamic estimates. When interpreting these quantities, it is useful to overlay the posterior uncertainty to distinguish genuine dynamics from estimation noise.

Note: this section focuses on dynamic entropy estimation ($S_\text{ratio}(t)$) and KL divergence as measures of pairwise interaction strength. Heat capacity, an equally informative thermodynamic quantity, would require additional computation and is not covered here.

Code Reference

Exercises

Run the cell below only if you are starting here without having completed Sections 1 and 2; it regenerates N, O, T, R, theta, p, spikes, and the fitted emd.

# Realign: run this cell if you are starting Section 5 without having run Sections 1 and 2
# Sets up: N, O, T, R, theta, p, spikes, emd (with thermodynamic quantities pre-computed)
N, O, T, R = 3, 2, 200, 100

transforms.initialise(N, O)
theta = synthesis.generate_thetas(N, O, T)
p = np.array([transforms.compute_p(theta[t]) for t in range(T)])
spikes = synthesis.generate_spikes(p, R)
emd = ssll.run(spikes, order=O, EM_Info=False)
energies.get_energies(emd)

print(f'Realign OK — N={N}, T={T}, R={R}, emd fitted and thermodynamic quantities computed')

Realign OK — N=3, T=200, R=100, emd fitted and thermodynamic quantities computed

Example: Compute thermodynamic quantities and plot the entropy ratio over time.

energies.get_energies(emd)

fig, ax = plt.subplots(figsize=(10, 3))
ax.plot(emd.S_ratio)
ax.set_xlabel('Time bin')
ax.set_ylabel(r'Entropy ratio $S_\mathrm{ratio}$')

Exercise: Plot emd.S1 and emd.S2 on the same axes. When is the gap between them largest?

Solution

fig, ax = plt.subplots(figsize=(10, 3))

ax.plot(emd.S1, label='S\u2081 (independent)')
ax.plot(emd.S2, label='S\u2082 (pairwise)')
ax.set_xlabel('Time bin')
ax.set_ylabel('Entropy')
ax.legend()

Exercise: Overlay the first coupling parameter emd.theta_s[:, N] and emd.S_ratio on a dual-axis plot. Do strong couplings correspond to high entropy ratios?

Solution

fig, ax1 = plt.subplots(figsize=(10, 3))
ax2 = ax1.twinx()
ax1.plot(emd.theta_s[:, N], 'b', label='θ₁₂')
ax2.plot(emd.S_ratio, 'r', alpha=0.7, label='S_ratio')
ax1.set_xlabel('Time bin')
ax1.set_ylabel('Coupling θ₁₂', color='b')
ax2.set_ylabel('Entropy ratio', color='r')
lines1, labels1 = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(lines1 + lines2, labels1 + labels2, loc='upper right')

Exercise: Plot emd.dkl over time. Print the time bin with maximum KL divergence.

Solution

fig, ax = plt.subplots(figsize=(10, 3))
ax.plot(emd.dkl)
ax.set_xlabel('Time bin')
ax.set_ylabel('KL divergence')
print(f'Maximum dkl at t={np.argmax(emd.dkl)}')

Maximum dkl at t=165

Exercise: Fit an order=1 model to the same data and call energies.get_energies; store the result as emd_ind. Verify that the entropy ratio is approximately zero for all time bins.

emd_ind = ssll.run(spikes, order=1, EM_Info=False)
energies.get_energies(emd_ind)
print(f'Order-1 model mean S_ratio: {emd_ind.S_ratio.mean():.4f}')

Order-1 model mean S_ratio: 0.0000

Exercise: Generate a population with strong mean coupling (mu2=1.5). Run ssll.run and energies.get_energies. Compare the mean emd.S_ratio to the default coupling case and confirm your expectation.

transforms.initialise(N, O)  # ssll.run(order=1) in the previous cell leaves transforms in O=1 state
theta_sc = synthesis.generate_thetas(N, O, T, mu2=1.5)
p_sc = np.array([transforms.compute_p(theta_sc[t]) for t in range(T)])
spikes_sc = synthesis.generate_spikes(p_sc, R)
emd_sc = ssll.run(spikes_sc, order=O, EM_Info=False)
energies.get_energies(emd_sc)
print(f'Default coupling mean S_ratio: {emd.S_ratio.mean():.4f}')
print(f'Strong coupling mean S_ratio:  {emd_sc.S_ratio.mean():.4f}')

Default coupling mean S_ratio: 0.0013
Strong coupling mean S_ratio:  0.1250