Transforming posteriors with coords: derived quantities, ranks, binning and re-indexing

E01-E06 treat the posterior as something you summarise. Here it is something you compute with. Almost every question a reader asks about a climate model - when is the hottest day of the year? has it moved? which year was the warmest? how often will it top 32 °C now? - is a function of the parameters, and the honest answer is that function evaluated on every draw, then summarised. The draws already carry names (chain, draw, station, k, part, year); the whole trick of this notebook is to give every other axis you invent along the way (day of year, decade, threshold, return period, rank) a name too, and then let xarray line everything up. No axis=2, no reshape, no loops.

The model is deliberately modest - a seasonal cycle, a trend, a random effect per year and a skewed daily scatter, fitted in about 20 seconds. Most of the notebook is what you do with it.

#1 · From a long table to a labelled Dataset

GHCN-Daily ships one file per station, one row per (date, element) - TMAX, TMIN, PRCP, snow, wind... - with values in tenths of a unit and a quality flag (blank = passed NOAA's checks). The shortest route to a labelled array is: index the table by the two things that will become dimensions, to_xarray(), and turn the element dimension into separate variables with to_dataset(dim="element").

Alignment is where labelled data bites. Sea-Tac starts in January 1948, Portland in April 1938, and both files run into the current, unfinished year. Every binary operation between them aligns first, and the default for arithmetic is an inner join: the decade of Portland data before 1948 disappears without a word. xr.align lets you choose, and join="exact" is the defensive setting when you believe the axes already match.

To analyse the stations together we stack them along a new named dimension. Passing a pd.Index with a name to xr.concat creates the station dimension and its labels in one go; join="outer" keeps every day of both and NaN-fills the gaps, which we then cut to the common complete years. (dtype=object matters under pandas 3: a list of strings becomes a pandas StringDtype index, which xarray keeps, and a later concat with ordinary NumPy string labels then fails with "Cannot interpret StringDtype as a data type".) Derived calendar quantities go on as non-dimension coordinates along time: they travel with every selection, and groupby can use them.

A daily series 28,490 steps long hides its structure. Re-indexing time by the pair (year, doy) and unstacking turns it into a year x day-of-year grid - the same numbers, now arranged so that the seasonal cycle runs across and the decades run up. Everything in this pipeline is by name: a 7-day centred rolling mean, the re-index, the facet by station.

Portland (right) is hotter in summer and its summer is longer than Sea-Tac's (left) - the width of the red band. Individual hot spells show up as short dark-red streaks, and the band is visibly wider and darker in the last decade or two at both stations. The grid is not quite a calendar: in leap years every day after 28 February moves one column right, and only leap years fill column 366. At a 7-day smoothing that shift is invisible, but it is the reason we model on doy with a 365.25-day period rather than on columns of this grid.

How independent are consecutive days? Anomalies from the day-of-year average are the classic xarray idiom: groupby("time.dayofyear") and subtract the grouped mean - xarray broadcasts each group's value back onto its days. xr.corr with shift then gives the autocorrelation at any lag, for both stations at once.

Tomorrow is correlated 0.67 with today, but a week apart the correlation is down to about 0.12. A model that treats days as independent would count one heat wave as seven pieces of evidence and report intervals far too narrow. The simplest honest fix for a notebook whose subject is not autocorrelation: fit on every 7th day (4,070 days per station) and keep the other six days of every week as held-out data for checking the model.

#2 · A modest model, written on the same labels

For station $s$, year $y$ and day of year $d$ with $\omega = 2\pi d/365.25$ and time $t$ in decades from 1987:

$$ \begin{aligned} \text{TMAX} &\sim \text{SkewNormal}(\text{mean}=m,\ \text{sd}=\sigma,\ \text{skew}=\alpha)\\ m &= a_s + b_s t + u_{s,y} + \sum_{k=1}^{3}\big(c_{s,k}^{\cos}\cos k\omega + c_{s,k}^{\sin}\sin k\omega\big) + t\,\big(d_s^{\cos}\cos\omega + d_s^{\sin}\sin\omega\big)\\ \log\sigma &= l_s + g_s^{\cos}\cos\omega + g_s^{\sin}\sin\omega, \qquad \alpha = e_s + h_s^{\cos}\cos\omega + h_s^{\sin}\sin\omega, \qquad u_{s,y}\sim N(0,\tau_s) \end{aligned} $$

Three harmonics give the seasonal cycle, $b$ is the warming per decade, $d$ lets the seasonal cycle change over time (summers may warm faster than winters), $u$ is a warm or cool year, and the daily scatter is skewed with a skewness that changes through the year: winter highs have a long cold tail (arctic outbreaks), spring and summer highs a long hot tail (heat waves).

The basis is a DataArray. The harmonics are written once, as a function of a labelled doy, with named k and part dimensions matching the model's coords. The model uses its values on the observed rows; later, the same function evaluated on any grid of days broadcasts against the posterior draws by name.

One modelling detail matters for sampling. PyMC's SkewNormal(mu, sigma, alpha) uses the direct parametrisation, where mu is a location that moves whenever the skewness moves. Written that way, this model sampled into different modes per chain (r_hat 2.4 on the intercept and the skewness harmonics). Parametrising by the mean and sd - converting to the direct parameters inside the model - removes the trade-off, and the mean is also the quantity we want to talk about. The same function builds a plain Normal version, which we will need in section 5.

The priors allow daily maxima from well below freezing to the low 40s °C - wider than the data, nowhere absurd. The prior on $b$ (sd 0.5 °C per decade) is generous: the observed warming is a few tenths of a degree per decade.

Both fits are clean under nutpie: r_hat at most 1.01, ESS at least about 1,500 for every parameter, no divergences, and under half a minute each. az.rhat on a group returns a DataTree node; the .to_dataset().to_dataarray() detour turns "every variable" into one array whose .max() is the worst r_hat anywhere.

The headline parameter is already readable: both airports warmed by about 0.23 °C per decade ($b$), about 1.8 °C over the 78 years. The skewness $\alpha$ is positive on average ($e \approx 1$) with a strong negative cosine term ($h^{\cos}$): negative skew in winter, positive in summer - just what the data suggested.

#3 · Names, not positions

The posterior of c has five dimensions. Here is the same number twice.

values[:, :, 0, 0, 0] silently means something else as soon as any step reorders the dimensions - and it would not even raise. From here on everything is selected, reduced and broadcast by name. We pool chains and draws into one sample dimension with az.extract (500 random draws are plenty for the summaries below and keep every array small).

Two label details. az.extract indexes sample by a (chain, draw) MultiIndex; two random extracts - say from two different models - carry different labels, and combining them later would outer-join 250 + 250 samples into 500 half-NaN ones (xarray currently warns that this default will become join="exact"). Draws from different fits are not paired anyway, so we replace the MultiIndex by a plain counter. And the station labels are re-attached as object strings so they match ds exactly.

#4 · Derived quantities are functions of draws

A derived quantity is any function of the parameters. Write it once, as an ordinary function of labelled arrays, and it works on one draw, on 500 draws, on any grid of days or years you pass in - xarray broadcasts by dimension name. The model's expected TMAX (the climatology - the average high for a given day in a given year's climate, before that particular year's anomaly) and the daily spread:

Passing a decade array and a doy array produced a (station, sample, decade, doy) array: 2 x 500 x 8 x 365, no loop. The warming per decade on each day of the year is another derived quantity - $b + d^{\cos}\cos\omega + d^{\sin}\sin\omega$ - and a .quantile over sample gives a new quantile dimension to plot bands from.

On the scale of a 15 °C seasonal cycle, the 1.8 °C of warming is a modest but clearly resolved gap between the dashed (1950s) and solid (2020s) curves: the 90% bands of the expected high are only a few tenths of a degree wide. The right panel shows what the curves hide: warming is not uniform through the year. At both stations the posterior median is highest in mid-summer (about 0.28-0.30 °C per decade) and lowest in mid-winter (about 0.17).

The pointwise bands of the two seasons overlap, which tempts the conclusion "no evidence of a difference". That reads the wrong quantity. The contrast summer rate minus winter rate is itself a derived quantity - a difference of two .sels, draw by draw - and its posterior is what answers the question:

The summer-minus-winter difference in warming rate has a posterior median of about 0.13 °C per decade at Sea-Tac and 0.10 at Portland, and it is positive in 98% and 92% of draws. Overlapping pointwise bands are not evidence of "no difference": the two rates share $b$ within each draw, and the draw-by-draw difference cancels that shared uncertainty.

#A posterior over a date

idxmax("doy") returns the label of the maximum along a named dimension, so on the climatology it gives the date of the hottest day of the average year - once per draw, decade and station, i.e. a posterior distribution over a date. On a whole-day grid that posterior is lumpy (it can only take integer values); .interp to a 0.1-day grid near the peak (cubic, by name) smooths it. Because the climatology here is an explicit function we could also evaluate it on the fine grid directly - and the answers agree to the grid step.

The same method answers a question with more movement: when does summer arrive? Take "summer" as the part of the year when the average high is at least 21 °C (70 °F). (clim >= 21.1) is a boolean array; idxmax on it returns the first day that is True (the onset), and the length of summer is just .sum("doy").

The peak of the average year sits on about 4-5 August at Sea-Tac and 7 August at Portland, and it has not moved (the decades differ by about a day, within a posterior sd of about 1.3 days). What has moved is everything around it: in the 1950s climate the average high first reached 21 °C around 22 June at Sea-Tac, in the 2020s climate around 3 June; the autumn end moved from about 10 to 21 September, and the summer grew from about 80 to 110 days. Portland's summer is longer at both ends and grew similarly.

A ridgeline makes the shift visible: one row per decade, the posterior density of the onset, peak and end dates on a calendar axis. Every row is date_post.sel(station=..., decade=...) - the draws of one cell of a labelled array.

At both stations the onset ridges (left) walk steadily earlier decade by decade and the end ridges (right) walk later, while the paler peak ridges in the middle stay on the same August date. The ridges are narrow - a few days - because they are the posterior of the average climate, not of any one year's weather. The onset moves further than the end because the seasonal curve is flatter in June than in September, so the same warming shifts the crossing of 21 °C by more days.

The warming itself, two ways. $b$ is the model's underlying trend. The latent annual level $a + b\,t + u_y$ also includes each year's anomaly, and .polyfit("year", 1) - run on a (station, sample, year) array - fits a line through every draw's 78 annual levels at once: the realised 1948-2025 warming, with its own posterior. .polyfit on the raw annual means is the no-model answer.

The three estimates agree to within 0.02 °C per decade: roughly 0.23 at both stations, with 90% intervals of about ±0.06. The two stations are indistinguishable (the posterior probability that Portland warmed faster is near a coin flip) - the difference of two labelled arrays selected by station is itself a posterior.

#5 · Exceedance probabilities, return periods and a check on held-out days

How likely is a hot day (TMAX ≥ 32 °C) on a given date in a given year? For each draw that is the upper tail of that day's skew-normal. Two xarray points here:

• scipy.stats.skewnorm.sf is not a NumPy ufunc, so calling it on DataArrays would return a bare NumPy array (labels gone). xr.apply_ufunc wraps it: inputs broadcast by name, the result keeps every dimension and coordinate.
• This time we include each year's anomaly $u$: post["u"] already has a year dimension, so adding it to a climatology evaluated on post["year"] lines up year by year - no index arithmetic.

We restrict to May-September and 250 draws to keep the array at 48 MB (station x sample x doy x year = 2 x 250 x 153 x 78).

Setting the skewness to zero turns the same function into the Normal model's answer: we just added two all-zero variables to the Normal posterior with xr.zeros_like, which copies the right dims and coords.

The posterior mean of P(hot day) on every (day, year) is a picture of the changing climate. Over it we put the actual hot days from the full daily record - most of which the model never saw (six days out of seven were held out). .where(..., drop=True) pulls them out as a short list.

(The colour scale is square-root so that Sea-Tac's small probabilities stay visible.) The probability surface is a ridge centred on late July to early August that darkens toward the present and flickers year to year with $u$ - hot years such as 1958 and 2015 are darker columns. The observed hot days scatter over the dark part of each station's ridge. Portland's ridge is much wider and darker: 32 °C is a common summer day there and a rare one at Sea-Tac.

A held-out check, by counting. Summing over doy gives the expected number of hot days per year; drawing one Bernoulli per day and summing gives the predictive count, which is what an observed count should fall inside. We do this for both models and put them side by side with xr.concat over a new model dimension, then compare with the counts in the full daily record.

This is the check that separates the two likelihoods. With a Normal daily scatter (grey) the model misses the hot tail: at Sea-Tac it expects about half the hot days actually seen (1.2 against 2.6 a year before 1986), and only 71% (Sea-Tac) and 60% (Portland) of years fall inside its 90% interval. The skew-normal (purple) gets the long-run rate about right up to 2015 at both stations (2.3 against 2.6 and 3.6 against 3.5 a year at Sea-Tac; 10.2 against 9.9 and 13.3 against 14.2 at Portland), but two failures remain that no amount of skewness fixes:

• only 79% and 68% of years fall inside its 90% interval: the observed counts are more variable than predicted, with misses on both sides and more above. Real hot days come in heat waves; our daily draws are independent, so the predictive count is too well-behaved;
• since about 2015 both stations have had many more hot days than even the skew-normal expects (7.2 against 4.7 a year at Sea-Tac, 22.1 against 16.6 at Portland): recent summers have warmed faster than a straight-line trend.

That is useful to know before quoting any number below: the model's daily tail is reasonable on average, its heat-wave clustering is absent, and extremes are the place it is least trustworthy. (E37 models the hottest day of the year directly with extreme-value theory, which is the right tool for record-breaking events.)

Return periods of a 35 °C day. Assuming independent days, the probability that a year has at least one day at or above 35 °C is $1-\prod_d (1-p_d)$ - a product over the doy dimension, computed in log space - and its reciprocal is the return period in years. We evaluate it at the mid-points of the two halves of the record (1967 and 2006) and at 2025, and set the observed return period of each half beside it: the number of years divided by the number of years that had such a day (groupby("year").any(), then a two-bin BinGrouper).

The change is about right and the level is not. In the data a 35 °C day came once every 3.2 years at Sea-Tac in 1948-1986 and every 1.9 years in 1987-2025; the model says 2.2 and 1.6 years at the mid-points of those periods (and 1.4 for the 2025 climate), so it makes such a year too common. This is the other face of the clustering above: the model treats ~150 warm-season days as independent chances, but hot days arrive together in a few heat waves, so a year has fewer independent chances than days. Clustering leaves the expected number of hot days alone but lowers the chance of at least one. At Portland a 35 °C day happens in nearly every year, in the data and in the model.

#6 · Ranks inside each draw: which year was the warmest?

"Rank the posterior means" answers the wrong question. The posterior probability that year $y$ was the warmest is the fraction of draws in which $y$ has the highest latent annual level. So we rank inside each draw, along year.

DataArray.rank would do it, but it needs the optional bottleneck package, which this environment does not have. xr.apply_ufunc with input_core_dims=[["year"]] lifts scipy.stats.rankdata instead: xarray moves year to the last axis, the function ranks along it for every (station, sample), and the result comes back labelled. idxmax gives the same "warmest year" answer without ranks.

At Sea-Tac 2015 is the warmest year with a posterior probability of about 0.53 (its raw annual mean is the highest too), 2016 comes next at 0.14 and the rest spreads thinly over other recent years. Portland's ranking is more open: 2015 at 0.27, then 2021, 2025 and 2024 at 0.13-0.16 each. The posterior order can differ from the raw one - Sea-Tac's 2014 is second by raw annual mean but gets only 0.04, less than 2020 (ninth by raw mean): a year's level is estimated from the days the model saw (one in seven), so years with close raw means share probability in ways the raw order cannot show.

#Warming stripes that show their uncertainty

Ed Hawkins' "warming stripes" colour one vertical stripe per year by its temperature anomaly. A posterior has a distribution per year, so here each stripe is split vertically into 20 cells coloured by the 2.5%, ..., 97.5% posterior quantiles of the level anomaly (relative to 1961-1990). A year that is well determined is a uniform stripe; an uncertain one fades from cool at the bottom to warm at the top. The quantile dimension produced by .quantile(...) is the vertical axis. The raw stripe for comparison sits above.

Both stations show the familiar blue-to-red progression, and the posterior stripes add two things. The colour changes within most stripes from bottom to top: the 95% interval of one year's level is about 1.4 °C wide on average, as large as many year-to-year differences in the middle decades - which is why the "warmest year" probabilities spread over several years. And the levels are shrunk toward the trend line, so the posterior stripes are calmer than the raw ones (compare Sea-Tac's 1958: deep red raw, lighter in the posterior median).

Why is a year uncertain when its annual mean is known exactly? The level is the year's underlying warmth with day-to-day weather removed, and the model saw one day in seven - roughly the number of effectively independent days a year contains, given the autocorrelation of section 1.

#7 · Anomalies, resampling and binning on the observed record

With a posterior climatology in hand, anomalies of every observed day are one line: groupby("time.dayofyear") minus a climatology indexed by dayofyear. We use the posterior-mean climatology of the 1961-1990 climate (the WMO reference period, evaluated at its mid-point) on days 1-366, built on a dimension named dayofyear so that it matches the group name.

The anomalies average about zero over the reference period, as they should. Now the time machinery, all on the labelled record:

• BinGrouper bins years into custom, labelled periods, and passing two groupers to groupby gives a (period x season) table in one call. Multi-key groupby needs real coordinates: a virtual "time.season" works as a single key, but as a keyword key it raises KeyError - which is why season was attached as a coord in section 1.
• Seasons come back in alphabetical order (DJF, JJA, MAM, SON): .sel them into the calendar order.

The data agree with the model's seasonal warming curve and add a warning. Every season has warmed, and the largest recent anomalies are in summer (JJA): +1.9 °C at Sea-Tac and +2.6 °C at Portland in 2020-25, against +0.6 to +1.5 °C in the other seasons. A straight-line trend with a gently changing seasonal cycle (our $b$ and $d$) cannot follow a jump that size in the last few summers - the same thing the hot-day check found.

#A climate spiral

Monthly anomalies come from .resample(time="MS").mean(). Single months at one station are noisy (a monthly sd of about 1.5 °C), so we average each calendar month within each period - the same BinGrouper as above, now crossed with month - and draw one closed loop per period in polar coordinates: angle = month, radius = anomaly, colour = period. A warming climate is a set of loops that drift outward.

The 1948-59 loop (darkest) runs inside or on the zero ring in almost every month; the 2020-25 loop (yellow) runs outside it in most months and bulges furthest in June-August, 2-3 °C above the 1961-1990 climate in July. The loops in between cross each other a lot: the drift across decades is steady, but month by month it is buried in natural variability, and this display shows both at once. The recent summer bulge is larger than the change in any other season, as in the table above.

#8 · Lifting a SciPy root-finder over draws with apply_ufunc

The T-year return level is the temperature $x$ that the hottest day of the year exceeds with probability $1/T$. Under the model (independent days, year effect $u \sim N(0,\tau)$)

$$ P(\max_d \text{TMAX}_d \le x) = \mathbb{E}_u \prod_d F_d(x - u) = 1 - 1/T, $$

which has no closed form: solve it with scipy.optimize.brentq, integrating over $u$ with a 7-point Gauss-Hermite rule. brentq takes scalars, so the solver below is written for one draw, station, year and return period, with the day-of-year curves as 1-D arrays. xr.apply_ufunc(..., input_core_dims=[["doy"], ["doy"], ["doy"], [], []], vectorize=True) does the rest: doy is consumed, everything else (station, sample, year, period) is broadcast and looped over, and the answer comes back labelled.

Speed matters inside a vectorised loop: with stats.skewnorm.logcdf in the objective this cell took about 150 s; the same skew-normal CDF written with the ufuncs special.ndtr and special.owens_t (checked equal) takes a few seconds.

The model's curves move up by about 1 °C from the 1967 to the 2006 climate (and another 0.5 °C to 2025) at every return period. But at both stations the observed annual maxima of each half of the record (circles) lie below the matching model curve beyond a return period of about two years, and further below the rarer the event: the model's hottest day of the year is too hot. Two assumptions we made knowingly cause it - independent days (clustered hot days give fewer independent chances, as in section 5) and a skew-normal tail that keeps going where real summer highs level off (E37's extreme-value fits find a finite upper limit). The consequence is worth saying plainly: taken at face value this model calls Sea-Tac's 42.2 °C of 2021 roughly a 1-in-40-year event in the 2006 climate, while E37, modelling annual maxima directly, finds it far outside what the earlier record supports. A daily model is fine for the bulk and for counts; extremes need an extreme-value model. Here the point was the mechanics: a scalar SciPy solver lifted over 2 stations x 150 draws x 3 years x 10 periods = 9,000 solves by one apply_ufunc call.

#9 · Summary: the moves in this notebook

#Try it yourself

1. Nights, not days. Refit with TMIN instead of TMAX (one word in rows). Is the warming per decade larger at night, is the seasonal pattern of warming the same, and does the first frost-free date (climatological TMIN crossing 0 °C in spring, via the same idxmax-on-a-mask trick) move earlier?
2. A better trend. Replace $b\,t$ by a smooth curve (for example a random walk over 5-year blocks, indexed by a block coord) and redo the hot-day check of section 5. Does the 2016-2025 excess disappear? Which of the derived quantities change?
3. Heat-wave clustering. Fit on every day of June-August with an AR(1) daily anomaly (or simulate an AR(1) with the observed lag-1 correlation of 0.67 in the predictive draws) and recompute the predictive hot-day counts. Does the 90% interval now cover the observed years?