Misspecified N parameter in Poisson process generator

[joaohbbittar]

N parameter (poisson_batch_size): incorrect documentation and undocumented Gamma process behaviour

Problem

The parameter N (exposed as poisson_batch_size in DescendingDrive__Pool) is incorrectly
documented in all three places where it appears. More importantly, there is a behavioural issue
that the documentation silently hides: N > 1 does not produce a Poisson process — it produces
a Gamma process.

Status: discussion only. This PR contains no code changes and no new files committed to the
repo. It exists to flag the issue, share the empirical evidence below, and get agreement from
maintainers on which resolution to take — and where a validation notebook like the one below
should live — before anyone implements it.

Empirical findings

This claim was tested in a standalone notebook,
test_poisson_batch_size_synchrony_and_gamma.ipynb (not committed here — see "Test notebook"
below for why, and how to get it). For N ∈ {1, 2, 5, 10, 20}, pooled ISI statistics from
simulated DescendingDrive__Pool populations match the Gamma(N, N) prediction closely: CV(ISI)
tracks 1/√N, and a KS test against Gamma(N, N) is not rejected while a KS test against
Exponential (the true-Poisson null) is strongly rejected for every N > 1. Only N = 1 passes
as Poisson. This confirms the derivation below empirically rather than just analytically.

How the parameter actually works

The Cython implementation generates each inter-spike threshold via the log-of-product method:

aux = 1.0
for _ in range(N):
    aux *= rand_uniform()
thres = -(1.0 / N) * log(aux)

Because each -log(Uᵢ) ~ Exp(1), the threshold is the average of N independent Exp(1)
variates, which follows an Erlang(N, N) = Gamma(shape=N, rate=N) distribution.

ISI statistics as a function of N

N	Threshold distribution	CV of ISI (predicted)
1	Exp(1) = Gamma(1, 1)	1.0 — true Poisson
2	Gamma(2, 2)	0.71
5	Gamma(5, 5)	0.45
N	Gamma(N, N)	1/√N
∞	δ(1)	0 — clock

The test notebook (see "Test notebook" below) confirms this prediction empirically for
N ∈ {1, 2, 5, 10, 20} (empirical CV(ISI) tracks 1/√N closely, KS test against Gamma(N, N)
not rejected, KS test against Exponential strongly rejected for N > 1) — exact figures,
statistics, and histogram/fit overlays are in the notebook rather than reproduced here.

The mean firing rate is preserved for all N, but ISI variance shrinks as 1/N. Increasing N
does not improve a Poisson approximation — it progressively replaces it with a Gamma process
trending toward perfectly regular firing.

Conflict with the existing Gamma pathway

The codebase already provides DD_Gamma / process_type="gamma" for Gamma-distributed ISIs.
DD with N > 1 implements the same process, silently and with an incompatible
parameterisation (here rate = shape = N always; DD_Gamma allows them to vary independently).
A caller setting poisson_batch_size=16 while using process_type="poisson" is unknowingly
requesting Gamma(16, 16) ISIs with CV ≈ 0.25.

Note — a possible (but empirically unsupported) motivation for N > 1.
One charitable reading of why N > 1 exists at all is that it was meant not to shape the ISI
distribution, but to advance the RNG state faster between spikes. In a typical network
simulation, many DD cells run in parallel, each seeded from a nearby integer (via
derive_subseed(class__ID, global__ID)). The xorshift64* generator used here has good global
statistical properties, but like all linear feedback shift register designs its output has
short-range structure — consecutive draws are not fully independent. When many generators are
seeded close together, their streams could in principle be correlated at short lags, introducing
spurious synchrony across the population. Consuming N draws per threshold does advance the RNG
state N steps, which would help with this — but the current formula can't separate that effect
from the distributional change, since folding all N draws into the threshold via
-(1/N) * log(aux) is exactly what turns Exp(1) into Gamma(N, N).The current seeding mechanism for neurons is probably already avoiding this issue.

This was tested directly and the artefact was not observed. For N = 1 (the setting the
docs currently call "true Poisson"), a population of 500 parallel DescendingDrive__Pool
neurons seeded from nearby integers was checked for excess synchrony against a jittered negative
control (population-averaged cross-correlogram and pairwise spike-count correlations). No
significant zero-lag CCG peak or positive shift in pairwise correlations was detected relative
to the jittered control at this population size and seed spacing.

So: plausible original intent, real theoretical mechanism, but not currently supported by
evidence at the scale tested. It's noted here only because it may explain how the current code
came to be, and because it could conceivably still matter at different population sizes,
seed-spacing schemes, or RNGs. Whether RNG decorrelation was the actual original intent is
unknown either way — there's no comment in the code to that effect — and this note should not
be read as a confirmed issue the way the Gamma-process bug above is.

---

Documentation errors

Three files need correction regardless of which behavioural fix is chosen.

cells.py — DD docstring:
Current description calls N a "maximum firing rate in Hz", which has no connection to the
implementation. N does not affect rate at all.

descending_drive.py — poisson_batch_size:
Current description says higher values "improve statistical accuracy". This is wrong: N=1 is
already an exact draw. Larger N changes the distribution, not the accuracy of sampling from it.

_poisson_process_generator.pyx:
Current description says higher values "may improve statistical properties" — undefined and
misleading for the same reason.

---

Possible resolutions — decision needed

Both options below require existing callers using N > 1 under process_type="poisson" to
change their call — neither is a free lunch on that front. What differs is where the surviving
capability lives, and how big the migration is.

Option A — Reject N > 1 for process_type="poisson"
Either raise a ValueError, or emit a UserWarning (and proceed with the existing, now-honestly-
documented Gamma(N, N) behaviour), when N > 1 is passed alongside process_type="poisson".
Correct all docstrings to describe the actual Gamma(N, N) behaviour. A caller hitting this has a
one-argument fix: switch to DD(process_type="gamma", N=16) (or, if warning rather than raising,
no fix is strictly required, but the warning makes clear they're getting Gamma ISIs rather than
Poisson). The error variant is the stricter, more honest option; the warning variant is softer and
non-breaking, but lets a caller keep silently relying on Gamma-under-a-Poisson-label if they
ignore the warning — which is close to the status quo this PR is trying to fix.

Option B — Deprecate N in DD entirely
Route all Gamma-like behaviour through DD_Gamma / process_type="gamma", which already
exists and is more flexible. DD would only ever generate a true Poisson process. A caller
hitting this has to switch classes: DD_Gamma(shape=16, rate=16) instead of
DD(process_type="poisson", N=16). Cleaner long-term API (no overloaded, partially-ignored N
parameter on DD), but a larger migration than Option A's same-class argument swap.

Optional addition, only if the decorrelation note above turns out to matter: a separate
rng_burn parameter could be introduced alongside whichever option is chosen above, defaulting
to 0 (non-breaking). It would discard extra RNG calls after an exact Exp(1) draw, advancing
the RNG state without touching the ISI distribution:

thres = -log(rand_uniform())
for _ in range(rng_burn):
    rand_uniform()  # advance RNG state, value discarded

This keeps CV = 1 (genuine Poisson) while giving independent control over how many RNG steps
are consumed per spike. Since the empirical findings above don't currently motivate it, this can
simply be skipped unless maintainers know of a regime where it's needed.

---

Files affected

File	Change
myogen/simulator/neuron/_cython/_poisson_process_generator.pyx	Correct N in Parameters and Attributes
myogen/simulator/neuron/cells.py	Correct N; remove false "maximum firing rate" claim
myogen/simulator/neuron/populations/descending_drive.py	Correct poisson_batch_size; add guard or deprecation per chosen option

None of the rows above are implemented in this PR; they're listed so a follow-up implementation
PR has a checklist once a resolution is agreed. The validation notebook used to produce the
"Empirical findings" above is not in this table since it isn't being committed — see "Test
notebook" below.

---

Test notebook

The numbers cited in "Empirical findings" above came from a notebook,
test_poisson_batch_size_synchrony_and_gamma.ipynb, that I wrote to check the Gamma-process
claim (and the decorrelation note) rather than just asserting them. It's not committed to this
PR:

• There's no established convention in this repo for where a validation/evidence notebook like
  this should live (there's no tests/notebooks/ or similar today), and
  picking a location unilaterally felt like the wrong call for a discussion-only PR. Happy to add
  it properly once maintainers have a preference — or to help establish that convention if there
  isn't one yet.
• It's exploratory/evidence-gathering, not a permanent regression test, so it may not belong in
  the repo at all depending on how you'd want to treat this kind of artifact long-term.

It's available on request — I'm happy to attach it directly to this PR/issue (e.g. as a
comment attachment) or share it however's easiest, so you can re-run or inspect it without
having to take the empirical claims above on faith. It runs real DescendingDrive__Pool
simulations via NEURON at a constant firing rate (isolating the generator's own statistics from
any drive-profile effects) and is already executed, with plots and summary tables, so it can be
read without re-running.

If maintainers decide on Option A or B above, a small permanent regression test (e.g. a
pytest-based KS-test assertion that N = 1 ISIs are Exponential) would likely be worth adding
separately in an implementation PR — happy to extract one from this notebook if useful.

[RaulSimpetru]

Thanks @joaohbbittar — your analysis nailed it. Confirmed empirically: at N = 5 the "Poisson" generator produced CV(ISI) ≈ 0.447 = 1/√5, matching Gamma(5, 5) exactly.

The fix is implemented and folded into #19 (commit 603b245):

• The generator now draws a single Exp(1) inter-spike threshold — a genuine Poisson process (CV = 1). poisson_batch_size is removed; regular/Gamma firing goes through process_type="gamma", shape=N.
• The three docstrings you flagged are corrected. In addition, the afferents' poisson_batch_size — which was really the Gamma shape of the AffIa/AffII/AffIb generators — is renamed to shape.
• Per your suggestion, this ships a permanent regression test rather than a throwaway notebook: tests/test_poisson_process.py proves the process is Poisson via statistical goodness-of-fit (KS vs Exponential, χ² on spike counts, Fano factor, serial independence), through both the public API and the raw generator. The Gamma path is validated symmetrically (KS fit + shape recovery + wrong-shape rejection).

Closing as superseded by #19. Thanks again for the careful writeup and the empirical validation — it made this a quick, confident fix.