Gaussian process Poisson regression

Poisson regression is suitable for estimating the latent rate of count data, such as the number of buses arriving at a stop within an hour or the number of car accidents per day. Here, we estimate the rate using a Gaussian process given synthetic count data. First, we draw the log rate \(\eta\) from a Gaussian process with squared exponential kernel \(k\) observed at discrete positions \(x\), i.e. the covariance of two observations is

\[ \mathrm{cov}\left(\eta(x),\eta(y)\right)=\sigma^2\exp\left(-\frac{d(x, y)^2}{2\ell^2}\right) + \epsilon\delta(x - y), \]

where \(d\left(\cdot,\cdot\right)\) is a distance measure, \(\sigma\) is the marginal scale of the process, and \(\ell\) is its correlation length. The “nugget” variance \(\epsilon\) ensures the covariance matrix is numerically positive semi-definite. Second, we sample counts \(y\) from a Poisson distribution with rate \(\lambda = \exp\eta\). Let’s generate a sample with \(n\) regularly-spaced observations and visualize it.

from gptools.util.kernels import DiagonalKernel, ExpQuadKernel, Kernel
import matplotlib as mpl
from matplotlib import pyplot as plt
import numpy as np


def simulate(x: np.ndarray, kernel: Kernel, mu: float = 0) -> dict:
    """Simulate a Gaussian process Poisson model with mean mu observed at x."""
    # Add an extra dimension to the vector x because the kernel expects an array
    # with shape (num_data_points, num_dimensions).
    X = x[:, None]
    cov = kernel.evaluate(X)
    eta = np.random.multivariate_normal(mu * np.ones_like(x), cov)
    rate = np.exp(eta)
    y = np.random.poisson(rate)
    # Return the results as a dictionary, including input arguments (we need to extract kernel
    # parameters from the `CompositeKernel`).
    return {"x": x, "X": X, "mu": mu, "y": y, "eta": eta, "y": y, "rate": rate, "n": x.size,
            "sigma": kernel.a.sigma, "length_scale": kernel.a.length_scale,
            "epsilon": kernel.b.epsilon}


def plot_sample(sample: dict, ax: mpl.axes.Axes = None) -> mpl.axes.Axes:
    """Visualize a sample of a Gaussian process Poisson model."""
    ax = ax or plt.gca()
    ax.plot(sample["x"], sample["rate"], label=r"rate $\lambda$")
    ax.scatter(sample["x"], sample["y"], marker=".", color="k", label="counts $y$",
               alpha=0.5)
    ax.set_xlabel("Location $x$")
    return ax


np.random.seed(0)
n = 64
x = np.arange(n)
kernel = ExpQuadKernel(sigma=1.2, length_scale=5, period=n) + DiagonalKernel(1e-3, n)
sample = simulate(x, kernel)
plot_sample(sample).legend()

<matplotlib.legend.Legend at 0x710acc261330>

We used a kernel with periodic boundary conditions by passing period=n to the kernel. The Gaussian process is thus defined on a ring with circumference \(n\).

To learn the latent rate using stan, we first define the data block of the program, i.e., the information we provide to the inference algorithm. Here, we assume that the marginal scale \(\sigma\), correlation length \(\ell\), and nugget variance \(\epsilon\) are known.

// Shared data definition for Poisson regression models.

// Kernel parameters.
real<lower=0> sigma, length_scale, epsilon;

// Information about observation locations and counts.
int n;
array [n] int y;
array [n] vector[1] X;

The model is simple: a multivariate normal prior for the log rate \(\eta\) and a Poisson observation model.

// Gaussian process with log link for Poisson observations using the default centered
// parameterization.

data {
    // Shared data definitions.
    #include data.stan
}

parameters {
    // Latent log-rate we seek to infer.
    vector[n] eta;
}

model {
    // Gaussian process prior and observation model. The covariance should technically have periodic
    // boundary conditions to match the generative model in the example.
    matrix[n, n] cov = add_diag(gp_exp_quad_cov(X, X, sigma, length_scale), epsilon);
    eta ~ multi_normal(zeros_vector(n), cov);
    y ~ poisson_log(eta);
}

Let’s compile the model, draw posterior samples, and visualize the result.

from cmdstanpy import CmdStanMCMC
from gptools.stan import compile_model
from gptools.util import Timer
from gptools.util.plotting import plot_band
import os


def sample_and_plot(stan_file: str, data: dict, return_fit: bool = False, **kwargs) -> CmdStanMCMC:
    """Draw samples from the posterior and visualize them."""
    # Set default parameters. We use a small number of samples during testing.
    niter = 10 if os.environ.get("CI") else 100
    defaults = {"iter_warmup": niter, "iter_sampling": niter, "chains": 1,
                "refresh": niter // 10 or None}
    defaults.update(kwargs)
    kwargs = defaults

    # Compile the model and draw posterior samples.
    model = compile_model(stan_file=stan_file)
    with Timer(f"sampled using {stan_file}"):
        fit = model.sample(data, **kwargs)

    # Visualize the result.
    ax = plot_sample(sample)
    plot_band(x, np.exp(fit.stan_variable("eta")), label="inferred rate", ax=ax)
    ax.legend()
    if return_fit:
        return fit


centered_fit = sample_and_plot("poisson_regression_centered.stan", sample, return_fit=True)

WARNING:root:cmdstan<2.36.0 had a bug in the evaluation of Matern 3/2 kernels with different length scales for different dimensions; see https://github.com/stan-dev/math/pull/3084 for details. Your model may yield unexpected results or crash if you use nearest-neighbor Gaussian processes with Matern 3/2 kernels. Your cmdstan version is 2.31; consider upgrading.

WARNING:cmdstanpy:CmdStanModel(compile=...) is deprecated and will be removed in the next major version. The constructor will always ensure a model has a compiled executable.
If you wish to force recompilation, use force_compile=True instead.

INFO:cmdstanpy:compiling stan file /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_centered.stan to exe file /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_centered

DEBUG:cmdstanpy:cmd: make STANCFLAGS+=--filename-in-msg=poisson_regression_centered.stan STANCFLAGS+=--include-paths=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/python/stan/gptools/stan,/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_centered
cwd: /home/docs/.cmdstan/cmdstan-2.31.0

DEBUG:cmdstanpy:Console output:

--- Translating Stan model to C++ code ---
bin/stanc --filename-in-msg=poisson_regression_centered.stan --include-paths=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/python/stan/gptools/stan,/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression --o=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_centered.hpp /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_centered.stan

--- Compiling, linking C++ code ---
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.9 -I stan/lib/stan_math/lib/boost_1.78.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS          -c -Wno-ignored-attributes   -x c++ -o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_centered.o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_centered.hpp
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.9 -I stan/lib/stan_math/lib/boost_1.78.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS                -Wl,-L,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb"      /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_centered.o src/cmdstan/main.o        -Wl,-L,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb"   stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_nvecserial.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_cvodes.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_idas.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_kinsol.a  stan/lib/stan_math/lib/tbb/libtbb.so.2 -o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_centered
rm -f /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_centered.o

INFO:cmdstanpy:compiled model executable: /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_centered

DEBUG:cmdstanpy:cmd: /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_centered info
cwd: None

DEBUG:cmdstanpy:input tempfile: /tmp/tmpikaxryv3/wc_a1btt.json

INFO:cmdstanpy:CmdStan start processing

DEBUG:cmdstanpy:idx 0

DEBUG:cmdstanpy:running CmdStan, num_threads: 1

DEBUG:cmdstanpy:CmdStan args: ['/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_centered', 'id=1', 'random', 'seed=59449', 'data', 'file=/tmp/tmpikaxryv3/wc_a1btt.json', 'output', 'file=/tmp/tmpikaxryv3/poisson_regression_centered0nrsbnc_/poisson_regression_centered-20250317164850.csv', 'refresh=10', 'method=sample', 'num_samples=100', 'num_warmup=100', 'algorithm=hmc', 'adapt', 'engaged=1']

INFO:cmdstanpy:CmdStan done processing.

DEBUG:cmdstanpy:runset
RunSet: chains=1, chain_ids=[1], num_processes=1
 cmd (chain 1):
	['/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_centered', 'id=1', 'random', 'seed=59449', 'data', 'file=/tmp/tmpikaxryv3/wc_a1btt.json', 'output', 'file=/tmp/tmpikaxryv3/poisson_regression_centered0nrsbnc_/poisson_regression_centered-20250317164850.csv', 'refresh=10', 'method=sample', 'num_samples=100', 'num_warmup=100', 'algorithm=hmc', 'adapt', 'engaged=1']
 retcodes=[0]
 per-chain output files (showing chain 1 only):
 csv_file:
	/tmp/tmpikaxryv3/poisson_regression_centered0nrsbnc_/poisson_regression_centered-20250317164850.csv
 console_msgs (if any):
	/tmp/tmpikaxryv3/poisson_regression_centered0nrsbnc_/poisson_regression_centered-20250317164850_0-stdout.txt

DEBUG:cmdstanpy:Chain 1 console:
method = sample (Default)
  sample
    num_samples = 100
    num_warmup = 100
    save_warmup = 0 (Default)
    thin = 1 (Default)
    adapt
      engaged = 1 (Default)
      gamma = 0.050000000000000003 (Default)
      delta = 0.80000000000000004 (Default)
      kappa = 0.75 (Default)
      t0 = 10 (Default)
      init_buffer = 75 (Default)
      term_buffer = 50 (Default)
      window = 25 (Default)
    algorithm = hmc (Default)
      hmc
        engine = nuts (Default)
          nuts
            max_depth = 10 (Default)
        metric = diag_e (Default)
        metric_file =  (Default)
        stepsize = 1 (Default)
        stepsize_jitter = 0 (Default)
    num_chains = 1 (Default)
id = 1 (Default)
data
  file = /tmp/tmpikaxryv3/wc_a1btt.json
init = 2 (Default)
random
  seed = 59449
output
  file = /tmp/tmpikaxryv3/poisson_regression_centered0nrsbnc_/poisson_regression_centered-20250317164850.csv
  diagnostic_file =  (Default)
  refresh = 10
  sig_figs = -1 (Default)
  profile_file = profile.csv (Default)
num_threads = 1 (Default)


Gradient evaluation took 0.000318 seconds
1000 transitions using 10 leapfrog steps per transition would take 3.18 seconds.
Adjust your expectations accordingly!


WARNING: There aren't enough warmup iterations to fit the
         three stages of adaptation as currently configured.
         Reducing each adaptation stage to 15%/75%/10% of
         the given number of warmup iterations:
           init_buffer = 15
           adapt_window = 75
           term_buffer = 10

Iteration:   1 / 200 [  0%]  (Warmup)
Iteration:  10 / 200 [  5%]  (Warmup)
Iteration:  20 / 200 [ 10%]  (Warmup)
Iteration:  30 / 200 [ 15%]  (Warmup)
Iteration:  40 / 200 [ 20%]  (Warmup)
Iteration:  50 / 200 [ 25%]  (Warmup)
Iteration:  60 / 200 [ 30%]  (Warmup)
Iteration:  70 / 200 [ 35%]  (Warmup)
Iteration:  80 / 200 [ 40%]  (Warmup)
Iteration:  90 / 200 [ 45%]  (Warmup)
Iteration: 100 / 200 [ 50%]  (Warmup)
Iteration: 101 / 200 [ 50%]  (Sampling)
Iteration: 110 / 200 [ 55%]  (Sampling)
Iteration: 120 / 200 [ 60%]  (Sampling)
Iteration: 130 / 200 [ 65%]  (Sampling)
Iteration: 140 / 200 [ 70%]  (Sampling)
Iteration: 150 / 200 [ 75%]  (Sampling)
Iteration: 160 / 200 [ 80%]  (Sampling)
Iteration: 170 / 200 [ 85%]  (Sampling)
Iteration: 180 / 200 [ 90%]  (Sampling)
Iteration: 190 / 200 [ 95%]  (Sampling)
Iteration: 200 / 200 [100%]  (Sampling)

 Elapsed Time: 4.098 seconds (Warm-up)
               4.185 seconds (Sampling)
               8.283 seconds (Total)

sampled using poisson_regression_centered.stan in 8.319 seconds

Stan’s Hamiltonian Monte Carlo sampler explores the posterior distribution well and recovers the underlying rate. But it takes its time about it–especially given how small the dataset is. The exploration is slow because adjacenct elements of the log rate are highly correlated under the posterior distribution as shown in the scatter plot below. The mean correlation between adjacent elements is 0.95. The posterior is highly correlated because the Gaussian process prior demands that adjacent points are close: if one moves up, the other must follow to keep the log rate \(\eta\) smooth. If the data are “strong”, i.e., there is more information in the likelihood than in the prior, the data can decorrelate the posterior: each log rate is primarily informed by the corresponding count rather than the prior.

fig, ax = plt.subplots()
etas = centered_fit.stan_variable("eta")
idx = 20
ax.scatter(etas[:, idx - 1], etas[:, idx]).set_edgecolor("w")
ax.set_aspect("equal")
ax.set_xlabel(fr"Log rate $\eta_{{{idx}}}$")
ax.set_ylabel(fr"Log rate $\eta_{{{idx + 1}}}$")
fig.tight_layout()

In this example, the data are relatively weak (recall that the variance of the Poisson distribution is equal to its mean such that, for small counts, the likelihood is not particularly informative). We can alleviate this problem by switching to a non-centered parameterization such that the parameters are uncorrelated under the prior distribution. In particular, we use a parameter \(z\sim\text{Normal}\left(0, 1\right)\) and transform the white noise such that it is a draw from the Gaussian process prior. This approach is a higher-dimensional analogue of sampling a standard normal distribution and multiplying by the standard deviation \(\sigma\) as opposed to sampling from a normal distribution with standard deviation \(\sigma\) directly. The log rate is

\[ \eta = Lz, \]

where \(L\) is the Cholesky decomposition of the covariance matrix \(\Sigma\) of the kernel. The Cholesky decomposition of the covariance matrix is the higher-dimensional analogue of taking the square route of the variance. Because the data are weak, they can only slightly correlate the parameters \(z\). The Stan program is shown below.

// Gaussian process with log link for Poisson observations.

data {
    #include data.stan
}

parameters {
    vector[n] z;
}

transformed parameters {
    vector[n] eta;
    {
        // Evaluate the covariance and apply the Cholesky decomposition to the white noise z. We
        // wrap the evaluation in braces because Stan only writes top-level variables to the output
        // CSV files, and we don't need to store the entire covariance matrix. The covariance should
        // technically have periodic boundary conditions to match the generative model in the
        // example.
        matrix[n, n] cov = add_diag(gp_exp_quad_cov(X, X, sigma, length_scale), epsilon);
        eta = cholesky_decompose(cov) * z;
    }
}

model {
    // White noise prior (implies eta ~ multi_normal(zeros_vector(n), cov)) and observation model.
    z ~ normal(0, 1);
    y ~ poisson_log(eta);
}

non_centered_fit = sample_and_plot("poisson_regression_non_centered.stan", sample, return_fit=True)

WARNING:root:cmdstan<2.36.0 had a bug in the evaluation of Matern 3/2 kernels with different length scales for different dimensions; see https://github.com/stan-dev/math/pull/3084 for details. Your model may yield unexpected results or crash if you use nearest-neighbor Gaussian processes with Matern 3/2 kernels. Your cmdstan version is 2.31; consider upgrading.

WARNING:cmdstanpy:CmdStanModel(compile=...) is deprecated and will be removed in the next major version. The constructor will always ensure a model has a compiled executable.
If you wish to force recompilation, use force_compile=True instead.

INFO:cmdstanpy:compiling stan file /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_non_centered.stan to exe file /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_non_centered

DEBUG:cmdstanpy:cmd: make STANCFLAGS+=--filename-in-msg=poisson_regression_non_centered.stan STANCFLAGS+=--include-paths=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/python/stan/gptools/stan,/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_non_centered
cwd: /home/docs/.cmdstan/cmdstan-2.31.0

DEBUG:cmdstanpy:Console output:

--- Translating Stan model to C++ code ---
bin/stanc --filename-in-msg=poisson_regression_non_centered.stan --include-paths=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/python/stan/gptools/stan,/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression --o=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_non_centered.hpp /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_non_centered.stan

--- Compiling, linking C++ code ---
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.9 -I stan/lib/stan_math/lib/boost_1.78.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS          -c -Wno-ignored-attributes   -x c++ -o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_non_centered.o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_non_centered.hpp
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.9 -I stan/lib/stan_math/lib/boost_1.78.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS                -Wl,-L,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb"      /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_non_centered.o src/cmdstan/main.o        -Wl,-L,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb"   stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_nvecserial.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_cvodes.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_idas.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_kinsol.a  stan/lib/stan_math/lib/tbb/libtbb.so.2 -o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_non_centered
rm -f /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_non_centered.o

INFO:cmdstanpy:compiled model executable: /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_non_centered

DEBUG:cmdstanpy:cmd: /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_non_centered info
cwd: None

DEBUG:cmdstanpy:input tempfile: /tmp/tmpikaxryv3/7_cy9xth.json

INFO:cmdstanpy:CmdStan start processing

DEBUG:cmdstanpy:idx 0

DEBUG:cmdstanpy:running CmdStan, num_threads: 1

DEBUG:cmdstanpy:CmdStan args: ['/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_non_centered', 'id=1', 'random', 'seed=93567', 'data', 'file=/tmp/tmpikaxryv3/7_cy9xth.json', 'output', 'file=/tmp/tmpikaxryv3/poisson_regression_non_centeredcvhfdl02/poisson_regression_non_centered-20250317164921.csv', 'refresh=10', 'method=sample', 'num_samples=100', 'num_warmup=100', 'algorithm=hmc', 'adapt', 'engaged=1']

INFO:cmdstanpy:CmdStan done processing.

DEBUG:cmdstanpy:runset
RunSet: chains=1, chain_ids=[1], num_processes=1
 cmd (chain 1):
	['/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_non_centered', 'id=1', 'random', 'seed=93567', 'data', 'file=/tmp/tmpikaxryv3/7_cy9xth.json', 'output', 'file=/tmp/tmpikaxryv3/poisson_regression_non_centeredcvhfdl02/poisson_regression_non_centered-20250317164921.csv', 'refresh=10', 'method=sample', 'num_samples=100', 'num_warmup=100', 'algorithm=hmc', 'adapt', 'engaged=1']
 retcodes=[0]
 per-chain output files (showing chain 1 only):
 csv_file:
	/tmp/tmpikaxryv3/poisson_regression_non_centeredcvhfdl02/poisson_regression_non_centered-20250317164921.csv
 console_msgs (if any):
	/tmp/tmpikaxryv3/poisson_regression_non_centeredcvhfdl02/poisson_regression_non_centered-20250317164921_0-stdout.txt

DEBUG:cmdstanpy:Chain 1 console:
method = sample (Default)
  sample
    num_samples = 100
    num_warmup = 100
    save_warmup = 0 (Default)
    thin = 1 (Default)
    adapt
      engaged = 1 (Default)
      gamma = 0.050000000000000003 (Default)
      delta = 0.80000000000000004 (Default)
      kappa = 0.75 (Default)
      t0 = 10 (Default)
      init_buffer = 75 (Default)
      term_buffer = 50 (Default)
      window = 25 (Default)
    algorithm = hmc (Default)
      hmc
        engine = nuts (Default)
          nuts
            max_depth = 10 (Default)
        metric = diag_e (Default)
        metric_file =  (Default)
        stepsize = 1 (Default)
        stepsize_jitter = 0 (Default)
    num_chains = 1 (Default)
id = 1 (Default)
data
  file = /tmp/tmpikaxryv3/7_cy9xth.json
init = 2 (Default)
random
  seed = 93567
output
  file = /tmp/tmpikaxryv3/poisson_regression_non_centeredcvhfdl02/poisson_regression_non_centered-20250317164921.csv
  diagnostic_file =  (Default)
  refresh = 10
  sig_figs = -1 (Default)
  profile_file = profile.csv (Default)
num_threads = 1 (Default)


Gradient evaluation took 0.000311 seconds
1000 transitions using 10 leapfrog steps per transition would take 3.11 seconds.
Adjust your expectations accordingly!


WARNING: There aren't enough warmup iterations to fit the
         three stages of adaptation as currently configured.
         Reducing each adaptation stage to 15%/75%/10% of
         the given number of warmup iterations:
           init_buffer = 15
           adapt_window = 75
           term_buffer = 10

Iteration:   1 / 200 [  0%]  (Warmup)
Iteration:  10 / 200 [  5%]  (Warmup)
Iteration:  20 / 200 [ 10%]  (Warmup)
Iteration:  30 / 200 [ 15%]  (Warmup)
Iteration:  40 / 200 [ 20%]  (Warmup)
Iteration:  50 / 200 [ 25%]  (Warmup)
Iteration:  60 / 200 [ 30%]  (Warmup)
Iteration:  70 / 200 [ 35%]  (Warmup)
Iteration:  80 / 200 [ 40%]  (Warmup)
Iteration:  90 / 200 [ 45%]  (Warmup)
Iteration: 100 / 200 [ 50%]  (Warmup)
Iteration: 101 / 200 [ 50%]  (Sampling)
Iteration: 110 / 200 [ 55%]  (Sampling)
Iteration: 120 / 200 [ 60%]  (Sampling)
Iteration: 130 / 200 [ 65%]  (Sampling)
Iteration: 140 / 200 [ 70%]  (Sampling)
Iteration: 150 / 200 [ 75%]  (Sampling)
Iteration: 160 / 200 [ 80%]  (Sampling)
Iteration: 170 / 200 [ 85%]  (Sampling)
Iteration: 180 / 200 [ 90%]  (Sampling)
Iteration: 190 / 200 [ 95%]  (Sampling)
Iteration: 200 / 200 [100%]  (Sampling)

 Elapsed Time: 0.438 seconds (Warm-up)
               0.437 seconds (Sampling)
               0.875 seconds (Total)

sampled using poisson_regression_non_centered.stan in 0.907 seconds

Sampling from the posterior is substantially faster. While the log rates \(\eta\) remain highly correlated (see left panel below), the posterior for the parameters \(z\) the Hamiltonian sampler is exploring are virtually uncorrelated (see right panel below). Takeaway: if the data are strong, use a centered parameterization. If the data are weak, use a non-centered parameterization.

fig, (ax1, ax2) = plt.subplots(1, 2)
etas = non_centered_fit.stan_variable("eta")
zs = non_centered_fit.stan_variable("z")
idx = 20
ax1.scatter(etas[:, idx - 1], etas[:, idx]).set_edgecolor("w")
ax1.set_aspect("equal")
ax1.set_xlabel(fr"Log rate $\eta_{{{idx}}}$")
ax1.set_ylabel(fr"Log rate $\eta_{{{idx + 1}}}$")

ax2.scatter(zs[:, idx - 1], zs[:, idx]).set_edgecolor("w")
ax2.set_aspect("equal")
ax2.set_xlabel(fr"Log rate $z_{{{idx}}}$")
ax2.set_ylabel(fr"Log rate $z_{{{idx + 1}}}$")

fig.tight_layout()

Nearest-neighbor Gaussian processes

Despite the non-centered parameterization being faster, its performance scales poorly with increasing sample size. Evaluating the Cholesky decomposition scales with \(n^3\) such that this approach is only feasible for relatively small datasets. However, intuitively, only local correlations are important for the Gaussian process (at least for the squared exponential kernel we used above). Points separated by a large distance are virtually uncorrelated but these negligible correlations nevertheless slow down our calculations. Fortunately, we can factorize the joint distribution to obtain a product of conditional distributions.

\[ p\left(\eta\mid \Sigma\right)=\prod_{i=1}^n p\left(\eta_{i}\mid\eta_{<i},\Sigma\right), \]

where \(\eta_{<i}\) denotes all parameters \(\left\{\eta_1, \ldots \eta_{i-1}\right\}\) preceeding \(\eta_i\). In a nearest neighbor Gaussian process, we approximate the full conditional by only conditioning on at most the \(k\) preceeding parameters, i.e., we neglect all correlations between parameters that are further apart than the \(k^\text{th}\) nearest neighbor. We thus need to decompose \(n\) square matrices each having \(k + 1\) rows. The computational cost scales as \(nk^3\)–a substantial saving if \(n\) is large compared with \(k\).

More generally, we can construct a Gaussian process on a directed acyclic graph using the factorization into conditionals. Here, we employ this general formulation and represent the nearest neighbor Gaussian process as a graph. The conditioning structure is illustrated below.

from gptools.util.graph import lattice_predecessors, predecessors_to_edge_index

k = 5
predecessors = lattice_predecessors((n,), k)

idx = 20
fig, ax = plt.subplots()
ax.scatter(predecessors[idx, -1] - np.arange(3) - 1, np.zeros(3), color="silver",
           label="ignored nodes")
ax.scatter(idx + np.arange(3) + 1, np.zeros(3), color="silver")
ax.scatter(predecessors[idx], np.zeros(k), label="conditioning nodes")
ax.scatter(idx, 0, label="target node")
ax.legend()
for j in predecessors[idx, :-1]:
    ax.annotate("", xy=(j, 0), xytext=(idx, 0),
                arrowprops={"arrowstyle": "-|>", "connectionstyle": "arc3,rad=.5"})
ax.set_ylim(-.25, .25)
ax.set_axis_off()
fig.tight_layout()

The function lattice_predecessors constructs a nearest neighbor matrix predecessors with shape (n, k + 1) such that the first k elements of the \(i^\text{th}\) row are the predecessors of node \(i\) and the last element is the node itself. The corresponding Stan program is shown below. The distribution gp_graph_exp_quad_cov_lpdf(vector, vector, array [] vector, real, real, array [,] int) encodes the Gaussian process.

// Graph gaussian process with log link for Poisson observations.

functions {
    #include gptools/util.stan
    #include gptools/graph.stan
}

data {
    #include data.stan
    // Information about the graph.
    int num_edges;
    array [2, num_edges] int edge_index;
}

transformed data {
    array [n] int degrees = out_degrees(n, edge_index);
}

parameters {
    vector[n] eta;
}

model {
    // Graph Gaussian process prior and observation model.
    eta ~ gp_graph_exp_quad_cov(zeros_vector(n), X, sigma, length_scale, edge_index, degrees,
                                epsilon);
    y ~ poisson_log(eta);
}

edges = predecessors_to_edge_index(predecessors)
sample.update({"edge_index": edges, "num_edges": edges.shape[1]})
sample_and_plot("poisson_regression_graph_centered.stan", sample)

WARNING:root:cmdstan<2.36.0 had a bug in the evaluation of Matern 3/2 kernels with different length scales for different dimensions; see https://github.com/stan-dev/math/pull/3084 for details. Your model may yield unexpected results or crash if you use nearest-neighbor Gaussian processes with Matern 3/2 kernels. Your cmdstan version is 2.31; consider upgrading.

WARNING:cmdstanpy:CmdStanModel(compile=...) is deprecated and will be removed in the next major version. The constructor will always ensure a model has a compiled executable.
If you wish to force recompilation, use force_compile=True instead.

INFO:cmdstanpy:compiling stan file /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_centered.stan to exe file /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_centered

DEBUG:cmdstanpy:cmd: make STANCFLAGS+=--filename-in-msg=poisson_regression_graph_centered.stan STANCFLAGS+=--include-paths=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/python/stan/gptools/stan,/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_centered
cwd: /home/docs/.cmdstan/cmdstan-2.31.0

DEBUG:cmdstanpy:Console output:

--- Translating Stan model to C++ code ---
bin/stanc --filename-in-msg=poisson_regression_graph_centered.stan --include-paths=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/python/stan/gptools/stan,/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression --o=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_centered.hpp /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_centered.stan

--- Compiling, linking C++ code ---
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.9 -I stan/lib/stan_math/lib/boost_1.78.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS          -c -Wno-ignored-attributes   -x c++ -o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_centered.o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_centered.hpp
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.9 -I stan/lib/stan_math/lib/boost_1.78.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS                -Wl,-L,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb"      /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_centered.o src/cmdstan/main.o        -Wl,-L,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb"   stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_nvecserial.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_cvodes.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_idas.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_kinsol.a  stan/lib/stan_math/lib/tbb/libtbb.so.2 -o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_centered
rm -f /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_centered.o

INFO:cmdstanpy:compiled model executable: /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_centered

DEBUG:cmdstanpy:cmd: /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_centered info
cwd: None

DEBUG:cmdstanpy:input tempfile: /tmp/tmpikaxryv3/ndfrnusj.json

INFO:cmdstanpy:CmdStan start processing

DEBUG:cmdstanpy:idx 0

DEBUG:cmdstanpy:running CmdStan, num_threads: 1

DEBUG:cmdstanpy:CmdStan args: ['/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_centered', 'id=1', 'random', 'seed=7930', 'data', 'file=/tmp/tmpikaxryv3/ndfrnusj.json', 'output', 'file=/tmp/tmpikaxryv3/poisson_regression_graph_centeredbqvr7jaw/poisson_regression_graph_centered-20250317164945.csv', 'refresh=10', 'method=sample', 'num_samples=100', 'num_warmup=100', 'algorithm=hmc', 'adapt', 'engaged=1']

INFO:cmdstanpy:CmdStan done processing.

DEBUG:cmdstanpy:runset
RunSet: chains=1, chain_ids=[1], num_processes=1
 cmd (chain 1):
	['/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_centered', 'id=1', 'random', 'seed=7930', 'data', 'file=/tmp/tmpikaxryv3/ndfrnusj.json', 'output', 'file=/tmp/tmpikaxryv3/poisson_regression_graph_centeredbqvr7jaw/poisson_regression_graph_centered-20250317164945.csv', 'refresh=10', 'method=sample', 'num_samples=100', 'num_warmup=100', 'algorithm=hmc', 'adapt', 'engaged=1']
 retcodes=[0]
 per-chain output files (showing chain 1 only):
 csv_file:
	/tmp/tmpikaxryv3/poisson_regression_graph_centeredbqvr7jaw/poisson_regression_graph_centered-20250317164945.csv
 console_msgs (if any):
	/tmp/tmpikaxryv3/poisson_regression_graph_centeredbqvr7jaw/poisson_regression_graph_centered-20250317164945_0-stdout.txt

DEBUG:cmdstanpy:Chain 1 console:
method = sample (Default)
  sample
    num_samples = 100
    num_warmup = 100
    save_warmup = 0 (Default)
    thin = 1 (Default)
    adapt
      engaged = 1 (Default)
      gamma = 0.050000000000000003 (Default)
      delta = 0.80000000000000004 (Default)
      kappa = 0.75 (Default)
      t0 = 10 (Default)
      init_buffer = 75 (Default)
      term_buffer = 50 (Default)
      window = 25 (Default)
    algorithm = hmc (Default)
      hmc
        engine = nuts (Default)
          nuts
            max_depth = 10 (Default)
        metric = diag_e (Default)
        metric_file =  (Default)
        stepsize = 1 (Default)
        stepsize_jitter = 0 (Default)
    num_chains = 1 (Default)
id = 1 (Default)
data
  file = /tmp/tmpikaxryv3/ndfrnusj.json
init = 2 (Default)
random
  seed = 7930
output
  file = /tmp/tmpikaxryv3/poisson_regression_graph_centeredbqvr7jaw/poisson_regression_graph_centered-20250317164945.csv
  diagnostic_file =  (Default)
  refresh = 10
  sig_figs = -1 (Default)
  profile_file = profile.csv (Default)
num_threads = 1 (Default)


Gradient evaluation took 0.000343 seconds
1000 transitions using 10 leapfrog steps per transition would take 3.43 seconds.
Adjust your expectations accordingly!


WARNING: There aren't enough warmup iterations to fit the
         three stages of adaptation as currently configured.
         Reducing each adaptation stage to 15%/75%/10% of
         the given number of warmup iterations:
           init_buffer = 15
           adapt_window = 75
           term_buffer = 10

Iteration:   1 / 200 [  0%]  (Warmup)
Iteration:  10 / 200 [  5%]  (Warmup)
Iteration:  20 / 200 [ 10%]  (Warmup)
Iteration:  30 / 200 [ 15%]  (Warmup)
Iteration:  40 / 200 [ 20%]  (Warmup)
Iteration:  50 / 200 [ 25%]  (Warmup)
Iteration:  60 / 200 [ 30%]  (Warmup)
Iteration:  70 / 200 [ 35%]  (Warmup)
Iteration:  80 / 200 [ 40%]  (Warmup)
Iteration:  90 / 200 [ 45%]  (Warmup)
Iteration: 100 / 200 [ 50%]  (Warmup)
Iteration: 101 / 200 [ 50%]  (Sampling)
Iteration: 110 / 200 [ 55%]  (Sampling)
Iteration: 120 / 200 [ 60%]  (Sampling)
Iteration: 130 / 200 [ 65%]  (Sampling)
Iteration: 140 / 200 [ 70%]  (Sampling)
Iteration: 150 / 200 [ 75%]  (Sampling)
Iteration: 160 / 200 [ 80%]  (Sampling)
Iteration: 170 / 200 [ 85%]  (Sampling)
Iteration: 180 / 200 [ 90%]  (Sampling)
Iteration: 190 / 200 [ 95%]  (Sampling)
Iteration: 200 / 200 [100%]  (Sampling)

 Elapsed Time: 4.346 seconds (Warm-up)
               4.255 seconds (Sampling)
               8.601 seconds (Total)

sampled using poisson_regression_graph_centered.stan in 8.635 seconds

The centered graph Gaussian process is indeed faster than the standard implementation. However, it suffers from the same challenges: the posterior for \(\eta\) is highly correlated. Fortunately, we can also construct a non-centered parameterization.

// Graph gaussian process with log link for Poisson observations.

functions {
    #include gptools/util.stan
    #include gptools/graph.stan
}

data {
    #include data.stan
    // Information about the graph.
    int num_edges;
    array [2, num_edges] int edge_index;
}

transformed data {
    array [n] int degrees = out_degrees(n, edge_index);
}

parameters {
    vector[n] z;
}

transformed parameters {
    // Transform white noise to a sample from the graph Gaussian process.
    vector[n] eta = gp_inv_graph_exp_quad_cov(
        z, zeros_vector(n), X, sigma, length_scale, edge_index, degrees, epsilon);
}

model {
    // White noise prior (implies eta ~ gp_graph_exp_quad_cov(...)) and observation model.
    z ~ normal(0, 1);
    y ~ poisson_log(eta);
}

edges = predecessors_to_edge_index(predecessors)
sample.update({"edge_index": edges, "num_edges": edges.shape[1]})
sample_and_plot("poisson_regression_graph_non_centered.stan", sample)

WARNING:root:cmdstan<2.36.0 had a bug in the evaluation of Matern 3/2 kernels with different length scales for different dimensions; see https://github.com/stan-dev/math/pull/3084 for details. Your model may yield unexpected results or crash if you use nearest-neighbor Gaussian processes with Matern 3/2 kernels. Your cmdstan version is 2.31; consider upgrading.

WARNING:cmdstanpy:CmdStanModel(compile=...) is deprecated and will be removed in the next major version. The constructor will always ensure a model has a compiled executable.
If you wish to force recompilation, use force_compile=True instead.

INFO:cmdstanpy:compiling stan file /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_non_centered.stan to exe file /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_non_centered

DEBUG:cmdstanpy:cmd: make STANCFLAGS+=--filename-in-msg=poisson_regression_graph_non_centered.stan STANCFLAGS+=--include-paths=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/python/stan/gptools/stan,/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_non_centered
cwd: /home/docs/.cmdstan/cmdstan-2.31.0

DEBUG:cmdstanpy:Console output:

--- Translating Stan model to C++ code ---
bin/stanc --filename-in-msg=poisson_regression_graph_non_centered.stan --include-paths=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/python/stan/gptools/stan,/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression --o=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_non_centered.hpp /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_non_centered.stan

--- Compiling, linking C++ code ---
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.9 -I stan/lib/stan_math/lib/boost_1.78.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS          -c -Wno-ignored-attributes   -x c++ -o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_non_centered.o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_non_centered.hpp
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.9 -I stan/lib/stan_math/lib/boost_1.78.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS                -Wl,-L,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb"      /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_non_centered.o src/cmdstan/main.o        -Wl,-L,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb"   stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_nvecserial.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_cvodes.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_idas.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_kinsol.a  stan/lib/stan_math/lib/tbb/libtbb.so.2 -o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_non_centered
rm -f /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_non_centered.o

INFO:cmdstanpy:compiled model executable: /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_non_centered

DEBUG:cmdstanpy:cmd: /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_non_centered info
cwd: None

DEBUG:cmdstanpy:input tempfile: /tmp/tmpikaxryv3/ckwtnceb.json

INFO:cmdstanpy:CmdStan start processing

DEBUG:cmdstanpy:idx 0

DEBUG:cmdstanpy:running CmdStan, num_threads: 1

DEBUG:cmdstanpy:CmdStan args: ['/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_non_centered', 'id=1', 'random', 'seed=77413', 'data', 'file=/tmp/tmpikaxryv3/ckwtnceb.json', 'output', 'file=/tmp/tmpikaxryv3/poisson_regression_graph_non_centeredf5ffkyt6/poisson_regression_graph_non_centered-20250317165017.csv', 'refresh=10', 'method=sample', 'num_samples=100', 'num_warmup=100', 'algorithm=hmc', 'adapt', 'engaged=1']

INFO:cmdstanpy:CmdStan done processing.

DEBUG:cmdstanpy:runset
RunSet: chains=1, chain_ids=[1], num_processes=1
 cmd (chain 1):
	['/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_graph_non_centered', 'id=1', 'random', 'seed=77413', 'data', 'file=/tmp/tmpikaxryv3/ckwtnceb.json', 'output', 'file=/tmp/tmpikaxryv3/poisson_regression_graph_non_centeredf5ffkyt6/poisson_regression_graph_non_centered-20250317165017.csv', 'refresh=10', 'method=sample', 'num_samples=100', 'num_warmup=100', 'algorithm=hmc', 'adapt', 'engaged=1']
 retcodes=[0]
 per-chain output files (showing chain 1 only):
 csv_file:
	/tmp/tmpikaxryv3/poisson_regression_graph_non_centeredf5ffkyt6/poisson_regression_graph_non_centered-20250317165017.csv
 console_msgs (if any):
	/tmp/tmpikaxryv3/poisson_regression_graph_non_centeredf5ffkyt6/poisson_regression_graph_non_centered-20250317165017_0-stdout.txt

DEBUG:cmdstanpy:Chain 1 console:
method = sample (Default)
  sample
    num_samples = 100
    num_warmup = 100
    save_warmup = 0 (Default)
    thin = 1 (Default)
    adapt
      engaged = 1 (Default)
      gamma = 0.050000000000000003 (Default)
      delta = 0.80000000000000004 (Default)
      kappa = 0.75 (Default)
      t0 = 10 (Default)
      init_buffer = 75 (Default)
      term_buffer = 50 (Default)
      window = 25 (Default)
    algorithm = hmc (Default)
      hmc
        engine = nuts (Default)
          nuts
            max_depth = 10 (Default)
        metric = diag_e (Default)
        metric_file =  (Default)
        stepsize = 1 (Default)
        stepsize_jitter = 0 (Default)
    num_chains = 1 (Default)
id = 1 (Default)
data
  file = /tmp/tmpikaxryv3/ckwtnceb.json
init = 2 (Default)
random
  seed = 77413
output
  file = /tmp/tmpikaxryv3/poisson_regression_graph_non_centeredf5ffkyt6/poisson_regression_graph_non_centered-20250317165017.csv
  diagnostic_file =  (Default)
  refresh = 10
  sig_figs = -1 (Default)
  profile_file = profile.csv (Default)
num_threads = 1 (Default)


Gradient evaluation took 0.000354 seconds
1000 transitions using 10 leapfrog steps per transition would take 3.54 seconds.
Adjust your expectations accordingly!


WARNING: There aren't enough warmup iterations to fit the
         three stages of adaptation as currently configured.
         Reducing each adaptation stage to 15%/75%/10% of
         the given number of warmup iterations:
           init_buffer = 15
           adapt_window = 75
           term_buffer = 10

Iteration:   1 / 200 [  0%]  (Warmup)
Iteration:  10 / 200 [  5%]  (Warmup)
Iteration:  20 / 200 [ 10%]  (Warmup)
Iteration:  30 / 200 [ 15%]  (Warmup)
Iteration:  40 / 200 [ 20%]  (Warmup)
Iteration:  50 / 200 [ 25%]  (Warmup)
Iteration:  60 / 200 [ 30%]  (Warmup)
Iteration:  70 / 200 [ 35%]  (Warmup)
Iteration:  80 / 200 [ 40%]  (Warmup)
Iteration:  90 / 200 [ 45%]  (Warmup)
Iteration: 100 / 200 [ 50%]  (Warmup)
Iteration: 101 / 200 [ 50%]  (Sampling)
Iteration: 110 / 200 [ 55%]  (Sampling)
Iteration: 120 / 200 [ 60%]  (Sampling)
Iteration: 130 / 200 [ 65%]  (Sampling)
Iteration: 140 / 200 [ 70%]  (Sampling)
Iteration: 150 / 200 [ 75%]  (Sampling)
Iteration: 160 / 200 [ 80%]  (Sampling)
Iteration: 170 / 200 [ 85%]  (Sampling)
Iteration: 180 / 200 [ 90%]  (Sampling)
Iteration: 190 / 200 [ 95%]  (Sampling)
Iteration: 200 / 200 [100%]  (Sampling)

 Elapsed Time: 0.594 seconds (Warm-up)
               0.591 seconds (Sampling)
               1.185 seconds (Total)

sampled using poisson_regression_graph_non_centered.stan in 1.219 seconds

The non-centered parameterization of the graph Gaussian process is even faster, and we have reduced the runtime by more than an order of magnitude for this small dataset.

Gaussian process using fast Fourier transforms

If, as in this example, observations are regularly spaced, we can evaluate the likelihood using the fast Fourier transform. Similarly, we can construct a non-centered parameterization by transforming Fourier coefficients with gp_inv_rfft.

// Graph gaussian process with log link for Poisson observations.

functions {
    #include gptools/util.stan
    #include gptools/fft1.stan
}

data {
    #include data.stan
}

parameters {
    vector[n] eta;
}

transformed parameters {
    // Evaluate covariance of the point at zero with everything else.
    vector[n %/% 2 + 1] cov_rfft = gp_periodic_exp_quad_cov_rfft(n, sigma, length_scale, n)
        + epsilon;
}

model {
    // Fourier Gaussian process and observation model.
    eta ~ gp_rfft(zeros_vector(n), cov_rfft);
    y ~ poisson_log(eta);
}

sample_and_plot("poisson_regression_fourier_centered.stan", sample)

WARNING:root:cmdstan<2.36.0 had a bug in the evaluation of Matern 3/2 kernels with different length scales for different dimensions; see https://github.com/stan-dev/math/pull/3084 for details. Your model may yield unexpected results or crash if you use nearest-neighbor Gaussian processes with Matern 3/2 kernels. Your cmdstan version is 2.31; consider upgrading.

WARNING:cmdstanpy:CmdStanModel(compile=...) is deprecated and will be removed in the next major version. The constructor will always ensure a model has a compiled executable.
If you wish to force recompilation, use force_compile=True instead.

INFO:cmdstanpy:compiling stan file /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_centered.stan to exe file /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_centered

DEBUG:cmdstanpy:cmd: make STANCFLAGS+=--filename-in-msg=poisson_regression_fourier_centered.stan STANCFLAGS+=--include-paths=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/python/stan/gptools/stan,/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_centered
cwd: /home/docs/.cmdstan/cmdstan-2.31.0

DEBUG:cmdstanpy:Console output:

--- Translating Stan model to C++ code ---
bin/stanc --filename-in-msg=poisson_regression_fourier_centered.stan --include-paths=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/python/stan/gptools/stan,/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression --o=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_centered.hpp /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_centered.stan

--- Compiling, linking C++ code ---
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.9 -I stan/lib/stan_math/lib/boost_1.78.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS          -c -Wno-ignored-attributes   -x c++ -o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_centered.o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_centered.hpp
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.9 -I stan/lib/stan_math/lib/boost_1.78.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS                -Wl,-L,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb"      /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_centered.o src/cmdstan/main.o        -Wl,-L,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb"   stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_nvecserial.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_cvodes.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_idas.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_kinsol.a  stan/lib/stan_math/lib/tbb/libtbb.so.2 -o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_centered
rm -f /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_centered.o

INFO:cmdstanpy:compiled model executable: /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_centered

DEBUG:cmdstanpy:cmd: /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_centered info
cwd: None

DEBUG:cmdstanpy:input tempfile: /tmp/tmpikaxryv3/48asbeem.json

INFO:cmdstanpy:CmdStan start processing

DEBUG:cmdstanpy:idx 0

DEBUG:cmdstanpy:running CmdStan, num_threads: 1

DEBUG:cmdstanpy:CmdStan args: ['/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_centered', 'id=1', 'random', 'seed=54856', 'data', 'file=/tmp/tmpikaxryv3/48asbeem.json', 'output', 'file=/tmp/tmpikaxryv3/poisson_regression_fourier_centeredswqez1pv/poisson_regression_fourier_centered-20250317165039.csv', 'refresh=10', 'method=sample', 'num_samples=100', 'num_warmup=100', 'algorithm=hmc', 'adapt', 'engaged=1']

INFO:cmdstanpy:CmdStan done processing.

DEBUG:cmdstanpy:runset
RunSet: chains=1, chain_ids=[1], num_processes=1
 cmd (chain 1):
	['/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_centered', 'id=1', 'random', 'seed=54856', 'data', 'file=/tmp/tmpikaxryv3/48asbeem.json', 'output', 'file=/tmp/tmpikaxryv3/poisson_regression_fourier_centeredswqez1pv/poisson_regression_fourier_centered-20250317165039.csv', 'refresh=10', 'method=sample', 'num_samples=100', 'num_warmup=100', 'algorithm=hmc', 'adapt', 'engaged=1']
 retcodes=[0]
 per-chain output files (showing chain 1 only):
 csv_file:
	/tmp/tmpikaxryv3/poisson_regression_fourier_centeredswqez1pv/poisson_regression_fourier_centered-20250317165039.csv
 console_msgs (if any):
	/tmp/tmpikaxryv3/poisson_regression_fourier_centeredswqez1pv/poisson_regression_fourier_centered-20250317165039_0-stdout.txt

DEBUG:cmdstanpy:Chain 1 console:
method = sample (Default)
  sample
    num_samples = 100
    num_warmup = 100
    save_warmup = 0 (Default)
    thin = 1 (Default)
    adapt
      engaged = 1 (Default)
      gamma = 0.050000000000000003 (Default)
      delta = 0.80000000000000004 (Default)
      kappa = 0.75 (Default)
      t0 = 10 (Default)
      init_buffer = 75 (Default)
      term_buffer = 50 (Default)
      window = 25 (Default)
    algorithm = hmc (Default)
      hmc
        engine = nuts (Default)
          nuts
            max_depth = 10 (Default)
        metric = diag_e (Default)
        metric_file =  (Default)
        stepsize = 1 (Default)
        stepsize_jitter = 0 (Default)
    num_chains = 1 (Default)
id = 1 (Default)
data
  file = /tmp/tmpikaxryv3/48asbeem.json
init = 2 (Default)
random
  seed = 54856
output
  file = /tmp/tmpikaxryv3/poisson_regression_fourier_centeredswqez1pv/poisson_regression_fourier_centered-20250317165039.csv
  diagnostic_file =  (Default)
  refresh = 10
  sig_figs = -1 (Default)
  profile_file = profile.csv (Default)
num_threads = 1 (Default)


Gradient evaluation took 5.1e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.51 seconds.
Adjust your expectations accordingly!


WARNING: There aren't enough warmup iterations to fit the
         three stages of adaptation as currently configured.
         Reducing each adaptation stage to 15%/75%/10% of
         the given number of warmup iterations:
           init_buffer = 15
           adapt_window = 75
           term_buffer = 10

Iteration:   1 / 200 [  0%]  (Warmup)
Iteration:  10 / 200 [  5%]  (Warmup)
Iteration:  20 / 200 [ 10%]  (Warmup)
Iteration:  30 / 200 [ 15%]  (Warmup)
Iteration:  40 / 200 [ 20%]  (Warmup)
Iteration:  50 / 200 [ 25%]  (Warmup)
Iteration:  60 / 200 [ 30%]  (Warmup)
Iteration:  70 / 200 [ 35%]  (Warmup)
Iteration:  80 / 200 [ 40%]  (Warmup)
Iteration:  90 / 200 [ 45%]  (Warmup)
Iteration: 100 / 200 [ 50%]  (Warmup)
Iteration: 101 / 200 [ 50%]  (Sampling)
Iteration: 110 / 200 [ 55%]  (Sampling)
Iteration: 120 / 200 [ 60%]  (Sampling)
Iteration: 130 / 200 [ 65%]  (Sampling)
Iteration: 140 / 200 [ 70%]  (Sampling)
Iteration: 150 / 200 [ 75%]  (Sampling)
Iteration: 160 / 200 [ 80%]  (Sampling)
Iteration: 170 / 200 [ 85%]  (Sampling)
Iteration: 180 / 200 [ 90%]  (Sampling)
Iteration: 190 / 200 [ 95%]  (Sampling)
Iteration: 200 / 200 [100%]  (Sampling)

 Elapsed Time: 0.285 seconds (Warm-up)
               0.251 seconds (Sampling)
               0.536 seconds (Total)

sampled using poisson_regression_fourier_centered.stan in 0.567 seconds

sample_and_plot("poisson_regression_fourier_non_centered.stan", sample)

WARNING:root:cmdstan<2.36.0 had a bug in the evaluation of Matern 3/2 kernels with different length scales for different dimensions; see https://github.com/stan-dev/math/pull/3084 for details. Your model may yield unexpected results or crash if you use nearest-neighbor Gaussian processes with Matern 3/2 kernels. Your cmdstan version is 2.31; consider upgrading.

WARNING:cmdstanpy:CmdStanModel(compile=...) is deprecated and will be removed in the next major version. The constructor will always ensure a model has a compiled executable.
If you wish to force recompilation, use force_compile=True instead.

INFO:cmdstanpy:compiling stan file /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered.stan to exe file /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered

DEBUG:cmdstanpy:cmd: make STANCFLAGS+=--filename-in-msg=poisson_regression_fourier_non_centered.stan STANCFLAGS+=--include-paths=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/python/stan/gptools/stan,/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered
cwd: /home/docs/.cmdstan/cmdstan-2.31.0

DEBUG:cmdstanpy:Console output:

--- Translating Stan model to C++ code ---
bin/stanc --filename-in-msg=poisson_regression_fourier_non_centered.stan --include-paths=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/python/stan/gptools/stan,/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression --o=/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered.hpp /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered.stan

--- Compiling, linking C++ code ---
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.9 -I stan/lib/stan_math/lib/boost_1.78.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS          -c -Wno-ignored-attributes   -x c++ -o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered.o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered.hpp
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O3 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.9 -I stan/lib/stan_math/lib/boost_1.78.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS                -Wl,-L,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb"      /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered.o src/cmdstan/main.o        -Wl,-L,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb" -Wl,-rpath,"/home/docs/.cmdstan/cmdstan-2.31.0/stan/lib/stan_math/lib/tbb"   stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_nvecserial.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_cvodes.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_idas.a stan/lib/stan_math/lib/sundials_6.1.1/lib/libsundials_kinsol.a  stan/lib/stan_math/lib/tbb/libtbb.so.2 -o /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered
rm -f /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered.o

INFO:cmdstanpy:compiled model executable: /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered

DEBUG:cmdstanpy:cmd: /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered info
cwd: None

DEBUG:cmdstanpy:input tempfile: /tmp/tmpikaxryv3/z7iof_u4.json

INFO:cmdstanpy:CmdStan start processing

DEBUG:cmdstanpy:idx 0

DEBUG:cmdstanpy:running CmdStan, num_threads: 1

DEBUG:cmdstanpy:CmdStan args: ['/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered', 'id=1', 'random', 'seed=92804', 'data', 'file=/tmp/tmpikaxryv3/z7iof_u4.json', 'output', 'file=/tmp/tmpikaxryv3/poisson_regression_fourier_non_centeredhgg5d4sh/poisson_regression_fourier_non_centered-20250317165102.csv', 'refresh=10', 'method=sample', 'num_samples=100', 'num_warmup=100', 'algorithm=hmc', 'adapt', 'engaged=1']

INFO:cmdstanpy:CmdStan done processing.

DEBUG:cmdstanpy:runset
RunSet: chains=1, chain_ids=[1], num_processes=1
 cmd (chain 1):
	['/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered', 'id=1', 'random', 'seed=92804', 'data', 'file=/tmp/tmpikaxryv3/z7iof_u4.json', 'output', 'file=/tmp/tmpikaxryv3/poisson_regression_fourier_non_centeredhgg5d4sh/poisson_regression_fourier_non_centered-20250317165102.csv', 'refresh=10', 'method=sample', 'num_samples=100', 'num_warmup=100', 'algorithm=hmc', 'adapt', 'engaged=1']
 retcodes=[0]
 per-chain output files (showing chain 1 only):
 csv_file:
	/tmp/tmpikaxryv3/poisson_regression_fourier_non_centeredhgg5d4sh/poisson_regression_fourier_non_centered-20250317165102.csv
 console_msgs (if any):
	/tmp/tmpikaxryv3/poisson_regression_fourier_non_centeredhgg5d4sh/poisson_regression_fourier_non_centered-20250317165102_0-stdout.txt

DEBUG:cmdstanpy:Chain 1 console:
method = sample (Default)
  sample
    num_samples = 100
    num_warmup = 100
    save_warmup = 0 (Default)
    thin = 1 (Default)
    adapt
      engaged = 1 (Default)
      gamma = 0.050000000000000003 (Default)
      delta = 0.80000000000000004 (Default)
      kappa = 0.75 (Default)
      t0 = 10 (Default)
      init_buffer = 75 (Default)
      term_buffer = 50 (Default)
      window = 25 (Default)
    algorithm = hmc (Default)
      hmc
        engine = nuts (Default)
          nuts
            max_depth = 10 (Default)
        metric = diag_e (Default)
        metric_file =  (Default)
        stepsize = 1 (Default)
        stepsize_jitter = 0 (Default)
    num_chains = 1 (Default)
id = 1 (Default)
data
  file = /tmp/tmpikaxryv3/z7iof_u4.json
init = 2 (Default)
random
  seed = 92804
output
  file = /tmp/tmpikaxryv3/poisson_regression_fourier_non_centeredhgg5d4sh/poisson_regression_fourier_non_centered-20250317165102.csv
  diagnostic_file =  (Default)
  refresh = 10
  sig_figs = -1 (Default)
  profile_file = profile.csv (Default)
num_threads = 1 (Default)


Gradient evaluation took 5e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.5 seconds.
Adjust your expectations accordingly!


WARNING: There aren't enough warmup iterations to fit the
         three stages of adaptation as currently configured.
         Reducing each adaptation stage to 15%/75%/10% of
         the given number of warmup iterations:
           init_buffer = 15
           adapt_window = 75
           term_buffer = 10

Iteration:   1 / 200 [  0%]  (Warmup)
Iteration:  10 / 200 [  5%]  (Warmup)
Iteration:  20 / 200 [ 10%]  (Warmup)
Iteration:  30 / 200 [ 15%]  (Warmup)
Iteration:  40 / 200 [ 20%]  (Warmup)
Iteration:  50 / 200 [ 25%]  (Warmup)
Iteration:  60 / 200 [ 30%]  (Warmup)
Iteration:  70 / 200 [ 35%]  (Warmup)
Iteration:  80 / 200 [ 40%]  (Warmup)
Iteration:  90 / 200 [ 45%]  (Warmup)
Iteration: 100 / 200 [ 50%]  (Warmup)
Iteration: 101 / 200 [ 50%]  (Sampling)
Iteration: 110 / 200 [ 55%]  (Sampling)
Iteration: 120 / 200 [ 60%]  (Sampling)
Iteration: 130 / 200 [ 65%]  (Sampling)
Iteration: 140 / 200 [ 70%]  (Sampling)
Iteration: 150 / 200 [ 75%]  (Sampling)
Iteration: 160 / 200 [ 80%]  (Sampling)
Iteration: 170 / 200 [ 85%]  (Sampling)
Iteration: 180 / 200 [ 90%]  (Sampling)
Iteration: 190 / 200 [ 95%]  (Sampling)
Iteration: 200 / 200 [100%]  (Sampling)

 Elapsed Time: 0.046 seconds (Warm-up)
               0.053 seconds (Sampling)
               0.099 seconds (Total)

sampled using poisson_regression_fourier_non_centered.stan in 0.140 seconds

Given the substantial performance improvements, we can readily increase the sample size as illustrated below for a dataset with more than a thousand observations. Sampling from the forward model takes a while because evaluating the periodic squared exponential kernel requires a series approximation.

x = np.arange(1024)
kernel = ExpQuadKernel(sigma=1.2, length_scale=15, period=x.size) + DiagonalKernel(1e-3, x.size)
sample = simulate(x, kernel)
plot_sample(sample).legend()

<matplotlib.legend.Legend at 0x710ab46447c0>

sample_and_plot("poisson_regression_fourier_non_centered.stan", sample)

WARNING:root:cmdstan<2.36.0 had a bug in the evaluation of Matern 3/2 kernels with different length scales for different dimensions; see https://github.com/stan-dev/math/pull/3084 for details. Your model may yield unexpected results or crash if you use nearest-neighbor Gaussian processes with Matern 3/2 kernels. Your cmdstan version is 2.31; consider upgrading.

WARNING:cmdstanpy:CmdStanModel(compile=...) is deprecated and will be removed in the next major version. The constructor will always ensure a model has a compiled executable.
If you wish to force recompilation, use force_compile=True instead.

DEBUG:cmdstanpy:found newer exe file, not recompiling

DEBUG:cmdstanpy:cmd: /home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered info
cwd: None

DEBUG:cmdstanpy:input tempfile: /tmp/tmpikaxryv3/glhk1e8b.json

INFO:cmdstanpy:CmdStan start processing

DEBUG:cmdstanpy:idx 0

DEBUG:cmdstanpy:running CmdStan, num_threads: 1

DEBUG:cmdstanpy:CmdStan args: ['/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered', 'id=1', 'random', 'seed=96484', 'data', 'file=/tmp/tmpikaxryv3/glhk1e8b.json', 'output', 'file=/tmp/tmpikaxryv3/poisson_regression_fourier_non_centeredwe9n9wty/poisson_regression_fourier_non_centered-20250317165103.csv', 'refresh=10', 'method=sample', 'num_samples=100', 'num_warmup=100', 'algorithm=hmc', 'adapt', 'engaged=1']

INFO:cmdstanpy:CmdStan done processing.

DEBUG:cmdstanpy:runset
RunSet: chains=1, chain_ids=[1], num_processes=1
 cmd (chain 1):
	['/home/docs/checkouts/readthedocs.org/user_builds/gptools-stan/checkouts/latest/docs/poisson_regression/poisson_regression_fourier_non_centered', 'id=1', 'random', 'seed=96484', 'data', 'file=/tmp/tmpikaxryv3/glhk1e8b.json', 'output', 'file=/tmp/tmpikaxryv3/poisson_regression_fourier_non_centeredwe9n9wty/poisson_regression_fourier_non_centered-20250317165103.csv', 'refresh=10', 'method=sample', 'num_samples=100', 'num_warmup=100', 'algorithm=hmc', 'adapt', 'engaged=1']
 retcodes=[0]
 per-chain output files (showing chain 1 only):
 csv_file:
	/tmp/tmpikaxryv3/poisson_regression_fourier_non_centeredwe9n9wty/poisson_regression_fourier_non_centered-20250317165103.csv
 console_msgs (if any):
	/tmp/tmpikaxryv3/poisson_regression_fourier_non_centeredwe9n9wty/poisson_regression_fourier_non_centered-20250317165103_0-stdout.txt

DEBUG:cmdstanpy:Chain 1 console:
method = sample (Default)
  sample
    num_samples = 100
    num_warmup = 100
    save_warmup = 0 (Default)
    thin = 1 (Default)
    adapt
      engaged = 1 (Default)
      gamma = 0.050000000000000003 (Default)
      delta = 0.80000000000000004 (Default)
      kappa = 0.75 (Default)
      t0 = 10 (Default)
      init_buffer = 75 (Default)
      term_buffer = 50 (Default)
      window = 25 (Default)
    algorithm = hmc (Default)
      hmc
        engine = nuts (Default)
          nuts
            max_depth = 10 (Default)
        metric = diag_e (Default)
        metric_file =  (Default)
        stepsize = 1 (Default)
        stepsize_jitter = 0 (Default)
    num_chains = 1 (Default)
id = 1 (Default)
data
  file = /tmp/tmpikaxryv3/glhk1e8b.json
init = 2 (Default)
random
  seed = 96484
output
  file = /tmp/tmpikaxryv3/poisson_regression_fourier_non_centeredwe9n9wty/poisson_regression_fourier_non_centered-20250317165103.csv
  diagnostic_file =  (Default)
  refresh = 10
  sig_figs = -1 (Default)
  profile_file = profile.csv (Default)
num_threads = 1 (Default)


Gradient evaluation took 0.000572 seconds
1000 transitions using 10 leapfrog steps per transition would take 5.72 seconds.
Adjust your expectations accordingly!


WARNING: There aren't enough warmup iterations to fit the
         three stages of adaptation as currently configured.
         Reducing each adaptation stage to 15%/75%/10% of
         the given number of warmup iterations:
           init_buffer = 15
           adapt_window = 75
           term_buffer = 10

Iteration:   1 / 200 [  0%]  (Warmup)
Iteration:  10 / 200 [  5%]  (Warmup)
Iteration:  20 / 200 [ 10%]  (Warmup)
Iteration:  30 / 200 [ 15%]  (Warmup)
Iteration:  40 / 200 [ 20%]  (Warmup)
Iteration:  50 / 200 [ 25%]  (Warmup)
Iteration:  60 / 200 [ 30%]  (Warmup)
Iteration:  70 / 200 [ 35%]  (Warmup)
Iteration:  80 / 200 [ 40%]  (Warmup)
Iteration:  90 / 200 [ 45%]  (Warmup)
Iteration: 100 / 200 [ 50%]  (Warmup)
Iteration: 101 / 200 [ 50%]  (Sampling)
Iteration: 110 / 200 [ 55%]  (Sampling)
Iteration: 120 / 200 [ 60%]  (Sampling)
Iteration: 130 / 200 [ 65%]  (Sampling)
Iteration: 140 / 200 [ 70%]  (Sampling)
Iteration: 150 / 200 [ 75%]  (Sampling)
Iteration: 160 / 200 [ 80%]  (Sampling)
Iteration: 170 / 200 [ 85%]  (Sampling)
Iteration: 180 / 200 [ 90%]  (Sampling)
Iteration: 190 / 200 [ 95%]  (Sampling)
Iteration: 200 / 200 [100%]  (Sampling)

 Elapsed Time: 1.767 seconds (Warm-up)
               1.402 seconds (Sampling)
               3.169 seconds (Total)

sampled using poisson_regression_fourier_non_centered.stan in 3.222 seconds