Hi all,
I’ve been experimenting with running Stan inside webR and wanted to share a working proof-of-concept.
Thanks to @seantalts’ work on stanli (the Stan interpreter/runtime) and @andrjohns’ work on stanr and his stanli-backend fork of brms, it is now fully possible to compile/interpret and sample Stan models completely client-side in the browser.
I set up Wasm builds of the brms fork and the stanli R package in R-universe and tested the workflow across three levels:
1. Setup & Core Runtime (stanli)
stanli provides the lightweight runtime/interpreter for Stan in WebAssembly.
Workaround for browser CORS: In webR, dynamically downloading the precompiled runtime viastanli::stanli_install()currently gets blocked by browser CORS restrictions. The simple workaround for now is to grabstanli-runtime-emscripten-wasm32.tar.gz, unpacklibstanli.so, and upload/mount it into the virtual root (/libstanli.so). Note that this isn’t needed for usingstanr, sincestanrbundles the runtime forstanlidirectly.
# Install stanli from my r-universe
install.packages("stanli", repos = "https://seantalts.r-universe.dev")
# Point to the uploaded runtime
# Sys.setenv(STANLI_RUNTIME = "/libstanli.so")
# Not relevant anymore, the wasm build includes the stanli runtime
library(stanli)
stopifnot(stanli_available())
Direct sampling with stanli:
model_code <- "
data {
int<lower=0> N;
vector[N] y;
}
parameters {
real mu;
real<lower=0> sigma;
}
model {
mu ~ normal(0, 10);
sigma ~ exponential(1);
y ~ normal(mu, sigma);
}
"
model_data <- list(
N = 8L,
y = c(28, 8, -3, 7, -1, 1, 18, 12)
)
model <- stanli_model(code = model_code, data = model_data)
fit <- sample_model(model, chains = 2, warmup = 250, samples = 2500, seed = 123)
Results
> summary(fit)
variable mean mcse_mean sd mcse_sd q5 q50 q95 ess_bulk ess_tail
1 mu 8.069115 0.05625743 2.668532 0.04914493 3.702832 8.108840 12.46221 2291.173 1948.553
2 sigma 7.739020 0.02743319 1.344503 0.02867524 5.810334 7.605246 10.16252 2536.218 2211.829
rhat
1 1.000607
2 1.002873
> stanli_diagnose(fit)
No divergent transitions.
No transitions saturated the maximum treedepth of 10.
E-BFMI is above 0.3 in every chain.
R-hat is below 1.01 for every parameter (worst 1.003, sigma).
Bulk ESS is at least 100 per chain for every parameter (worst 2291, mu).
Tail ESS is at least 100 per chain for every parameter (worst 1949, mu).
No problems detected.
2. Modern Interface with stanr
Andrew’s stanr package gives a familiar cmdstanr-style API on top of stanli, with seamless integration into posterior and full diagnostics (note: stanr bundles the stanli runtime directly, so manual mounting or the stanli R package isn’t strictly required here):
install.packages("stanr", repos = c("https://andrjohns.r-universe.dev", "https://repo.r-wasm.org/"))
library(stanr)
mod <- stan_model(
code = model_code,
model_name = "normal_location_scale",
backend = "stanli"
)
fit <- mod$sample(
data = model_data,
chains = 4,
iter_warmup = 500,
iter_sampling = 1000,
seed = 123
)
## Prints sampling as usual
Gradient evaluation took 0.000315 seconds
1000 transitions using 10 leapfrog steps per transition would take 3.15 seconds.
Adjust your expectations accordingly!
Gradient evaluation took 1.5e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
Adjust your expectations accordingly!
Gradient evaluation took 1e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
Adjust your expectations accordingly!
Gradient evaluation took 1e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
Adjust your expectations accordingly!
Chain [1] Iteration: 1 / 1500 [ 0%] (Warmup)
Chain [1] Iteration: 100 / 1500 [ 6%] (Warmup)
Chain [1] Iteration: 200 / 1500 [ 13%] (Warmup)
Chain [1] Iteration: 300 / 1500 [ 20%] (Warmup)
Chain [1] Iteration: 400 / 1500 [ 26%] (Warmup)
Chain [1] Iteration: 500 / 1500 [ 33%] (Warmup)
Chain [1] Iteration: 501 / 1500 [ 33%] (Sampling)
Chain [1] Iteration: 600 / 1500 [ 40%] (Sampling)
Chain [1] Iteration: 700 / 1500 [ 46%] (Sampling)
Chain [1] Iteration: 800 / 1500 [ 53%] (Sampling)
Chain [1] Iteration: 900 / 1500 [ 60%] (Sampling)
Chain [1] Iteration: 1000 / 1500 [ 66%] (Sampling)
Chain [1] Iteration: 1100 / 1500 [ 73%] (Sampling)
Chain [1] Iteration: 1200 / 1500 [ 80%] (Sampling)
Chain [1] Iteration: 1300 / 1500 [ 86%] (Sampling)
Chain [1] Iteration: 1400 / 1500 [ 93%] (Sampling)
Chain [1] Iteration: 1500 / 1500 [100%] (Sampling)
Elapsed Time: 0.079 seconds (Warm-up)
0.088 seconds (Sampling)
0.167 seconds (Total)
Chain [2] Iteration: 1 / 1500 [ 0%] (Warmup)
Chain [2] Iteration: 100 / 1500 [ 6%] (Warmup)
Chain [2] Iteration: 200 / 1500 [ 13%] (Warmup)
Chain [2] Iteration: 300 / 1500 [ 20%] (Warmup)
Chain [2] Iteration: 400 / 1500 [ 26%] (Warmup)
Chain [2] Iteration: 500 / 1500 [ 33%] (Warmup)
Chain [2] Iteration: 501 / 1500 [ 33%] (Sampling)
Chain [2] Iteration: 600 / 1500 [ 40%] (Sampling)
Chain [2] Iteration: 700 / 1500 [ 46%] (Sampling)
Chain [2] Iteration: 800 / 1500 [ 53%] (Sampling)
Chain [2] Iteration: 900 / 1500 [ 60%] (Sampling)
Chain [2] Iteration: 1000 / 1500 [ 66%] (Sampling)
Chain [2] Iteration: 1100 / 1500 [ 73%] (Sampling)
Chain [2] Iteration: 1200 / 1500 [ 80%] (Sampling)
Chain [2] Iteration: 1300 / 1500 [ 86%] (Sampling)
Chain [2] Iteration: 1400 / 1500 [ 93%] (Sampling)
Chain [2] Iteration: 1500 / 1500 [100%] (Sampling)
Elapsed Time: 0.069 seconds (Warm-up)
0.103 seconds (Sampling)
0.172 seconds (Total)
Chain [3] Iteration: 1 / 1500 [ 0%] (Warmup)
Chain [3] Iteration: 100 / 1500 [ 6%] (Warmup)
Chain [3] Iteration: 200 / 1500 [ 13%] (Warmup)
Chain [3] Iteration: 300 / 1500 [ 20%] (Warmup)
Chain [3] Iteration: 400 / 1500 [ 26%] (Warmup)
Chain [3] Iteration: 500 / 1500 [ 33%] (Warmup)
Chain [3] Iteration: 501 / 1500 [ 33%] (Sampling)
Chain [3] Iteration: 600 / 1500 [ 40%] (Sampling)
Chain [3] Iteration: 700 / 1500 [ 46%] (Sampling)
Chain [3] Iteration: 800 / 1500 [ 53%] (Sampling)
Chain [3] Iteration: 900 / 1500 [ 60%] (Sampling)
Chain [3] Iteration: 1000 / 1500 [ 66%] (Sampling)
Chain [3] Iteration: 1100 / 1500 [ 73%] (Sampling)
Chain [3] Iteration: 1200 / 1500 [ 80%] (Sampling)
Chain [3] Iteration: 1300 / 1500 [ 86%] (Sampling)
Chain [3] Iteration: 1400 / 1500 [ 93%] (Sampling)
Chain [3] Iteration: 1500 / 1500 [100%] (Sampling)
Elapsed Time: 0.073 seconds (Warm-up)
0.094 seconds (Sampling)
0.167 seconds (Total)
Chain [4] Iteration: 1 / 1500 [ 0%] (Warmup)
Chain [4] Iteration: 100 / 1500 [ 6%] (Warmup)
Chain [4] Iteration: 200 / 1500 [ 13%] (Warmup)
Chain [4] Iteration: 300 / 1500 [ 20%] (Warmup)
Chain [4] Iteration: 400 / 1500 [ 26%] (Warmup)
Chain [4] Iteration: 500 / 1500 [ 33%] (Warmup)
Chain [4] Iteration: 501 / 1500 [ 33%] (Sampling)
Chain [4] Iteration: 600 / 1500 [ 40%] (Sampling)
Chain [4] Iteration: 700 / 1500 [ 46%] (Sampling)
Chain [4] Iteration: 800 / 1500 [ 53%] (Sampling)
Chain [4] Iteration: 900 / 1500 [ 60%] (Sampling)
Chain [4] Iteration: 1000 / 1500 [ 66%] (Sampling)
Chain [4] Iteration: 1100 / 1500 [ 73%] (Sampling)
Chain [4] Iteration: 1200 / 1500 [ 80%] (Sampling)
Chain [4] Iteration: 1300 / 1500 [ 86%] (Sampling)
Chain [4] Iteration: 1400 / 1500 [ 93%] (Sampling)
Chain [4] Iteration: 1500 / 1500 [100%] (Sampling)
Elapsed Time: 0.067 seconds (Warm-up)
0.089 seconds (Sampling)
0.156 seconds (Total)
# Summaries & Posterior draws
> fit$summary(variables = c("mu", "sigma"))
# A tibble: 2 × 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 mu 8.17 8.20 2.70 2.70 3.77 12.5 1.00 3176. 2325.
2 sigma 7.78 7.64 1.36 1.30 5.84 10.2 0.999 3037. 2627.
>
> fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0
$num_max_treedepth
[1] 0 0 0 0
$ebfmi
[1] 1.1990484 0.9455147 0.9590174 1.1145598
>
> draws <- fit$draws(variables = c("mu", "sigma"))
>
> draws
# A draws_array: 1000 iterations, 4 chains, and 2 variables
, , mu
1 2 3 4
1 10.2 5.4 9.5 6.3
2 6.2 9.0 7.9 10.3
3 12.8 7.7 10.0 9.9
4 10.8 6.6 8.2 10.6
5 11.6 8.6 10.0 15.6
, , sigma
1 2 3 4
1 6.3 7.7 7.5 4.8
2 7.3 6.0 9.9 9.8
3 6.8 6.9 6.7 9.4
4 9.7 7.7 8.9 9.5
5 10.8 6.5 8.1 7.0
# ... with 995 more iterations
3. Full High-Level Modeling with brms
The best part: Andrew’s stanli backend branch for brms works right out of the box in webR:
install.packages("rstan", repos = c("https://stan-dev.r-universe.dev", "https://repo.r-wasm.org/"))
# brms fork built for wasm
install.packages("brms", repos = c("https://staffanbetner.r-universe.dev", "https://repo.r-wasm.org/"))
library(brms)
fit <- brm(
formula = time ~ age * sex,
data = kidney,
backend = "stanli",
seed = 1234
)
# similar sampling output as above
summary(fit)
Result:
Family: gaussian
Links: mu = identity
Formula: time ~ age * sex
Data: kidney (Number of observations: 76)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 24.67 86.89 -147.60 196.37 1.00 1604 1985
age 0.68 1.91 -3.09 4.39 1.00 1585 1804
sexfemale 163.66 101.59 -35.78 357.08 1.00 1572 1978
age:sexfemale -2.43 2.21 -6.79 1.96 1.00 1510 1833
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 128.73 10.80 109.72 151.57 1.00 2404 2594
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
- Current friction points / Next steps:
Finding a smoother way to distributelibstanli.sofor the standalonestanlipackage in webR without hitting CORS restrictions (bundling the runtime directly in its Wasm build, similar to howstanrdoes it).- Chains run sequentially in single-threaded WebAssembly unless parallelized over Web Workers.
- 32-bit wasm memory limits for larger models.
Massive credit to @seantalts for developing stanli and @andrjohns for stanr and the brms integration!
Disclaimer: English is not my first language. The code, experiments, and findings are my own, but I used an LLM to assist with formatting and structuring the post for clarity.