PHPstats

Modern statistical computing for PHP

A collection of native PHP packages and extensions for modern statistical workflows.

Load, transform, and summarize data from your app using data frames. Work with probability distributions, fit statistical models, and draw beautiful charts without leaving PHP.

PHPstats is in public beta. Feedback and bug reports are welcome.

Explore the packages Get started

Summarize and chart data Query a CSV file with Polars Fit a logistic regression Work with probability distributions

use Phibble\DataFrame;use function Phplyr\{col, floor_date, group_by};use function Phplyr\{mean, mutate, summarize}; $temps = DataFrame::fromCsv('daily-temperatures.csv');$monthly = $temps |>    mutate(month: floor_date(col('date'), 'month')) |>    group_by('month', 'city') |>    summarize(temp_c: mean('temp_c')); $chart = \Phiz\Plot::of($monthly)    ->aes(x: 'month', y: 'temp_c', color: 'city')    ->line()    ->scaleX(format: '%b')    ->labs(x: 'Month', y: 'Mean temperature (°C)', color: 'City'); echo $monthly;echo $chart->embed();

use Phibble\DataFrame;use function Phplyr\{col, floor_date, mean, phplyr}; $temps = DataFrame::fromCsv('daily-temperatures.csv');$monthly = phplyr($temps)    ->mutate(month: floor_date(col('date'), 'month'))    ->groupBy('month', 'city')    ->summarize(temp_c: mean('temp_c'))    ->get(); $chart = \Phiz\Plot::of($monthly)    ->aes(x: 'month', y: 'temp_c', color: 'city')    ->line()    ->scaleX(format: '%b')    ->labs(x: 'Month', y: 'Mean temperature (°C)', color: 'City'); echo $monthly;echo $chart->embed();

36 × 3

# Groups:   month [12]
   month      city     temp_c
   <date>     <chr>     <dbl>
 1 2025-01-01 London     4.09
 2 2025-01-01 New York   1.04
 3 2025-01-01 Sydney    22.4 
 4 2025-02-01 London     5.45
 5 2025-02-01 New York   1.68
 6 2025-02-01 Sydney    22.3 
 7 2025-03-01 London     7.58
 8 2025-03-01 New York   6.10
 9 2025-03-01 Sydney    20.4 
10 2025-04-01 London    10.3 
# ℹ 26 more rows

use function Phpolars\{col, scan_csv}; $quarter = col('date')->dt()->quarter()->alias('quarter'); $quarters = scan_csv('daily-temperatures.csv', tryParseDates: true)    ->groupBy(['city', $quarter])    ->agg(        col('temp_c')->min()->alias('coldest'),        col('temp_c')->mean()->round(1)->alias('mean'),        col('temp_c')->max()->alias('warmest'),    )    ->sort(['city', 'quarter'])    ->collect(); echo $quarters;

shape: (12, 5)
┌──────────┬─────────┬─────────┬──────┬─────────┐
│ city     ┆ quarter ┆ coldest ┆ mean ┆ warmest │
│ ---      ┆ ---     ┆ ---     ┆ ---  ┆ ---     │
│ str      ┆ i8      ┆ f64     ┆ f64  ┆ f64     │
╞══════════╪═════════╪═════════╪══════╪═════════╡
│ London   ┆ 1       ┆ 0.3     ┆ 5.7  ┆ 13.5    │
│ London   ┆ 2       ┆ 5.5     ┆ 14.3 ┆ 20.7    │
│ London   ┆ 3       ┆ 11.2    ┆ 17.4 ┆ 23.3    │
│ London   ┆ 4       ┆ 1.2     ┆ 8.9  ┆ 16.8    │
│ New York ┆ 1       ┆ -5.9    ┆ 3.0  ┆ 11.5    │
│ …        ┆ …       ┆ …       ┆ …    ┆ …       │
│ New York ┆ 4       ┆ -1.6    ┆ 8.2  ┆ 19.3    │
│ Sydney   ┆ 1       ┆ 17.0    ┆ 21.7 ┆ 29.3    │
│ Sydney   ┆ 2       ┆ 9.8     ┆ 16.7 ┆ 24.3    │
│ Sydney   ┆ 3       ┆ 8.9     ┆ 14.8 ┆ 20.7    │
│ Sydney   ┆ 4       ┆ 13.9    ┆ 20.4 ┆ 26.5    │
└──────────┴─────────┴─────────┴──────┴─────────┘

use App\Models\Order;use Phibble\DataFrame;use function Phit\fit_generalized_linear; $orders = DataFrame::fromQuery(Order::query());$fit = fit_generalized_linear('returned ~ channel + items',    $orders, family: 'binomial'); echo $fit;

Generalized linear model fit by maximum likelihood
Family: binomial (logit)
Formula: returned ~ channel + items

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)
(Intercept)        -3.25258    0.20925 -15.544  < 2e-16
channelMarketplace  1.26342    0.21756   5.807 6.35e-09
channelWeb          0.43824    0.20916   2.095   0.0362
items               0.12544    0.05557   2.257   0.0240

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1328.3 on 2399 degrees of freedom
Residual deviance: 1285.7 on 2396 degrees of freedom
AIC: 1293.7

use Phibble\DataFrame;use function Phit\dist_gamma_pdf_array; $x = range(0.0, 20.0, 0.25);$rows = [];foreach ([1.0, 2.0, 4.0] as $k) {    $pdf = dist_gamma_pdf_array($x, $k, scale: 2.0);    foreach ($x as $i => $v) {        $rows[] = ['x' => $v, 'k' => "$k", 'density' => $pdf[$i]];    }}$gamma = DataFrame::fromRows($rows); echo \Phiz\Plot::of($gamma)    ->aes(x: 'x', y: 'density', color: 'k')    ->line()    ->labs(x: 'x', y: 'Density', color: 'Shape')    ->embed();

Designed to work together

PHPstats is a family of interoperable PHP packages for data processing, statistical analysis, modeling, and visualization. Inspired by R's tidyverse, they share data structures and API idioms for a cohesive statistical computing experience in PHP.

Read the table with phpolars

$schools = \Phpolars\scan_csv('schools.csv')    ->select(['school', 'study_hours', 'score'])    ->collect(); // A data frame prints as a table.echo $schools->head(5);

First five rows

shape: (5, 3)
┌───────────┬─────────────┬────────┐
│ school    ┆ study_hours ┆ score  │
│ ---       ┆ ---         ┆ ---    │
│ str       ┆ f64         ┆ f64    │
╞═══════════╪═════════════╪════════╡
│ school-01 ┆ 11.58       ┆ 78.257 │
│ school-01 ┆ 2.7         ┆ 67.242 │
│ school-01 ┆ 10.02       ┆ 88.223 │
│ school-01 ┆ 3.71        ┆ 68.546 │
│ school-01 ┆ 4.38        ┆ 60.137 │
└───────────┴─────────────┴────────┘

Fit a mixed model with phit

$data = $schools->toColumns(); $fit = \Phit\fit_mixed_gaussian(    'score ~ study_hours + (1 | school)',    $data,); // A fit prints an R-style summary.echo $fit; $data['fitted'] = $fit['fitted'];

Linear mixed model fit by REML
Formula: score ~ study_hours + (1 | school)

REML criterion at convergence: 1583.3

Random effects:
Groups   Name        Variance Std.Dev.
school   (Intercept)    20.87    4.568
Residual                29.09    5.394
Number of obs: 250, groups: school, 12

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  56.4929     1.5408   36.66   <2e-16
study_hours   2.3428     0.1089   21.52   <2e-16

Chart the fitted lines with phiz

// Map columns to aesthetics and add layers, as in ggplot2.$chart = \Phiz\Plot::of($data)    ->aes(x: 'study_hours', y: 'score', color: 'school')    ->point(alpha: 0.4)    ->line(aes: ['y' => 'fitted'], alpha: 0.8)    ->guides(color: 'none')    ->labs(x: 'Study hours', y: 'Test score'); // Embeds the chart in the page.echo $chart->embed();

Test scores and fitted school lines

Test scores for 250 students in 12 schools, loaded from a CSV with phpolars, modeled with a Gaussian mixed model in phit, and visualized with phiz.

View the complete example

Packages

Choose a package

Load CSV files and database queries into typed tables you can inspect, filter, sort, and join. Similar to R's tibble and readr.

Filter, group, summarize, join, and reshape tables as a pipeline of steps. Similar to R's dplyr, tidyr, and janitor.

Process large tables quickly in the Polars engine, from CSV, Parquet, or SQL, through a PHP extension. Similar to Python's Polars.

Run PHP in Jupyter notebooks and Quarto documents, with tables and charts shown inline. Similar to R's IRkernel.

Fit linear, generalized linear, and mixed models, and compute probabilities and summary statistics. Similar to R's stats and lme4.

Build charts layer by layer from your tables, and save or embed them as Vega-Lite. Similar to R's ggplot2.

Fit Bayesian models written in Stan and summarize their posterior draws. Its API is modeled on R's CmdStanR.

Data & Workflow

Statistics, Modeling, & Visualization

Example: Load an Eloquent query with typed columns

use App\Models\Order;use Phibble\DataFrame; $orders = DataFrame::fromQuery(Order::query(), types: [    'channel' => 'factor',    'coupon' => 'bool',    'returned' => 'bool',]); echo $orders->glimpse();

Example output

Data frame glimpse

Rows: 2,400
Columns: 8
$ id        <int> 79201, 79202, 79203, 79204, 79205, 79206, 79207, 79208, 7920…
$ placed_on <date> 2025-01-11, 2025-01-13, 2025-01-15, 2025-01-18, 2025-01-19,…
$ channel   <fct> Web, App, App, Web, Web, Web, App, Web, Marketplace, Web, We…
$ customer  <chr> "Returning", "Returning", "New", "Returning", "Returning", "…
$ items     <int> 1, 1, 2, 1, 3, 1, 1, 3, 1, 5, 2, 2, 1, 6, 5, 1, 1, 1, 2, 2, …
$ coupon    <lgl> FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE,…
$ total     <dbl> 42.19, 46.79, 71.76, 40.65, 79.93, 31.87, 36.50, 69.93, 33.7…
$ returned  <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALS…

View full example

Example: Filter, group, and summarize orders

use App\Models\Order; use function Phplyr\{col, filter, group_by, mean, n, round, sum, summarize}; $orders = \Phibble\DataFrame::fromQuery(Order::query()); $byChannel = $orders |>    filter(col('returned')->eq(0)) |>    group_by('channel') |>    summarize(        orders: n(),        revenue: sum('total'),        average: round(mean('total'), 2),    ); echo $byChannel;

use App\Models\Order; use function Phplyr\{col, mean, n, phplyr, round, sum}; $orders = \Phibble\DataFrame::fromQuery(Order::query()); $byChannel = phplyr($orders)    ->filter(col('returned')->eq(0))    ->groupBy('channel')    ->summarize(        orders: n(),        revenue: sum('total'),        average: round(mean('total'), 2),    )    ->get(); echo $byChannel;

Example output

Data frame preview · 3 × 4

  channel     orders revenue average
  <chr>        <int>   <dbl>   <dbl>
1 App            689  40039.    58.1
2 Marketplace    420  19346.    46.1
3 Web           1101  57111.    51.9

View full example

Example: Aggregate database rows with expressions

use Illuminate\Support\Facades\DB; use function Phpolars\{col, len, read_database}; $orders = read_database(    'select channel, cast(total as double precision) as total     from orders where not returned',    DB::connection()->getPdo(),); $byChannel = $orders    ->groupBy('channel')    ->agg(        len()->alias('orders'),        col('total')->sum()->round(2)->alias('revenue'),    )    ->sort('revenue', descending: true); echo $byChannel;

Example output

Polars table

shape: (3, 3)
┌─────────────┬────────┬──────────┐
│ channel     ┆ orders ┆ revenue  │
│ ---         ┆ ---    ┆ ---      │
│ str         ┆ u32    ┆ f64      │
╞═════════════╪════════╪══════════╡
│ Web         ┆ 1101   ┆ 57111.08 │
│ App         ┆ 689    ┆ 40039.21 │
│ Marketplace ┆ 420    ┆ 19346.19 │
└─────────────┴────────┴──────────┘

View full example

Example: A PHP cell in a Quarto document

---title: "Revenue by channel"jupyter: php--- ```{php}#| echo: falseuse function Phplyr\{arrange, desc, group_by, mean, n, round, sum, summarize};``` ```{php}$orders = \Phibble\DataFrame::fromRows(\App\Data\Orders::generate()); $orders |>    group_by('channel') |>    summarize(        orders: n(),        revenue: round(sum('total'), 2),        items: round(mean('items'), 2),    ) |>    arrange(desc('revenue'));```

Example output

Rendered notebook cell · 3 × 4

	channel <chr>	orders <int>	revenue <dbl>	items <dbl>
1	Web	1185	61735.9	1.83
2	App	723	42157.12	1.85
3	Marketplace	492	22817.18	1.74

View full example

Example: Fit a linear model

use App\Models\Track; $tracks = \Phibble\DataFrame::fromQuery(Track::query()); $fit = \Phit\fit_ols('streams ~ playlist_adds + label', $tracks); echo \Phibble\DataFrame::fromColumns([    'term' => $fit['columns'],    'estimate' => $fit['coefficients'],    'std_error' => $fit['std_errors'],    'p_value' => $fit['p_values'],]);

Example output

Coefficients · 3 × 4

  term             estimate std_error   p_value
  <chr>               <dbl>     <dbl>     <dbl>
1 (Intercept)         -2.20    0.655  8.64e-  4
2 playlist_adds        1.07    0.0210 2.53e-194
3 labelMajor label     8.07    0.755  5.05e- 24

View full example

Example: Points with a fitted line per group

use App\Models\Track; $tracks = \Phibble\DataFrame::fromQuery(Track::query()); $chart = \Phiz\Plot::of($tracks)    ->aes(x: 'playlist_adds', y: 'streams', color: 'label')    ->point(alpha: 0.4)    ->smooth(method: 'lm')    ->labs(        x: 'Playlist adds (thousands)',        y: 'Streams (millions)',        color: null,    ); echo $chart->embed();

Example output

Chart

View full example

Example: Sample a hierarchical model and summarize it

use App\Models\StudentScore;use CmdStan\StanData; use function CmdStan\cmdstan_model; $scores = \Phibble\DataFrame::fromQuery(StudentScore::query()); $data = StanData::fromFormula('score ~ study_hours + (1 | school)', $scores); $fit = cmdstan_model('school_intercepts.stan')    ->sample(data: $data, chains: 4, thin: 4, seed: 20261016); echo $fit->summaryTable(    variables: ['beta', 'sigma', 'tau'],    measures: ['mean', 'sd', 'quantile2', 'rhat', 'ess_bulk'],);

Example output

Posterior summary · 4 × 7

  variable  mean    sd    q5   q95  rhat ess_bulk
  <chr>    <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>
1 beta[1]  56.5  1.67  53.7  59.1   1.01     481.
2 beta[2]   2.34 0.108  2.16  2.52  1.00     991.
3 sigma     5.41 0.243  5.03  5.84  1.00     946.
4 tau       5.03 1.25   3.37  7.43  1.00     752.

View full example

Verified against established libraries

use function Phit\{    dist_gamma_cdf,    dist_gamma_quantile,    dist_normal_cdf,    dist_normal_quantile,}; $values = [    dist_normal_cdf(-8.0),    dist_normal_quantile(0.975),    dist_gamma_cdf(12.0, shape: 2.0, scale: 3.0),    dist_gamma_quantile(0.999, shape: 2.0, scale: 3.0),];

Normal and gamma probabilities and quantiles

Normal CDF at −8: phit6.22 × 10^-16 R 6.22 × 10^-16 SciPy (Python) 6.22 × 10^-16
Normal quantile at 0.975: phit1.96 R 1.96 SciPy (Python) 1.96
Gamma CDF at 12: phit0.908 R 0.908 SciPy (Python) 0.908
Gamma quantile at 0.999: phit27.7 R 27.7 SciPy (Python) 27.7

Distributions: phit's values and each library's values
Quantity	phit	R	SciPy (Python)
Normal CDF at −8	6.22 × 10^-16	6.22 × 10^-16	6.22 × 10^-16
Normal quantile at 0.975	1.96	1.96	1.96
Gamma CDF at 12	0.908	0.908	0.908
Gamma quantile at 0.999	27.7	27.7	27.7

Standard normal, and gamma with shape 2 and scale 3.

R: pnorm, qnorm, pgamma, and qgamma.

Python: SciPy's norm and gamma.

Distributions verification record

use Phibble\DataFrame;use function Phibble\Expr\{col, mean, n}; $data = 'https://raw.githubusercontent.com/phpstats/datasets/main/csv';$weather = DataFrame::fromCsv("$data/weather.csv"); $byAirport = $weather    ->groupBy(['origin'])    ->summarize([        'hours' => n(),        'temp' => mean(col('temp'), na_rm: true),        'humid' => mean(col('humid'), na_rm: true),    ]);

use function Phpolars\{col, len, read_csv}; $data = 'https://raw.githubusercontent.com/phpstats/datasets/main/csv';$weather = read_csv(    "$data/weather.csv",    columns: ['origin', 'temp', 'humid'],); $byAirport = $weather    ->groupBy('origin')    ->agg(        len()->alias('hours'),        col('temp')->mean()->alias('temp'),        col('humid')->mean()->alias('humid'),    )    ->sort('origin');

CSV processing of airport weather

EWR hours: phibble8,703 phpolars 8,703 dplyr (R) 8,703 pandas (Python) 8,703
EWR mean temperature: phibble55.55 phpolars 55.55 dplyr (R) 55.55 pandas (Python) 55.55
EWR mean humidity: phibble63.06 phpolars 63.06 dplyr (R) 63.06 pandas (Python) 63.06
JFK hours: phibble8,706 phpolars 8,706 dplyr (R) 8,706 pandas (Python) 8,706
JFK mean temperature: phibble54.47 phpolars 54.47 dplyr (R) 54.47 pandas (Python) 54.47
JFK mean humidity: phibble65.21 phpolars 65.21 dplyr (R) 65.21 pandas (Python) 65.21
LGA hours: phibble8,706 phpolars 8,706 dplyr (R) 8,706 pandas (Python) 8,706
LGA mean temperature: phibble55.76 phpolars 55.76 dplyr (R) 55.76 pandas (Python) 55.76
LGA mean humidity: phibble59.32 phpolars 59.32 dplyr (R) 59.32 pandas (Python) 59.32

CSV processing: phibble's values and each library's values
Quantity	phibble	phpolars	dplyr (R)	pandas (Python)
EWR hours	8,703	8,703	8,703	8,703
EWR mean temperature	55.55	55.55	55.55	55.55
EWR mean humidity	63.06	63.06	63.06	63.06
JFK hours	8,706	8,706	8,706	8,706
JFK mean temperature	54.47	54.47	54.47	54.47
JFK mean humidity	65.21	65.21	65.21	65.21
LGA hours	8,706	8,706	8,706	8,706
LGA mean temperature	55.76	55.76	55.76	55.76
LGA mean humidity	59.32	59.32	59.32	59.32

R: readr and dplyr.

Python: pandas.

Data: weather, 26,115 hourly readings at three New York City airports in 2013 (source)

use Phibble\DataFrame;use function Phit\fit_ols; $data = 'https://raw.githubusercontent.com/phpstats/datasets/main/csv';$mtcars = DataFrame::fromCsv("$data/mtcars.csv"); $fit = fit_ols('mpg ~ wt + hp', $mtcars);

Linear regression of fuel economy

(Intercept): phit37.23 lm (R) 37.23 statsmodels (Python) 37.23
wt: phit−3.88 lm (R) −3.88 statsmodels (Python) −3.88
hp: phit−0.03 lm (R) −0.03 statsmodels (Python) −0.03
Residual SD: phit2.59 lm (R) 2.59 statsmodels (Python) 2.59
R²: phit0.83 lm (R) 0.83 statsmodels (Python) 0.83

OLS: phit's values and each library's values
Quantity	phit	lm (R)	statsmodels (Python)
(Intercept)	37.23	37.23	37.23
wt	−3.88	−3.88	−3.88
hp	−0.03	−0.03	−0.03
Residual SD	2.59	2.59	2.59
R²	0.83	0.83	0.83

Model: mpg ~ wt + hp

R: lm().

Python: ols().

Data: mtcars, 32 cars from Motor Trend, 1974 (source)

OLS regression verification record

use Phibble\DataFrame;use function Phit\fit_generalized_linear; $data = 'https://raw.githubusercontent.com/phpstats/datasets/main/csv';$mroz = DataFrame::fromCsv("$data/mroz.csv"); $fit = fit_generalized_linear(    'lfp ~ k5 + k618 + age + wc + hc + lwg + inc',    $mroz,    family: 'binomial',);

Logistic regression of labor force participation

(Intercept): phit3.18 glm (R) 3.18 statsmodels (Python) 3.18
k5: phit−1.46 glm (R) −1.46 statsmodels (Python) −1.46
k618: phit−0.06 glm (R) −0.06 statsmodels (Python) −0.06
age: phit−0.06 glm (R) −0.06 statsmodels (Python) −0.06
wcyes: phit0.81 glm (R) 0.81 statsmodels (Python) 0.81
hcyes: phit0.11 glm (R) 0.11 statsmodels (Python) 0.11
lwg: phit0.60 glm (R) 0.60 statsmodels (Python) 0.60
inc: phit−0.03 glm (R) −0.03 statsmodels (Python) −0.03
Log-likelihood: phit−452.63 glm (R) −452.63 statsmodels (Python) −452.63

Logistic: phit's values and each library's values
Quantity	phit	glm (R)	statsmodels (Python)
(Intercept)	3.18	3.18	3.18
k5	−1.46	−1.46	−1.46
k618	−0.06	−0.06	−0.06
age	−0.06	−0.06	−0.06
wcyes	0.81	0.81	0.81
hcyes	0.11	0.11	0.11
lwg	0.60	0.60	0.60
inc	−0.03	−0.03	−0.03
Log-likelihood	−452.63	−452.63	−452.63

Model: lfp ~ k5 + k618 + age + wc + hc + lwg + inc

R: glm(family = binomial).

Python: logit().

Data: mroz, 753 married women in 1975 (source)

Logistic regression verification record

use Phibble\DataFrame;use function Phit\fit_generalized_linear; $data = 'https://raw.githubusercontent.com/phpstats/datasets/main/csv';$mroz = DataFrame::fromCsv("$data/mroz.csv"); $fit = fit_generalized_linear(    'lfp ~ k5 + k618 + age + wc + hc + lwg + inc',    $mroz,    family: 'binomial',    link: 'probit',);

Probit regression of labor force participation

(Intercept): phit1.92 glm (R) 1.92 statsmodels (Python) 1.92
k5: phit−0.87 glm (R) −0.87 statsmodels (Python) −0.87
k618: phit−0.04 glm (R) −0.04 statsmodels (Python) −0.04
age: phit−0.04 glm (R) −0.04 statsmodels (Python) −0.04
wcyes: phit0.49 glm (R) 0.49 statsmodels (Python) 0.49
hcyes: phit0.06 glm (R) 0.06 statsmodels (Python) 0.06
lwg: phit0.37 glm (R) 0.37 statsmodels (Python) 0.37
inc: phit−0.02 glm (R) −0.02 statsmodels (Python) −0.02
Log-likelihood: phit−452.69 glm (R) −452.69 statsmodels (Python) −452.69

Probit: phit's values and each library's values
Quantity	phit	glm (R)	statsmodels (Python)
(Intercept)	1.92	1.92	1.92
k5	−0.87	−0.87	−0.87
k618	−0.04	−0.04	−0.04
age	−0.04	−0.04	−0.04
wcyes	0.49	0.49	0.49
hcyes	0.06	0.06	0.06
lwg	0.37	0.37	0.37
inc	−0.02	−0.02	−0.02
Log-likelihood	−452.69	−452.69	−452.69

Model: lfp ~ k5 + k618 + age + wc + hc + lwg + inc

R: glm(family = binomial("probit")).

Python: probit().

Data: mroz, 753 married women in 1975 (source)

Probit regression verification record

use Phibble\DataFrame;use function Phit\fit_mixed_gaussian; $data = 'https://raw.githubusercontent.com/phpstats/datasets/main/csv';$sleepstudy = DataFrame::fromCsv("$data/sleepstudy.csv"); $fit = fit_mixed_gaussian(    'Reaction ~ Days + (1 | Subject)',    $sleepstudy,);

Linear mixed model of reaction times

(Intercept): phit251.41 lme4 (R) 251.41 statsmodels (Python) 251.41
Days: phit10.47 lme4 (R) 10.47 statsmodels (Python) 10.47
Subject variance: phit1,378.18 lme4 (R) 1,378.18 statsmodels (Python) 1,378.18
Residual variance: phit960.46 lme4 (R) 960.46 statsmodels (Python) 960.46
Log-likelihood: phit−893.23 lme4 (R) −893.23 statsmodels (Python) −893.23

Mixed Gaussian: phit's values and each library's values
Quantity	phit	lme4 (R)	statsmodels (Python)
(Intercept)	251.41	251.41	251.41
Days	10.47	10.47	10.47
Subject variance	1,378.18	1,378.18	1,378.18
Residual variance	960.46	960.46	960.46
Log-likelihood	−893.23	−893.23	−893.23

Model: Reaction ~ Days + (1 | Subject)

R: lmer().

Python: mixedlm().

Data: sleepstudy, 18 subjects over 10 days of sleep restriction (source)

Gaussian mixed model verification record

use Phibble\DataFrame;use function Phit\fit_mixed_logistic; $data = 'https://raw.githubusercontent.com/phpstats/datasets/main/csv';$cbpp = DataFrame::fromCsv(    "$data/cbpp.csv",    types: ['period' => 'f'],); $fit = fit_mixed_logistic(    'cbind(incidence, size - incidence) ~ period + (1 | herd)',    $cbpp,);

Mixed logistic regression of pleuropneumonia

(Intercept): phit−1.40 lme4 (R) −1.40 GPBoost (Python) −1.40
period2: phit−0.99 lme4 (R) −0.99 GPBoost (Python) −0.99
period3: phit−1.13 lme4 (R) −1.13 GPBoost (Python) −1.13
period4: phit−1.58 lme4 (R) −1.58 GPBoost (Python) −1.58
Herd variance: phit0.41 lme4 (R) 0.41 GPBoost (Python) 0.41
Log-likelihood: phit−92.03 lme4 (R) −92.03 GPBoost (Python) −92.03

Mixed logistic: phit's values and each library's values
Quantity	phit	lme4 (R)	GPBoost (Python)
(Intercept)	−1.40	−1.40	−1.40
period2	−0.99	−0.99	−0.99
period3	−1.13	−1.13	−1.13
period4	−1.58	−1.58	−1.58
Herd variance	0.41	0.41	0.41
Log-likelihood	−92.03	−92.03	−92.03

Model: cbind(incidence, size - incidence) ~ period + (1 | herd)

R: glmer().

Python: GPBoost.

Data: cbpp, 15 cattle herds over four periods (source)

Mixed logistic model verification record

use function CmdStan\cmdstan_model; $model = cmdstan_model('schools.stan'); $fit = $model->sample(    data: [        'J' => 8,        'y' => [28, 8, -3, 7, -1, 1, 18, 12],        'sigma' => [15, 10, 16, 11, 9, 11, 10, 18],    ],    seed: 20261101,    chains: 4,    iter_warmup: 1000,    iter_sampling: 1000,);

Bayesian hierarchical model of coaching effects

mu: CmdStanPHP4.43 CmdStanR (R) 4.43 CmdStanPy (Python) 4.43
tau: CmdStanPHP3.57 CmdStanR (R) 3.57 CmdStanPy (Python) 3.57
theta[1]: CmdStanPHP6.13 CmdStanR (R) 6.13 CmdStanPy (Python) 6.13
theta[2]: CmdStanPHP4.94 CmdStanR (R) 4.94 CmdStanPy (Python) 4.94
theta[3]: CmdStanPHP4.02 CmdStanR (R) 4.02 CmdStanPy (Python) 4.02
theta[4]: CmdStanPHP4.76 CmdStanR (R) 4.76 CmdStanPy (Python) 4.76
theta[5]: CmdStanPHP3.66 CmdStanR (R) 3.66 CmdStanPy (Python) 3.66
theta[6]: CmdStanPHP4.03 CmdStanR (R) 4.03 CmdStanPy (Python) 4.03
theta[7]: CmdStanPHP6.39 CmdStanR (R) 6.39 CmdStanPy (Python) 6.39
theta[8]: CmdStanPHP4.98 CmdStanR (R) 4.98 CmdStanPy (Python) 4.98

Mixed Gaussian (Bayesian): CmdStanPHP's values and each library's values
Quantity	CmdStanPHP	CmdStanR (R)	CmdStanPy (Python)
mu	4.43	4.43	4.43
tau	3.57	3.57	3.57
theta[1]	6.13	6.13	6.13
theta[2]	4.94	4.94	4.94
theta[3]	4.02	4.02	4.02
theta[4]	4.76	4.76	4.76
theta[5]	3.66	3.66	3.66
theta[6]	4.03	4.03	4.03
theta[7]	6.39	6.39	6.39
theta[8]	4.98	4.98	4.98

Model: y ~ normal(mu + tau * theta, sigma)

R: CmdStanR.

Python: CmdStanPy.

All three run the same CmdStan executable, and their draws are identical.

Data: schools, eight schools, Rubin (1981) (source)

Bayesian interfaces verification record

Performant, powerful, and flexible

$normalCdf = dist_normal_cdf_array($x);$normalQuantile = dist_normal_quantile_array($p);$gammaCdf = dist_gamma_cdf_array($positive, 2.5, 1.5);$gammaQuantile = dist_gamma_quantile_array($p, 2.5, 1.5);

normal_cdf <- pnorm(x)normal_quantile <- qnorm(p)gamma_cdf <- pgamma(positive, shape = 2.5, scale = 1.5)gamma_quantile <- qgamma(p, shape = 2.5, scale = 1.5)

normal_cdf = stats.norm.cdf(x)normal_quantile = stats.norm.ppf(p)gamma_cdf = stats.gamma.cdf(positive, 2.5, scale=1.5)gamma_quantile = stats.gamma.ppf(p, 2.5, scale=1.5)

$normalCdf = dist_normal_cdf_array($x);$normalQuantile = dist_normal_quantile_array($p);$gammaCdf = dist_gamma_cdf_array($positive, 2.5, 1.5);$gammaQuantile = dist_gamma_quantile_array($p, 2.5, 1.5);

Medians of 7 runs, each in a fresh process.

2 AMD Zen 5 cores, 8 GB RAM, Linux, PHP 8.5.12, recorded 2026-11-01.

Memory is the peak anonymous memory: the most memory the process allocated at once, including the runtime, the data, and native allocations, but not file pages it maps.

The full record gives versions, thread settings, and where timing starts and stops.

View the complete benchmark

Normal CDF and quantile, gamma (shape 2.5, scale 1.5) CDF and quantile, 1,000,000 inputs each

Distributions: median elapsed time and peak anonymous memory of 7 runs
Implementation	Elapsed time	Peak anonymous memory
phit	850 ms	113 MiB
R	860 ms	102 MiB
SciPy (Python)	403 ms	191 MiB
phit-polyfill	19.4 s	114 MiB

For easier installation, use phit-polyfill. For faster computation, use the phit extension.

use function Phpolars\{col, len, read_csv}; $summary = read_csv(    $path,    columns: [        'placed_on', 'region', 'status', 'quantity', 'amount',    ],)    ->filter(col('status')->eq('completed'))    ->groupBy([        'region',        col('placed_on')->str()->slice(0, 7)->alias('month'),    ])    ->agg(        len()->alias('orders'),        col('quantity')->sum()->alias('units'),        col('amount')->sum()->alias('revenue'),        col('amount')->mean()->alias('mean_amount'),    )    ->sort(['region', 'month']); $rows = $summary->rows();

use Phibble\DataFrame;use function Phibble\Expr\{col, mean, n, str_sub, sum}; $summary = DataFrame::fromCsv(    $path,    types: ['placed_on' => 'c'],)    ->filter(col('status')->eq('completed'))    ->mutate(['month' => str_sub(col('placed_on'), 1, 7)])    ->groupBy(['region', 'month'])    ->summarize([        'orders' => n(),        'units' => sum(col('quantity')),        'revenue' => sum(col('amount'), na_rm: true),        'mean_amount' => mean(col('amount'), na_rm: true),    ])    ->arrange(['region', 'month']); $rows = array_map(array_values(...), $summary->toRows());

library(readr)library(dplyr) summary <- read_csv(  path,  col_select = c(    placed_on, region, status, quantity, amount  ),  col_types = cols(placed_on = col_character()),  progress = FALSE,  num_threads = 2) |>  filter(status == "completed") |>  group_by(region, month = substr(placed_on, 1, 7)) |>  summarise(    orders = n(),    units = sum(quantity),    revenue = sum(amount, na.rm = TRUE),    mean_amount = mean(amount, na.rm = TRUE),    .groups = "drop"  ) |>  arrange(region, month)

import pandas orders = pandas.read_csv(    path,    usecols=[        "placed_on", "region", "status", "quantity", "amount",    ],)completed = orders[orders["status"] == "completed"]summary = (    completed    .assign(month=completed["placed_on"].str.slice(0, 7))    .groupby(["region", "month"], sort=True)    .agg(        orders=("status", "size"),        units=("quantity", "sum"),        revenue=("amount", "sum"),        mean_amount=("amount", "mean"),    )    .reset_index())

Medians of 7 runs, each in a fresh process.

2 AMD Zen 5 cores, 8 GB RAM, Linux, PHP 8.5.12, recorded 2026-11-01.

Memory is the peak anonymous memory: the most memory the process allocated at once, including the runtime, the data, and native allocations, but not file pages it maps.

The full record gives versions, thread settings, and where timing starts and stops.

View the complete benchmark

Order-events CSV, 20,000,000 rows, 1 GB

CSV processing: median elapsed time and peak anonymous memory of 7 runs
Implementation	Elapsed time	Peak anonymous memory
phpolars	3.0 s	4.4 GiB
phibble (phit reader)	22.1 s	4.5 GiB
phibble (PHP reader)	33.3 s	4.8 GiB
tibble (R)	7.2 s	2.0 GiB
pandas (Python)	6.3 s	3.3 GiB

use function Phpolars\{col, len, scan_csv}; $summary = scan_csv($path)    ->filter(col('status')->eq('completed'))    ->groupBy([        'region',        col('placed_on')->str()->slice(0, 7)->alias('month'),    ])    ->agg(        len()->alias('orders'),        col('quantity')->sum()->alias('units'),        col('amount')->sum()->alias('revenue'),        col('amount')->mean()->alias('mean_amount'),    )    ->sort(['region', 'month'])    ->collect(engine: 'streaming'); $rows = $summary->rows();

use Phibble\DataFrame;use function Phibble\Expr\{col, is_na, n, str_sub, sum}; $totals = [];$add = function (DataFrame $chunk) use (&$totals) {    $part = $chunk        ->filter(col('status')->eq('completed'))        ->mutate([            'month' => str_sub(col('placed_on'), 1, 7),        ])        ->groupBy(['region', 'month'])        ->summarize([            'orders' => n(),            'units' => sum(col('quantity')),            'revenue' => sum(col('amount'), na_rm: true),            'amounts' => sum(is_na(col('amount'))->not()),        ]);    foreach ($part->toRows() as $r) {        $key = $r['region']."\0".$r['month'];        $t = $totals[$key] ?? [            'orders' => 0, 'units' => 0,            'revenue' => 0.0, 'amounts' => 0,        ];        foreach ($t as $field => $total) {            $t[$field] = $total + $r[$field];        }        $totals[$key] = $t;    }};DataFrame::fromCsvChunked(    $path,    $add,    chunkSize: 100_000,    types: ['placed_on' => 'c'],);ksort($totals, SORT_STRING);$rows = [];foreach ($totals as $key => $t) {    [$region, $month] = explode("\0", $key);    $mean = $t['amounts'] > 0        ? $t['revenue'] / $t['amounts']        : null;    $rows[] = [        $region, $month,        $t['orders'], $t['units'], $t['revenue'], $mean,    ];}

$in = fopen($path, 'r');$header = fgetcsv($in, escape: '');$groups = [];while (($row = fgetcsv($in, escape: '')) !== false) {    [, $placedOn, $region, , $status, $quantity, $amount]        = $row;    if ($status !== 'completed') {        continue;    }    $key = $region."\0".substr($placedOn, 0, 7);    $g = $groups[$key] ?? [        'orders' => 0, 'units' => 0,        'revenue' => 0.0, 'amounts' => 0,    ];    $g['orders']++;    $g['units'] += (int) $quantity;    if ($amount !== '') {        $g['revenue'] += (float) $amount;        $g['amounts']++;    }    $groups[$key] = $g;}fclose($in);ksort($groups, SORT_STRING);$rows = [];foreach ($groups as $key => $g) {    [$region, $month] = explode("\0", $key);    $mean = $g['amounts'] > 0        ? $g['revenue'] / $g['amounts']        : null;    $rows[] = [        $region, $month,        $g['orders'], $g['units'], $g['revenue'], $mean,    ];}

library(readr)library(dplyr) parts <- read_csv_chunked(  path,  DataFrameCallback$new(function(chunk, pos) {    chunk |>      filter(status == "completed") |>      group_by(region, month = substr(placed_on, 1, 7)) |>      summarise(        orders = n(),        units = sum(quantity),        revenue = sum(amount, na.rm = TRUE),        amounts = sum(!is.na(amount)),        .groups = "drop"      )  }),  chunk_size = 100000,  col_types = cols_only(    placed_on = col_character(),    region = col_character(),    status = col_character(),    quantity = col_integer(),    amount = col_double()  ),  progress = FALSE)summary <- parts |>  group_by(region, month) |>  summarise(    orders = sum(orders),    units = sum(units),    revenue = sum(revenue),    mean_amount = sum(revenue) / sum(amounts),    .groups = "drop"  ) |>  arrange(region, month)

import pandas parts = []chunks = pandas.read_csv(    path,    usecols=[        "placed_on", "region", "status", "quantity", "amount",    ],    chunksize=100_000,)for chunk in chunks:    completed = chunk[chunk["status"] == "completed"]    parts.append(        completed        .assign(            month=completed["placed_on"].str.slice(0, 7),            amounts=completed["amount"].notna(),        )        .groupby(["region", "month"])        .agg(            orders=("status", "size"),            units=("quantity", "sum"),            revenue=("amount", "sum"),            amounts=("amounts", "sum"),        )    )totals = (    pandas.concat(parts)    .groupby(level=["region", "month"])    .sum())summary = (    totals    .assign(mean_amount=totals["revenue"] / totals["amounts"])    .reset_index()    .sort_values(["region", "month"]))

Medians of 7 runs, each in a fresh process.

2 AMD Zen 5 cores, 8 GB RAM, Linux, PHP 8.5.12, recorded 2026-11-01.

Memory is the peak anonymous memory: the most memory the process allocated at once, including the runtime, the data, and native allocations, but not file pages it maps.

The full record gives versions, thread settings, and where timing starts and stops.

View the complete benchmark

Order-events CSV, 20,000,000 rows, 1 GB

CSV processing: median elapsed time and peak anonymous memory of 7 runs
Implementation	Elapsed time	Peak anonymous memory
phpolars	1.9 s	40 MiB
phibble (phit reader)	17.3 s	71 MiB
phibble (PHP reader)	29.5 s	59 MiB
plain PHP	38.8 s	7 MiB
tibble (R)	7.7 s	161 MiB
pandas (Python)	6.9 s	102 MiB

$fit = fit_ols(    'y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10',    $data,);

fit <- lm(  y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10,  data = d)table <- coef(summary(fit))

fit = smf.ols(    'y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10',    data=d,).fit()coefficients, standard_errors = fit.params, fit.bse

$fit = fit_ols(    'y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10',    $data,);

Medians of 7 runs, each in a fresh process.

2 AMD Zen 5 cores, 8 GB RAM, Linux, PHP 8.5.12, recorded 2026-11-01.

Memory is the peak anonymous memory: the most memory the process allocated at once, including the runtime, the data, and native allocations, but not file pages it maps.

The full record gives versions, thread settings, and where timing starts and stops.

View the complete benchmark

OLS, 1,000,000 rows, 10 predictors

OLS: median elapsed time and peak anonymous memory of 7 runs
Implementation	Elapsed time	Peak anonymous memory
phit	742 ms	799 MiB
lm (R)	1.3 s	1.2 GiB
statsmodels (Python)	357 ms	568 MiB
phit-polyfill	5.2 s	1.0 GiB

For easier installation, use phit-polyfill. For faster computation, use the phit extension.

$fit = fit_generalized_linear(    'y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10',    $data,    family: 'binomial',);

fit <- glm(  y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10,  data = d,  family = binomial)table <- coef(summary(fit))

fit = smf.logit(    'y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10',    data=d,).fit(method='newton', disp=0)coefficients, standard_errors = fit.params, fit.bse

$fit = fit_generalized_linear(    'y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10',    $data,    family: 'binomial',);

Medians of 7 runs, each in a fresh process.

2 AMD Zen 5 cores, 8 GB RAM, Linux, PHP 8.5.12, recorded 2026-11-01.

Memory is the peak anonymous memory: the most memory the process allocated at once, including the runtime, the data, and native allocations, but not file pages it maps.

The full record gives versions, thread settings, and where timing starts and stops.

View the complete benchmark

Logistic regression, 1,000,000 rows, 10 predictors, 0/1 response

Logistic: median elapsed time and peak anonymous memory of 7 runs
Implementation	Elapsed time	Peak anonymous memory
phit	1.7 s	801 MiB
glm (R)	2.5 s	1.2 GiB
statsmodels (Python)	495 ms	485 MiB
phit-polyfill	20.8 s	1.4 GiB

For easier installation, use phit-polyfill. For faster computation, use the phit extension.

$fit = fit_generalized_linear(    'y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10',    $data,    family: 'binomial',    link: 'probit',);

fit <- glm(  y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10,  data = d,  family = binomial("probit"))table <- coef(summary(fit))

fit = smf.probit(    'y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10',    data=d,).fit(method='newton', disp=0)coefficients, standard_errors = fit.params, fit.bse

$fit = fit_generalized_linear(    'y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10',    $data,    family: 'binomial',    link: 'probit',);

Medians of 7 runs, each in a fresh process.

2 AMD Zen 5 cores, 8 GB RAM, Linux, PHP 8.5.12, recorded 2026-11-01.

Memory is the peak anonymous memory: the most memory the process allocated at once, including the runtime, the data, and native allocations, but not file pages it maps.

The full record gives versions, thread settings, and where timing starts and stops.

View the complete benchmark

Probit regression, 1,000,000 rows, 10 predictors, 0/1 response

Probit: median elapsed time and peak anonymous memory of 7 runs
Implementation	Elapsed time	Peak anonymous memory
phit	3.0 s	801 MiB
glm (R)	2.7 s	1.2 GiB
statsmodels (Python)	635 ms	485 MiB
phit-polyfill	43.3 s	1.4 GiB

For easier installation, use phit-polyfill. For faster computation, use the phit extension.

$fit = fit_mixed_gaussian(    'y ~ x1 + x2 + x3 + (1 | g)',    $data,);

fit <- lmer(  y ~ x1 + x2 + x3 + (1 | g),  data = d,  REML = TRUE)table <- coef(summary(fit))

fit = smf.mixedlm(    'y ~ x1 + x2 + x3',    d,    groups=d['g'],).fit(reml=True)coefficients = fit.fe_paramsstandard_errors = fit.bse_fe

$fit = fit_mixed_gaussian(    'y ~ x1 + x2 + x3 + (1 | g)',    $data,);

Medians of 7 runs, each in a fresh process.

2 AMD Zen 5 cores, 8 GB RAM, Linux, PHP 8.5.12, recorded 2026-11-01.

Memory is the peak anonymous memory: the most memory the process allocated at once, including the runtime, the data, and native allocations, but not file pages it maps.

The full record gives versions, thread settings, and where timing starts and stops.

View the complete benchmark

Gaussian random intercept by REML, 20,000 rows in 200 groups, 3 predictors

Mixed Gaussian: median elapsed time and peak anonymous memory of 7 runs
Implementation	Elapsed time	Peak anonymous memory
phit	15 ms	16 MiB
lme4 (R)	82 ms	244 MiB
statsmodels (Python)	143 ms	106 MiB
phit-polyfill	104 ms	25 MiB

For easier installation, use phit-polyfill. For faster computation, use the phit extension.

$fit = fit_mixed_logistic(    'y ~ x1 + x2 + x3 + (1 | g)',    $data,);

fit <- glmer(  y ~ x1 + x2 + x3 + (1 | g),  data = d,  family = binomial,  control = glmerControl(tolPwrss = 1e-10))table <- coef(summary(fit))

model = gpb.GPModel(    group_data=d['g'],    likelihood='bernoulli_logit',)model.set_optim_params(params={    'optimizer_cov': 'lbfgs',    'delta_rel_conv': 1e-12,})model.fit(y=y, X=X)table = model.get_coef(std_err=True)

$fit = fit_mixed_logistic(    'y ~ x1 + x2 + x3 + (1 | g)',    $data,);

Medians of 7 runs, each in a fresh process.

2 AMD Zen 5 cores, 8 GB RAM, Linux, PHP 8.5.12, recorded 2026-11-01.

Memory is the peak anonymous memory: the most memory the process allocated at once, including the runtime, the data, and native allocations, but not file pages it maps.

The full record gives versions, thread settings, and where timing starts and stops.

View the complete benchmark

Binomial logit random intercept by Laplace, 20,000 rows in 200 groups, 3 predictors, 0/1 response

Mixed logistic: median elapsed time and peak anonymous memory of 7 runs
Implementation	Elapsed time	Peak anonymous memory
phit	298 ms	18 MiB
lme4 (R)	1.7 s	254 MiB
GPBoost (Python)	89 ms	111 MiB
phit-polyfill	3.0 s	38 MiB

For easier installation, use phit-polyfill. For faster computation, use the phit extension.

$model = cmdstan_model(exeFile: $exe);$fit = $model->sample(    data: $data,    seed: 20261101,    chains: 4,    parallel_chains: 4,    iter_warmup: 1000,    iter_sampling: 1000,    output_dir: $out,    show_messages: false,);$summary = $fit->summary(    ['mu', 'tau', 'theta'],    measures: [        'mean', 'sd', 'quantile2', 'rhat', 'ess_bulk',    ],);

model <- cmdstan_model(exe_file = exe)fit <- model$sample(  data = data,  seed = 20261101,  chains = 4,  parallel_chains = 4,  iter_warmup = 1000,  iter_sampling = 1000,  output_dir = out,  show_messages = FALSE,  refresh = 0)summary <- fit$summary(  c("mu", "tau", "theta"),  "mean",  "sd",  ~quantile(.x, probs = c(0.05, 0.95)),  "rhat",  "ess_bulk")

model = cmdstanpy.CmdStanModel(exe_file=exe)fit = model.sample(    data=data,    seed=20261101,    chains=4,    parallel_chains=4,    iter_warmup=1000,    iter_sampling=1000,    output_dir=out,    show_progress=False,)summary = fit.summary(percentiles=(5, 95))

Medians of 7 runs, each in a fresh process.

2 AMD Zen 5 cores, 8 GB RAM, Linux, PHP 8.5.12, recorded 2026-11-01.

Memory is the peak anonymous memory: the most memory the process allocated at once, including the runtime, the data, and native allocations, but not file pages it maps.

The full record gives versions, thread settings, and where timing starts and stops.

View the complete benchmark

Eight schools, 4 chains × 2,000 iterations

Mixed Gaussian (Bayesian): median elapsed time and peak anonymous memory of 7 runs
Implementation	Elapsed time	Peak anonymous memory
CmdStanPHP	133 ms	11 MiB
CmdStanR (R)	446 ms	145 MiB
CmdStanPy (Python)	113 ms	53 MiB

Install the packages Browse examples