Installer avec Codex ou Claude Copiez ce prompt, collez-le dans Codex, Claude ou un autre assistant, puis laissez-le vérifier la page du skill et l'installer pour vous.
Une commande directe contourne le prompt de vérification. Examinez la source avant de l'exécuter.
This skill covers core quantitative methods for complex systems research: power-law
distribution testing, long-range correlations (Hurst exponent), fractal geometry,
agent-based modelling with Mesa, sample entropy, bond percolation, and
information-theoretic measures.
The powerlaw package fits power-law distributions using maximum-likelihood estimation
and performs likelihood-ratio tests against alternative heavy-tailed distributions.
import powerlaw
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
deffit_power_law(data: np.ndarray, xmin: float = None, discrete: bool = False) -> dict:
"""
Fit a power-law distribution to empirical data using MLE.
Parameters
----------
data : np.ndarray — positive values (e.g., degree sequence, word frequencies)
xmin : float or None — minimum value for fitting; estimated automatically if None
discrete : bool — True for integer data (e.g., degree counts)
Returns
-------
dict with keys: alpha, xmin, sigma (standard error), n_tail, fit object
"""
data = np.asarray(data, dtype=float)
data = data[data > 0]
fit = powerlaw.Fit(data, xmin=xmin, discrete=discrete, verbose=False)
return {
"alpha": fit.alpha,
"xmin": fit.xmin,
"sigma": fit.sigma,
"n_tail": int(np.sum(data >= fit.xmin)),
"n_total": (data),
: fit,
}
() -> pd.DataFrame:
data = np.asarray(data, dtype=)
data = data[data > ]
fit = powerlaw.Fit(data, discrete=discrete, verbose=)
alternatives = [, , , ]
records = []
alt alternatives:
:
R, p = fit.distribution_compare(, alt, normalized_ratio=)
preferred = (R > p < ) (alt (R < p < ) )
records.append({: alt, : (R, ), : (p, ), : preferred})
Exception exc:
records.append({: alt, : np.nan, : np.nan, : })
pd.DataFrame(records)
():
fit = fit_result[]
ax :
fig, ax = plt.subplots(figsize=(, ))
fit.plot_ccdf(ax=ax, color=, linewidth=, label=)
fit.power_law.plot_ccdf(ax=ax, color=, linestyle=,
label=)
fit.lognormal.plot_ccdf(ax=ax, color=, linestyle=,
label=)
ax.set_xscale(); ax.set_yscale()
ax.set_xlabel(); ax.set_ylabel()
ax.set_title(title)
ax.legend(fontsize=)
ax
len
"fit"
def
test_power_law_vs_alternatives
data: np.ndarray, discrete: bool = False
"""
Test whether a power law is a better fit than lognormal, exponential,
stretched exponential (Weibull), and truncated power law.
Uses the Vuong likelihood-ratio test (p < 0.05 favours the first distribution
if R > 0, the second if R < 0).
Returns
-------
pd.DataFrame with columns: alternative, R (log-likelihood ratio), p_value, preferred
"""
"""
Plot empirical CCDF with the fitted power-law (and lognormal comparison) overlay.
"""
"fit"
if
is
None
6
5
"steelblue"
0.5
"Empirical CCDF"
"red"
"--"
f"Power law α={fit.alpha:.2f}"
"green"
":"
"Lognormal"
"log"
"log"
"x"
"P(X ≥ x)"
9
return
2. Hurst Exponent via R/S Analysis
The Hurst exponent H characterises long-range correlations:
H ≈ 0.5: random walk / white noise
H > 0.5: persistent (trending) process
H < 0.5: anti-persistent (mean-reverting) process
defcompute_hurst(ts: np.ndarray, min_window: int = 8) -> dict:
"""
Estimate the Hurst exponent using rescaled range (R/S) analysis.
Parameters
----------
ts : np.ndarray — 1-D time series
min_window : int — smallest window size (must be ≥ 4)
Returns
-------
dict with keys: hurst, intercept, windows, rs_values, r_squared
"""
ts = np.asarray(ts, dtype=float)
n = len(ts)
# Generate window sizes as powers of 2
max_power = int(np.log2(n)) - 1
windows = [2 ** p for p inrange(int(np.log2(min_window)), max_power + 1)]
rs_values = []
for w in windows:
rs_per_window = []
for start inrange(0, n - w + 1, w):
sub = ts[start: start + w]
mean_sub = np.mean(sub)
deviation = np.cumsum(sub - mean_sub)
R = np.max(deviation) - np.min(deviation)
S = np.std(sub, ddof=1)
if S > 0:
rs_per_window.append(R / S)
if rs_per_window:
rs_values.append(np.mean(rs_per_window))
else:
rs_values.append(np.nan)
windows_arr = np.array(windows, dtype=float)
rs_arr = np.array(rs_values, dtype=float)
valid = ~np.isnan(rs_arr) & (rs_arr > 0)
log_w = np.log(windows_arr[valid])
log_rs = np.log(rs_arr[valid])
coeffs = np.polyfit(log_w, log_rs, 1)
hurst = coeffs[0]
rs_pred = np.polyval(coeffs, log_w)
ss_res = np.sum((log_rs - rs_pred) ** 2)
ss_tot = np.sum((log_rs - np.mean(log_rs)) ** 2)
r2 = 1 - ss_res / ss_tot if ss_tot > 0else0.0return {
"hurst": round(hurst, 4),
"intercept": coeffs[1],
"windows": windows_arr[valid].tolist(),
"rs_values": rs_arr[valid].tolist(),
"r_squared": round(r2, 4),
}
3. Fractal Dimension via Box-Counting
defestimate_fractal_dimension(array_2d: np.ndarray, threshold: float = None) -> dict:
"""
Estimate the fractal (Hausdorff) dimension of a binary 2-D pattern via
box-counting.
Parameters
----------
array_2d : np.ndarray — 2-D float or binary array
threshold : float or None — binarisation threshold (default: mean of array)
Returns
-------
dict with keys: fractal_dimension, r_squared, box_sizes, counts
"""
arr = np.asarray(array_2d, dtype=float)
if threshold isNone:
threshold = np.mean(arr)
binary = (arr >= threshold).astype(int)
# Pad to next power of 2 for clean box counting
max_dim = max(binary.shape)
next_pow2 = 2 ** int(np.ceil(np.log2(max_dim)))
padded = np.zeros((next_pow2, next_pow2), dtype=int)
padded[: binary.shape[0], : binary.shape[1]] = binary
box_sizes = []
counts = []
size = next_pow2
while size >= 2:
# Count non-empty boxes of given size
n_boxes_per_side = next_pow2 // size
reshaped = padded.reshape(n_boxes_per_side, size, n_boxes_per_side, size)
box_sums = reshaped.sum(axis=(1, 3))
count = np.sum(box_sums > 0)
if count > 0:
box_sizes.append(size)
counts.append(count)
size //= 2
log_s = np.log(1.0 / np.array(box_sizes, dtype=float))
log_c = np.log(np.array(counts, dtype=float))
coeffs = np.polyfit(log_s, log_c, 1)
fd = coeffs[0]
c_pred = np.polyval(coeffs, log_s)
ss_res = np.sum((log_c - c_pred) ** 2)
ss_tot = np.sum((log_c - np.mean(log_c)) ** 2)
r2 = 1 - ss_res / ss_tot if ss_tot > 0else0.0return {
"fractal_dimension": round(fd, 4),
"r_squared": round(r2, 4),
"box_sizes": box_sizes,
"counts": counts,
}
4. Agent-Based Modeling: Schelling Segregation
import random
from mesa import Agent, Model
from mesa.time import RandomActivation
from mesa.space import SingleGrid
from mesa.datacollection import DataCollector
classSchellingAgent(Agent):
"""A single household agent in the Schelling segregation model."""def__init__(self, unique_id: int, model: "SchellingModel", agent_type: int):
super().__init__(unique_id, model)
self.agent_type = agent_type # 0 or 1 (two groups)self.is_happy = Falsedefstep(self):
neighbours = self.model.grid.get_neighbors(self.pos, moore=True, include_center=False)
ifnot neighbours:
self.is_happy = Truereturn
same_type = sum(1for n in neighbours if n.agent_type == self.agent_type)
self.is_happy = same_type / len(neighbours) >= self.model.tolerance
ifnotself.is_happy:
self.model.grid.move_to_empty(self)
classSchellingModel(Model):
"""Schelling (1971) segregation model."""def__init__(self, n: int = 20, density: float = 0.9, tolerance: float = 0.3, seed: int = 42):
super().__init__()
self.n = n
self.density = density
self.tolerance = tolerance
self.schedule = RandomActivation(self)
self.grid = SingleGrid(n, n, torus=True)
self.running = Trueself.datacollector = DataCollector(
model_reporters={
"pct_happy": lambda m: sum(a.is_happy for a in m.schedule.agents) / m.schedule.get_agent_count(),
"segregation_index": lambda m: m.compute_segregation(),
}
)
agent_id = 0for cell inself.grid.coord_iter():
_, x, y = cell
ifself.random.random() < density:
agent_type = self.random.choice([0, 1])
agent = SchellingAgent(agent_id, self, agent_type)
self.grid.place_agent(agent, (x, y))
self.schedule.add(agent)
agent_id += 1defcompute_segregation(self) -> float:
"""
Moran's I-style segregation index:
fraction of same-type neighbours averaged over all agents.
"""
fractions = []
for agent inself.schedule.agents:
neighbours = self.grid.get_neighbors(agent.pos, moore=True, include_center=False)
if neighbours:
same = sum(1for n in neighbours if n.agent_type == agent.agent_type)
fractions.append(same / len(neighbours))
returnfloat(np.mean(fractions)) if fractions else0.5defstep(self):
self.datacollector.collect(self)
self.schedule.step()
pct_happy = sum(a.is_happy for a inself.schedule.agents) / self.schedule.get_agent_count()
if pct_happy >= 1.0:
self.running = Falsedefrun_schelling(
n: int = 20,
tolerance: float = 0.3,
steps: int = 50,
density: float = 0.9,
seed: int = 42,
) -> pd.DataFrame:
"""
Run the Schelling segregation model and return step-level metrics.
Returns
-------
pd.DataFrame with columns: step, pct_happy, segregation_index
"""
model = SchellingModel(n=n, density=density, tolerance=tolerance, seed=seed)
for _ inrange(steps):
ifnot model.running:
break
model.step()
df = model.datacollector.get_model_vars_dataframe().reset_index()
df.columns = ["step", "pct_happy", "segregation_index"]
return df
5. Sample Entropy
defcompute_sample_entropy(ts: np.ndarray, m: int = 2, r: float = None) -> float:
"""
Compute sample entropy (SampEn) of a 1-D time series.
SampEn measures regularity: lower values indicate more self-similar series.
Parameters
----------
ts : np.ndarray — 1-D time series
m : int — template length (embedding dimension)
r : float or None — tolerance (default: 0.2 × std(ts))
Returns
-------
float — sample entropy (nan if undefined)
"""
ts = np.asarray(ts, dtype=float)
N = len(ts)
if r isNone:
r = 0.2 * np.std(ts, ddof=1)
def_count_matches(template_len: int) -> int:
count = 0for i inrange(N - template_len):
template = ts[i: i + template_len]
for j inrange(N - template_len):
if i == j:
continue
candidate = ts[j: j + template_len]
if np.max(np.abs(template - candidate)) <= r:
count += 1return count
B = _count_matches(m)
A = _count_matches(m + 1)
if B == 0or A == 0:
return np.nan
returnfloat(-np.log(A / B))
6. Bond Percolation Threshold
import networkx as nx
defpercolation_threshold_bond(
G: nx.Graph,
n_trials: int = 30,
p_values: np.ndarray = None,
) -> dict:
"""
Estimate the bond percolation threshold of a graph by measuring the relative
size of the giant component as a function of edge occupation probability p.
Parameters
----------
G : nx.Graph
n_trials : int — Monte Carlo repetitions per p value
p_values : np.ndarray or None — occupation probabilities (default: 0 to 1 in 20 steps)
Returns
-------
dict with keys: p_values, giant_fraction_mean, giant_fraction_std,
estimated_threshold
"""if p_values isNone:
p_values = np.linspace(0.0, 1.0, 21)
N = G.number_of_nodes()
edges = list(G.edges())
giant_means = []
giant_stds = []
for p in p_values:
fractions = []
for _ inrange(n_trials):
kept = [e for e in edges if np.random.rand() < p]
subgraph = nx.Graph()
subgraph.add_nodes_from(G.nodes())
subgraph.add_edges_from(kept)
components = sorted(nx.connected_components(subgraph), key=len, reverse=True)
if components:
fractions.append(len(components[0]) / N)
else:
fractions.append(0.0)
giant_means.append(np.mean(fractions))
giant_stds.append(np.std(fractions))
# Estimate threshold as the p where giant fraction crosses 0.5 × max
half_max = max(giant_means) / 2.0
threshold_idx = next(
(i for i, gf inenumerate(giant_means) if gf >= half_max),
len(p_values) - 1,
)
threshold = float(p_values[threshold_idx])
return {
"p_values": p_values,
"giant_fraction_mean": np.array(giant_means),
"giant_fraction_std": np.array(giant_stds),
"estimated_threshold": threshold,
}
7. Information-Theoretic Measures
from scipy.stats import entropy as scipy_entropy
from scipy.special import entr
defcompute_information_measures(x: np.ndarray, y: np.ndarray = None, bins: int = 20) -> dict:
"""
Compute Shannon entropy, mutual information, and normalised mutual information
for one or two continuous signals.
Parameters
----------
x : np.ndarray — first signal
y : np.ndarray or None — second signal (required for MI)
bins : int — number of histogram bins for density estimation
Returns
-------
dict with keys: entropy_x, entropy_y (if y given), mutual_info, nmi
"""
x = np.asarray(x, dtype=float)
px, _ = np.histogram(x, bins=bins, density=True)
px = px / px.sum()
px = px[px > 0]
Hx = float(scipy_entropy(px, base=2))
result = {"entropy_x": round(Hx, 4)}
if y isnotNone:
y = np.asarray(y, dtype=float)
py, _ = np.histogram(y, bins=bins, density=True)
py = py / py.sum()
py = py[py > 0]
Hy = float(scipy_entropy(py, base=2))
result["entropy_y"] = round(Hy, 4)
# Joint entropy
pxy, _, _ = np.histogram2d(x, y, bins=bins, density=True)
pxy_flat = pxy.ravel()
pxy_flat = pxy_flat / pxy_flat.sum()
pxy_flat = pxy_flat[pxy_flat > 0]
Hxy = float(scipy_entropy(pxy_flat, base=2))
result["joint_entropy"] = round(Hxy, 4)
MI = Hx + Hy - Hxy
NMI = 2 * MI / (Hx + Hy) if (Hx + Hy) > 0else0.0
result["mutual_info"] = round(MI, 4)
result["nmi"] = round(NMI, 4)
return result
8. Examples
Example A — Test if a Citation Network Degree Distribution Follows a Power Law
import networkx as nx
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# --- Generate a synthetic scale-free citation network (Barabasi-Albert) ---
np.random.seed(42)
G = nx.barabasi_albert_graph(n=2000, m=3, seed=42)
degrees = np.array([d for _, d in G.degree()])
print(f"Network: {G.number_of_nodes()} nodes, {G.number_of_edges()} edges")
print(f"Degree range: [{degrees.min()}, {degrees.max()}], mean={degrees.mean():.2f}")
# Fit power law
pl_result = fit_power_law(degrees, discrete=True)
print(f"\nPower-law fit:")
print(f" alpha = {pl_result['alpha']:.4f}")
print(f" xmin = {pl_result['xmin']:.0f}")
print(f" sigma = {pl_result['sigma']:.4f}")
print(f" n_tail = {pl_result['n_tail']} / {pl_result['n_total']}")
# Compare against alternatives
comparison = test_power_law_vs_alternatives(degrees, discrete=True)
print("\nDistribution comparison (vs power law):")
print(comparison.to_string(index=False))
# Visualise
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
plot_power_law(pl_result, title="Citation Network Degree Distribution", ax=axes[0])
# Percolation analysis on the same network
perc = percolation_threshold_bond(G, n_trials=20)
axes[1].plot(perc["p_values"], perc["giant_fraction_mean"], color="steelblue", marker="o")
axes[1].fill_between(
perc["p_values"],
perc["giant_fraction_mean"] - perc["giant_fraction_std"],
perc["giant_fraction_mean"] + perc["giant_fraction_std"],
alpha=0.3, color="steelblue",
)
axes[1].axvline(perc["estimated_threshold"], linestyle="--", color="red",
label=f"p_c ≈ {perc['estimated_threshold']:.2f}")
axes[1].set_xlabel("Bond occupation probability p")
axes[1].set_ylabel("Relative size of giant component")
axes[1].set_title("Bond Percolation on BA Network")
axes[1].legend()
plt.tight_layout()
plt.savefig("/tmp/citation_network_analysis.png", dpi=150)
plt.show()
print(f"\nEstimated percolation threshold: p_c ≈ {perc['estimated_threshold']:.3f}")
Example B — Schelling Segregation Simulation with Varying Tolerance
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
TOLERANCES = [0.2, 0.3, 0.4, 0.5, 0.6, 0.7]
N_GRID = 25
STEPS = 80
results = []
for tol in TOLERANCES:
df_run = run_schelling(n=N_GRID, tolerance=tol, steps=STEPS, seed=0)
df_run["tolerance"] = tol
results.append(df_run)
df_all = pd.concat(results, ignore_index=True)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# Plot 1: % happy over time for each tolerancefor tol, grp in df_all.groupby("tolerance"):
axes[0].plot(grp["step"], grp["pct_happy"], label=f"τ={tol}")
axes[0].set_xlabel("Step")
axes[0].set_ylabel("Fraction of Happy Agents")
axes[0].set_title("Schelling Model: Happiness vs Time")
axes[0].legend(title="Tolerance", fontsize=8)
axes[0].set_ylim(0, 1.05)
# Plot 2: final segregation index vs tolerance
final_seg = (
df_all.groupby("tolerance")
.apply(lambda g: g.sort_values("step").iloc[-1]["segregation_index"])
.reset_index()
)
final_seg.columns = ["tolerance", "final_segregation"]
axes[1].plot(final_seg["tolerance"], final_seg["final_segregation"], marker="o", color="tomato")
axes[1].set_xlabel("Tolerance Threshold τ")
axes[1].set_ylabel("Final Segregation Index")
axes[1].set_title("Schelling Model: Segregation vs Tolerance")
axes[1].set_ylim(0, 1.0)
plt.tight_layout()
plt.savefig("/tmp/schelling_analysis.png", dpi=150)
plt.show()
# Print summary table
summary = df_all.groupby("tolerance").agg(
final_happy=("pct_happy", "last"),
final_segregation=("segregation_index", "last"),
steps_to_converge=("step", "max"),
).reset_index()
print("\nSchelling Model Summary:")
print(summary.to_string(index=False))
# Hurst exponent on the happiness time-series for τ=0.5
ts_happy = df_all[df_all["tolerance"] == 0.5]["pct_happy"].values
hurst_result = compute_hurst(ts_happy)
print(f"\nHurst exponent of happiness time-series (τ=0.5): H = {hurst_result['hurst']:.4f}")
print(f"R² of log-log fit: {hurst_result['r_squared']:.4f}")
# Sample entropy of the same series
se = compute_sample_entropy(ts_happy, m=2)
print(f"Sample entropy: {se:.4f}")
# Information measures comparing tolerant vs strict agents
ts_tol = df_all[df_all["tolerance"] == 0.7]["pct_happy"].values
ts_strict = df_all[df_all["tolerance"] == 0.2]["pct_happy"].values
min_len = min(len(ts_tol), len(ts_strict))
info = compute_information_measures(ts_tol[:min_len], ts_strict[:min_len])
print(f"\nInformation measures (tolerant vs strict happiness series):")
print(f" H(tolerant) = {info['entropy_x']:.4f} bits")
print(f" H(strict) = {info['entropy_y']:.4f} bits")
print(f" MI = {info['mutual_info']:.4f} bits")
print(f" NMI = {info['nmi']:.4f}")
9. Tips and Gotchas
Power-law vs lognormal: Real-world data rarely follows a pure power law. The
Vuong test frequently returns "inconclusive". Report both the fit quality and the
comparison result; never claim power law on visual inspection alone.
Hurst exponent bias: Short series (N < 500) give biased H estimates. Use
detrended fluctuation analysis (DFA) for non-stationary series; the nolds package
provides nolds.hurst_rs and nolds.dfa.
Schelling grid size: A 20×20 grid with density=0.9 has ~360 agents; results are
stochastic. Average over ≥10 seeds before drawing conclusions about the tolerance
threshold.
Sample entropy runtime: The naive O(N²) implementation above is slow for N > 1000.
Use antropy.sample_entropy (C extension) for large time series.
Mesa 2.x API: RandomActivation and SingleGrid are in mesa.time and
mesa.space respectively. Mesa 3.x reorganised these; pin mesa>=2.1,<3 or
update the imports accordingly.
Percolation on directed graphs: nx.connected_components only works on undirected
graphs. For directed networks use nx.strongly_connected_components.
