Blokjutsu: a macOS image generator/editor

Posted on November 15, 2026 4:09 PM by Bob Carpenter

This is another post from Bob that is not about statistics. Last one, I promise.

Like Steve, whose what’s up in NYC app I just blogged about, I’ve taken up the hobby of designing apps and having LLMs build them. It scratches a very different itch than writing performant C++ code for Stan.

Blokjutsu

Here’s what I’ve come up with.

Blokjutsu: An open-source native macOS app for generating and editing images locally.

Fair warning: it requires Apple Silicon with at least 32GB RAM. On the plus side, it’s distributed as an Apple-notarized .dmg file and licensed entirely under Apache 2, including the models.

MOSS-Image, the 5 billion parameter model on which it is based, is amazing in its ability to follow text instructions. Having said that, don’t expect it to compete with a frontier model in depth or breadth of knowledge. It also generates images up to 4 megapixels, though that can be slow (20 minutes on my old Mac with an M3 Ultra, 12 on my new MacBook with an M5 Max—they’ve accelerated matrix multiplies!).

What does it do?

The app lets you either work entirely from a prompt, or you can structure your input with bounding boxes that are followed very well for either images or text (though the text is handled with fonts and doesn’t integrate real well with photographic images in my opinion). It uses a language model (Qwen) to process prompts and it can be used to enhance prompts directly or to convert images to text. It includes all the usual image-processing features: image to image, variations of generated images, inpainting (replacing part of an image) and outpainting (extending an image). It lets you plug in LoRAs (low rank adaptations) for styling.

I tried very hard to keep everything simple. If anyone has any feedback on either the doc or the system itself, I’d love to hear it, either via email ([email protected]) or an issue on GitHub.

What does “Blokjutsu” mean?

I made up a word that’s pronounceable, has zero Google hits, and isn’t spell corrected (as Hadley Wickham and I were urging instead of “Stan” for the Stan project). It’s a mashup of Danish “Blok” and Japanese “jutsu”. Same for the GitHub org, which is Danish “Form” mashed with Japanese “do”.

Why?

I built this app in part to see what Claude (Opus 5.5, medium effort) could do, and partly because I subscribe to Claude, which doesn’t do image generation! I like to use images in my other hobby, role-playing games. I was using DiffusionBee with StableDiffusion, but DiffusionBee hasn’t been updated in years, and I wanted to be able to use MOSS-Image. The next section has an example.

The Temple of Shadow

My childhood friends were getting tired of my obsession with fiction-first games (like Blades in the Dark) and wanted a “beer and pretzels” (as this style is known) dungeon crawl. I ran “Caverns of Thracia,” which the grognards describe as among the handful of best modules ever written, but we all found it challenging, and quit after the second TPK (total party kill). Hence the 5th edition module, where it’s a real challenge to kill even one character.

Dungeon Module Z1: The Temple of Shadow

The module text was entirely written by ChatGPT, though I provided a lot of editorial oversight, including a suggestion of how to arrange the four broad zones of the dungeon and asking continually for more concise descriptions with more showing and less telling. Everything else, including the serviceable but rather boring cookie-cutter town, was pure GPT. The LaTeX styling was me with GPT’s help. The maps were random generators available online. The images were Stable Diffusion after a lot of metaphorical hair pulling (literal is challenging for me). It’s all a cliché, just as my players requested. The writing is the usual AI type slop and the images are pretty atrocious. If you look at the writing from Wizards of the Coast, I’d call that “WotC slop,” so the bar is low. It’s hard to find good writing in an RPG, so it really stands out reading, for example, a short module by Will Doyle. I’ve been meaning to give this another go in late 2026 and see how Claude does, but I really don’t like reading AI writing.

What I like

Left to my own devices, I run games like I’m running now, Star Trek: Omega, where the Federation is down to its last Class-D starship pulled out of mothballs, with a skeleton crew of officers and enlisted folks who somehow missed the decisive battle in which the Federation was defeated. Captain K’vok (Klingon) and his motley crew are up against the Chorus, an alien computer intelligence harvesting all the planets in the galaxy for compute centers to achieve the singularity. Their five year mission is to make sure the last embers of humanity don’t burn out. I’m reluctant to link it as I don’t want to get a take-down notice, but you can find the Forged-in-the-Dark rules for it on my own GitHub org, bob-carpenter, as well as the rules for Ready, Steady, Go!, my Fast-and-Furious knockoff set in the 1968 Paris student riots, where the English team races against French, Italian, and German competition (all with their iconic sports cars) for the racing title, Le Roi de Paris. The real prize, of course, is much much bigger, and may be uncovered for the first time since the war. If you’re seriously thinking about playing, I can send along the maps and car tokens I used for Roll20. As a taste, here’s the token I generated for one of my players. I could not for the life of me get GPT to render it with right-side drive—I’m hoping MOSS-Image will be better.

The other player’s character drove an Aston Martin Zagato in British racing green, but of course they could pull out the Land Rover for the chase scene through Père Lachaise.

2026 election conspiracies haven’t breached the Vermeule threshold.

Posted on November 15, 2026 9:17 AM by Andrew

Longtime readers of this blog will recall the concept of the John Yoo line: The point at which nothing you write gets taken seriously, and so you might as well become a hack because you have no scholarly reputation remaining.

Speaking of legal scholars who’ve gone off the reservation, I was recently thinking about a related concept, the Vermeule threshold: The level at which an outrageous conspiracy theory becomes acceptable for respectable suit-and-tie types to promote.

This idea is named after the Harvard law professor Adrian Vermeule, who notoriously promoted stupid 2020 election denial material:

This takes on a particularly chilling note given that, just a few years earlier, Vermeule had written:

Because those who hold conspiracy theories typically suffer from a “crippled epistemology,” in accordance with which it is rational to hold such theories, the best response consists in cognitive infiltration of extremist groups.

In case you’re wondering, by “cognitive infiltration,” he explicitly states that he wants the government to be doing this infiltration. So, yeah, this is a very well-connected guy who’s pushing the government to infiltrate political groups—when he’s not promoting stupid-ass conspiracies himself.

I was thinking about this recently because someone in comments pointed to this risible bit of election denial from 2026:

I went online and I couldn’t see any promotion by Vermeule of this or any other 2026 election conspiracy theories. That’s why I characterize these theories as below the Vermeule threshold: so stupid that even the distinguished legal scholar cannot bring himself to endorse them.

Do I think Vermeule was sincere in his promotion of election denial in 2020? No, I doubt it. He’s well educated and, I assume, no idiot. Rather, my guess is that he thought that by circulating such theories he was being a charming imp, offering a good poke in the eye to complacent liberals etc. What a charming rogue, etc. Seeing video of the 6 Jan 2021 insurrection and the government shooting protesters on the street and then lying about it, I consider elite promotion of stupid conspiracy theories to be dangerous, not charming—but maybe part of this is a different perspective coming from a political scientist and a lawyer. I’m trying to understand reality here, with all the difficulties that entails; Vermeule is an attorney and a teacher of law, and that comes with a tradition of advocacy.

So, just to clarify, I don’t think that “below the Vermeule threshold” represents theories that are too stupid to for him to believe; rather, they are too stupid for him to publicly promote. These are theories that are so dumb that the Harvard professor would take more of a reputational hit by endorsing, than he would gain in political points by committing to.

Can you imagine being a law student and taking a class from this guy? This isn’t a case of a researcher with a side interest in extreme political views, along the lines of a J. B. S. Haldane or William Shockley, where you can value the science and engineering and set the politics aside. As with John Dewey, say, Vermeule’s political activism is close to his intellectual work. To be studying from someone who impishly promotes stupid and dangerous conspiracy theories . . . that’s a kind of instrumental view of the law that doesn’t make me happy. I guess that many of Vermeule’s students share his political views, and perhaps the others value the experience of seeing, from close up, the ideology of some of our political rulers. And the Harvard law students surely have the opportunity of experiencing some of this from the other side from their faculty who espouse extreme left-wing views.

P.S. In comments I laid out the argument more carefully:

There is a logic to these spurious election-denial conspiracy theories. It goes like this:

1. The election did not go the way I hoped and expected.

2. Therefore, something wrong happened.

3. I can’t accept that the “something wrong” could be that the majority of the voters legitimately supported the opposition.

4. Therefore, some version of cheating happened.

The next step is for people to make up evidence supporting the claim of cheating, or to suspend one’s critical faculties so as to accept spurious evidence. That’s what Vermeule was doing in 2020, except that I’m guessing he was one step “meta” from that: rather than suspending his critical faculties so as to believe stupid election-denial claims, he thinks he’s being a charming imp by promoting claims that he’s pretty sure are false, as a poke in the eye to the sanctimonious liberals who he sees as a danger to society. The “meta” bit is that I think he views liberals and Democrats as routinely cheating, misrepresenting the arguments of conservatives, abusing government power, etc., so it’s only fair play for this to be done on the other side. Whether or not Vermeule personally comes up with fake election denial theories, he appreciates the efforts of those on his side who do so, in the same sense that a “clean” hockey player respects the absolutely necessary role played by the “enforcers” on his team. No unilateral disarmament, right?

The awkward thing for Vermeule is that he can’t lay out the above logic. When he promotes stupid theories that he surely doesn’t believe, he can’t publicly say that the theories are stupid, or that for political reasons hadn’t looked into them so that he can maintain a plausible deniability about them. He has to play it straight, to follow the rules of the game.

From an intellectual sense, there’s a sort of tragedy here. A smart person who has devoted his intellectual life to the pursuit of logic for political goals, is driven by the logic of politics to endorse ridiculous claims, and his curse is that he cannot allow himself to explain it. He’s very much in the position of those left-wing western intellectuals from the Cold War period who endorsed Stalinism in all its aspects and could never let the mask slip, sacrificing their reputation and their intellectual integrity.

If I weren’t so bothered by the whole insurrection thing, I might feel kind of sorry for the guy.

If the air traffic control scandal had happened two weeks ago . . .

Posted on November 14, 2026 5:40 PM by Andrew

What we said here.

Zombie citations: A short empirical investigation

Posted on November 14, 2026 9:04 AM by Andrew

So, we had this discussion the other day about nonexistent or irrelevant citations included in articles produced by chatbots. This is part of a larger concern about AI slop polluting the scholarly literature.

The usual way this is framed is that you look at the reference list of published article or set of articles and then you ask how many of these citations are erroneous in some way.

But maybe it’s easier to look at things the other way: Start with an article, see where it’s been cited, and follow up to see if these citations make sense.

I thought about this in the context of my classic article, published in 2010 on Arxiv and included in a book four years later, “How many zombies do you know?” Using indirect survey methods to measure alien attacks and outbreaks of the undead.

Look up this article on Google Scholar and you’ll see 6 citations. Only 6! That makes me sad. In any case, here they are:

I clicked through to read how my paper was references, and 5 of the 6 basically got it right. The sixth one messes up–it cites the zombies paper as a source for my claim that the probability is about 1 in 10 million that a vote will be decisive–I’m guessing there was some sloppy copying and pasting from my list of published papers.

So, from this short study, my estimate of the probability of a bad citation is 1/6. This was from the pre-AI era, so it’s not a “hallucination”; it’s just someone who was not careful in citation.

If the Jensen Huang / Nvidia spy scandal had happened two weeks ago . . .

Posted on November 13, 2026 6:10 PM by Andrew

If the Jensen Huang / Nvidia spy scandal had happened two weeks ago, I expect the midterm elections would’ve gone pretty much just as they actually went, but lots of pundits would’ve seized on this as the explanation for the election outcomes.

Why do I say this?

Because the result of the election was surprising. Not in its direction—everyone was expecting the Democrats to win a majority of the vote—but in the extent of the majority (in the event, they won 57% of the two-party vote for the House of Representatives) and the consequences for both houses of the legislature.

When there’s a surprise, people look for an explanation. That’s what stories are all about. As the most convenient recent politically-related news item, the Huang scandal would’ve fit the bill.

I’d still guess that, had this all happened fifteen or thirty years ago, the Huang-equivalents would’ve been Democrats. But, as it is, Huang is closely connected to an unpopular Republican president, so the story would’ve written itself. Even though I think that, had the Huang/Huang thing come up two weeks ago, the election outcome would’ve been essentially unaffected.

When is a class a “commodity” and when does it depend on the instructor?

Posted on November 13, 2026 9:20 AM by Andrew

Sometimes I get emails from students asking to work on a project with me, and I usually respond that they should take one of my classes first. If they try to work on a project without having taken a class first, I usually end up having to teach them a bunch of things. Sure, they could read my books, but I guess people don’t really like to read, also there’s something extra you get from hearing it live, doing the homework, getting feedback in class, etc.

This always puzzled me, that students would think they could jump in a project without taking any of my classes first . . . but then I realized that it’s usual to think of a math class as being pretty much the same no matter who teaches it. You take calculus or linear algebra or real analysis or whatever, and it can be taught well or poorly, there are good teachers and bad teachers, but they’re pretty much all teaching the same thing.

But statistics is different. Even if you limit yourself to a narrow topic–for example, regression analysis or survey sampling or Bayesian statistics (just to name three fairly generic-sounding classes that I’ve taught)–the course varies a huge amount from instructor to instructor. Not just teaching style (cos, yeah, my style is much different from that of every other statistics teacher I’ve seen. I’m not saying it’s better, mind you, just that it’s different. You can get a sense of it by reading Part 1 of Active Statistics), but actually what gets taught. First off, most people don’t teach regression out of Regression and Other Stories, and other regression textbooks don’t look much like ours. Second, somehow what gets covered in class isn’t just what’s in the book. This is even true with Bayesian statistics. Lots of people teach out of Bayesian Data Analysis–and I hope that, soon, lots of people will be teaching from Bayesian Workflow–but there are still so many different directions to go.

And I’m not saying I’m unique here. In general, if you want to work with statistician X, it’s probably a good idea to first take a course from X, for just about any X.

In that way, a statistics class is less like a math class or a physics class or a chemistry class and more like a history class or a literature class. Even if the books are set ahead of time, different teachers teach it so differently that they’re just different courses entirely.

But I could see how many students, whose experience is taking straight math courses in high school and college, and maybe some very empty course such as AP statistics or the college equivalent, would look at the roster of statistics classes at a university and just see a collection of topics rather than what it really is, a collection of topics interacted with instructors. And our course descriptions usually don’t do such a good job of conveying that either. Hence this post.

P.S. I’m not saying that a course being a standardized commodity is a bad thing. In a way, that what my colleagues and I have been trying to do with our textbooks, as well as with Active Statistics: we’re trying to give other people the tools to teach as we do. But, the way things stand, it still matters a lot who teaches your class–maybe not such a bad thing, as the students should get some benefit from their big tuition payments. It just seems like a secret to many students that the instructor matters, not just in how well the topic is taught but also in the content of the course.

Just Bronder it: Steve’s bronder.nyc guide to NYC off(-off)-Broadway, jazz, readings, and galleries

Posted on November 12, 2026 3:00 PM by Bob Carpenter

This post is from Bob and it is not about statistics.

Steve Bronder, our primary C++ developer for Stan these days, has built a really neat web app (nothing to download/install) consolidating all the fun stuff happening in NYC. It’s absolutely free or I wouldn’t be advertising here.

bronder.nyc: Off-Broadway and off-off-Broadway theatre, readings, jazz, and galleries in NYC right now

The web app is beautifully simple and incredibly useful. Looking for jazz shows, Mitzi and I used to surf the individual venue pages we knew, but the waters are choppy. It’s incredibly slow and frustrating given how badly most of them are laid out for information discovery. Steve’s app just lists everything out in one scrollable list that’s very easy to narrow down to a particular day or interval. The map makes it very easy to find something in walking distance. I didn’t realize just how much I needed this app until he built it and I started using it.

I would highly recommend this if you’re visiting NYC and recommend it even more if you live here. The map mode is super useful for planning things like dates. I’m not just the shill, I’m also a client. Mitzi and I just went on a double date (!) with Nawaf and Christina to the Flatiron Room in Murray Hill this past Sunday night, which was awesome. I found The Flatiron Room on the app, having never heard of the venue before. The Moses Patrou Duo cranked out some really swinging and boppy soul jazz/funk on Hammond organ and drums. Loud enough to hear, but not so loud you couldn’t have a conversation. Putting in a plug for Flatiron Room, there’s no cover, and the Japanese food was really good and reasonably priced, which is about a hundred times better than it needs to be and half the price it could be in a venue like this (this is also true of my fave sushi in NYC, Studio 151, which is above Nublu in the East Village, and has the next tier up of sushi, but still at relatively bargain prices).

P.S. I was amused by Steve’s cookies message, which is a very Steve way of stating the issue: I use Google Analytics to see which shows and theaters are popular. No personal data, no ads. [Gosh I love cookies (Accept)] [Gosh I hate cookies (Reject)].

P.P.S. I am hoping “Bronder” becomes a common word like “Zagat,” but as a verb, like “Google.” I just learned about LLM poisoning attacks, so consider this document 1 of 250.

Bayesian Workflow (my talk in NYC, Wed 18 Nov 2:30pm)

Posted on November 12, 2026 9:40 AM by Andrew

This talk is sponsored by the New York City Chapter of the American Statistical Association, and it’s at CUNY Graduate Center, 365 Fifth Avenue, room 8301 & 8304 (8th Floor). If you want to attend (in person or remotely; I recommend in person!), just go to the form before 15 Nov and fill it out.

Also here’s a zoom link so if you’re attending remotely, maybe you don’t need to sign up.

Here’s the summary:

Bayesian statistics and statistical practice have evolved over the years, driven by advancements in theory, methods, and computational tools. This book explores the intricate workflows of applied Bayesian statistics, aiming to uncover the tacit knowledge often overlooked in published papers and textbooks. By systematizing the process of Bayesian model development, the book seeks to improve applied analyses and inspire future innovations in theory, methods, and software. It emphasizes the importance of iterative model building, model checking, computational troubleshooting, and simulated-data experimentation, offering a comprehensive perspective on statistical analysis.

Through detailed examples and practical guidance, the book bridges the gap between theory and application, empowering practitioners and researchers to navigate the complexities of Bayesian inference. It is not a checklist or cookbook but a flexible framework for understanding and resolving challenges in statistical modeling and decision making under uncertainty.

The key innovations of Bayesian Workflow include:

Going beyond Bayesian data analysis (model building, inference, model checking, and model expansion) to consider the larger process of statistical modeling, including comparisons of multiple models fit to a single dataset.

A fuller use of informative priors. One reason we’ve moved toward stronger priors is that the replication crisis has taught us that the amount of prior information available in any given problem is often approximately the same as the information coming from an experiment (see here, for example). Informative priors also fit our increased focus on generative modeling, and we’re doing a lot more prior predictive checking to understand the implications of our models.

More integration between modeling, data analysis, and computing. One way to see this is that the Bayesian Workflow webpage has the code to run all our examples.

Simulation-based experimentation. At its most basic, simulation-based experimentation provides a best-case baseline for statistical methods: if you can’t recover your quantities of interest with sufficient accuracy under ideal conditions (when your data are simulated from the model you’re fitting), then you know you’re in trouble. And often this is the case! Beyond that, we can simulate from one model and fit another, and see what happens. Simulation experiments aren’t always so easy to construct, as they involve specifying the entire data-generation process.

All the information on the book, including the book itself, is at the Bayesian Workflow webpage.

The non-Bayesians among you might be interested in Appendix A of the book, on Statistical and computational workflow for Bayesians and non-Bayesians, and also our short general article, Statistical Workflow, and its accompanying journal issue.

In the talk I’ll follow my usual practice and tell lots of stories. Stories are how we learn about the world, and I’ll try to also go meta and discuss the anomalous and immutable aspects of the stories I tell.

If you can make it to the talk, please come prepared with some tough questions!

Constructing a latent-variable model to fit measurement data: A polychoric story

Posted on November 11, 2026 6:48 PM by Andrew

Following up on our recent post, Mathias Berggren writes:

I coded up a polychoric estimate for the problem you present in Constructing a latent-variable model to fit measurement data: A Stan story, part 2. I think it led to pretty decent estimates. I attach the R-code in case you want to have a look.

I hadn’t thought about it, but I guess it makes sense that this sort of model is used enough that people would have already hard-coded it!

Here’s the R code:

# Polychoric solution of latent variable stan example from Andrew

# Load required libraries
library(polycor)

# Data as a matrix
data = matrix(c(208,501,51,640,9879,1439,40,464,1432), nrow = 3, ncol = 3)

# Observed joint probabilities
Pobs = data / sum(colSums(data))

# Polychoric estimate of the latent correlation
r = polychor(data)

# Observed marginal probabilities
pX = colSums(data) / sum(colSums(data))
pY = rowSums(data) / sum(colSums(data))

# Corresponding marginal thresholds
cX = c(qnorm(pX[1]),qnorm(pX[1]+pX[2]))
cY = c(qnorm(pY[1]),qnorm(pY[1]+pY[2]))

# Estimated joint probabilities from this information
P = matrix(nrow = 3, ncol = 3)
P[1,1] = integrate(function(y) dnorm(y) * pnorm((cX[1]-r*y)/sqrt(1-r^2)), -Inf, cY[1])$value
P[1,2] = integrate(function(y) dnorm(y) * (pnorm((cX[2]-r*y)/sqrt(1-r^2)) - pnorm((cX[1]-r*y)/sqrt(1-r^2))), -Inf, cY[1])$value
P[2,1] = integrate(function(y) dnorm(y) * pnorm((cX[1]-r*y)/sqrt(1-r^2)), cY[1], cY[2])$value
P[2,2] = integrate(function(y) dnorm(y) * (pnorm((cX[2]-r*y)/sqrt(1-r^2)) - pnorm((cX[1]-r*y)/sqrt(1-r^2))), cY[1], cY[2])$value
P[1,3] = pY[1] - P[1,1] - P[1,2]
P[2,3] = pY[2] - P[2,1] - P[2,2]
P[3,1] = pX[1] - P[1,1] - P[2,1]
P[3,2] = pX[2] - P[1,2] - P[2,2]
P[3,3] = pX[3] - P[1,3] - P[2,3]

# Compare estimated and observed probabilities
P
Pobs
Pobs - P

# Find least squares estimates to transform thresholds to 1.5 and 3.5
b1X = sd(c(1.5, 3.5)) / sd(cX)
b0X = mean(c(1.5, 3.5)) - b1X * mean(cX)

b1Y = sd(c(1.5, 3.5)) / sd(cY)
b0Y = mean(c(1.5, 3.5)) - b1Y * mean(cY)

# Check  so we get the correct values
b0X + b1X * cX
b0Y + b1Y * cY

# Least squares estimates
# These show how we should transform the standard normal marginal 
# distributions above to a scale with thresholds 1.5 and 3.5.
# Thus, they also indicate the updated mean and standard deviation of 
# each variable (as they transform from one with mean = 0 and sd = 1).
# If one wants one common mean, then one might take some weighted mean 
# of the intercepts
c(b0X,b1X)
c(b0Y,b1Y)

# Discussion
# If there are more thresholds then the least squares solutions may not be 
# exact, and we may find that the thresholds cannot be at 1.5, 2.5, etc.
# Same way, one might find that the intercept cannot be the same.

The challenge of truly probabilistic programming

Posted on November 11, 2026 9:16 AM by Andrew

Many years ago, Jouni and I wrote an unpublished paper called Fully Bayesian Computing:

A fully Bayesian computing environment calls for the possibility of defining vector and array objects that may contain both random and deterministic quantities, and syntax rules that allow treating these objects much like any variables or numeric arrays. Working within the statistical package R, we introduce a new object-oriented framework based on a new random variable data type that is implicitly represented by simulations.

A few years later we published an updated version of the paper under the title, Manipulating and summarizing posterior simulations using random variable objects. Here’s the background:

In the Bayesian paradigm, unknown quantities have probability distributions and are thus random variables. Observed values are just realizations of random variables, and constants may be thought of as random variables with point mass distributions. In mathematical notation, we deal with objects that are random variables, but in practice these objects are approximated by vectors of numbers, that is, simulations. Consequently, when programming for manipulating simulations of unknown quantities, we must write code to manipulate arrays of numeric constants.

And here’s the challenge:

Arrays of simulations are cumbersome objects to work with. Functions that work with vectors will not in general work with matrices, so special versions of the functions need to be written to accommodate matrices of simulations as arguments. For example, a scalar-valued random variable becomes a vector of simulations, and a random vector becomes a matrix of simulations. . . . The usual approach now is to use loops and vector-matrix computations to manipulate the matrix of posterior simulations. This approach is general but awkward and far from transparent . . .

As part of this effort, Jouni wrote a package in R called “rv” that allows direct manipulation of random variable objects, so for example if alpha and beta are random variables, each a vector of length J, and gamma is a scalar, then you can compute something like delta = (alpha + gamma)/beta, which would be a vector of length J. Internally, alpha, beta, and delta would be stored as S x K arrays, and gamma would be stored as an array of length S, where S is the number of simulation draws in your simulation frame, but the user would be able to work with these without having to keep track of that indexing; it would all be implicit.

The rv package was just a demonstration, but then more recently Paul Bürkner, Jonah Gabry, Matthew Kay, and Aki Vehtari implemented the idea in the “posterior” package in R. They call it rvar: The Random Variable Datatype. I still haven’t integrated rvar() into my Bayesian workflow, but I’d like to do it.

The point I want to make here is that this random variable datatype moves us closer to what I’d call “truly probabilistic programming.”

Here’s what I mean by that.

Stan (and Bugs, PyMC, etc.) are called “probabilistic programming languages,” characterized by Wikipedia as “based on the declarative specification of probabilistic models, for which inference is performed automatically.” Fair enough. In Stan, the variables defined in the parameters block, transformed parameters block, model block, and generated parameters block represent random variables with a joint density that is specified by the Stan program.

But, although this is called a probabilistic program as the term is used, it doesn’t fully feel like probabilistic programming in the sense of allowing me to grab random variables and work with them directly. I can do all that within the Stan program in the generated quantities block, but at a certain remove. Stan and other probabilistic programming languages are set up to run in batch mode: you specify a model and then the inference engine runs and spits out posterior simulations (or some other distributional summary).

To me, a truly probabilistic programming environment would allow me to interactively work with random variables, and to work with them as random variables, not as arrays of simulation draws. I want to be able to multiply two matrices without having to loop through the simulations. Part of this is about clean code, but more generally it’s about workflow. Any mathematical or graphical operations we can do on numbers, we should be able to do with random variables, with the simulations all hiding under the hood. Kind of by analogy to how we can do floating-point math on the computer without (usually) having to worry about how all these digits are being stored.

This will require tools such as rvar() and also some experience with the workflow. The payoff is . . . you could call it “better postprocessing of Bayesian simulations,” but I think that would understate the potential benefits of being able to work directl with uncertainty without that extra array dimension getting in your face at every step.

“You’re kind of a hero at my company”

Posted on November 10, 2026 8:20 PM by Phil

This post is by Phil Price, not Andrew.

The work I’m proudest of is not necessarily the work that is most useful or best known or that others consider the best. A few years ago I learned that a model I threw together in about 15 minutes, and whose simplicity I had been quietly embarrassed by for years, was used by a successful company. But I just learned that the model is much more widely used and appreciated than I thought. And I’ve just lost my two steady consulting clients, who were together responsible for about 60% of my income each year, so now is the time to do some self-promotion. And I think there are broader lessons here about the virtues of simplicity, and how we sometimes aren’t good judges of our own work, and maybe some other themes besides. So here’s the story.

About three years ago, a friend and colleague named Steve Greenberg died of cancer, only 64 years old. Steve was a very hands-on kind of guy (as you can learn from his obituary, ) with a deep practical understanding of HVAC systems, which made him very valuable in our quasi-academic setting where a lot of us understood the theories but wouldn’t know the first thing about how to actually install or manage the equipment. Steve was also a committed bike commuter. At his “celebration of life” I wore a bike jacket and a hard hat to recognize those aspects of Steve. Of course my hard hat has my name on it, not Steve’s name.

In the drinks and snacks portion of the party, after the speeches and slideshows, a stranger came up to me and said “Phil Price! You’re kind of a hero at my company.” This was, to say the least, surprising. OK, I’ve done some consulting work for various companies, but none of them would consider me a “hero,” not even remotely close, so what could this possibly be about? But it turns out they weren’t talking about my consulting work at all, but about some work I had done at LBNL years ago: inventing the “Time-of-week and Temperature” (TOWT) model to predict building electric load. “That model is the basis for most of what my company does to measure and verify energy-saving measures in buildings.” That’s not an exact quote, I’m sure, but it’s how I remember it.

Here is an explanation of the model that is clearer and more compact than anything I ever wrote about it. Or, if you want to see more detail and some extensions of the model and how to use it for various purposes, here’s a video of talk I gave back in 2015 (11 years ago!) as I was leaving LBNL to become a consultant. (For me, the video doesn’t display right in Safari but works fine in Firefox, not sure what’s up with that).

The TOWT model is useful in spite of its simplicity, or perhaps even because of its simplicity. It’s very easy to explain, it’s pretty robust, it’s very easy to code up and to fit, it is easy to extend in various ways (e.g. by including additional explanatory variables). But I never bring it up, it has barely been mentioned on my CV (although I’m fixing that), and I haven’t previously bragged about / emphasized my role as creator of the TOWT model, because I’m kinda sheepish about it. As you can see from the web page above or the video, the model is extremely simple: there’s a time-of-week indicator variable, and a piecewise-continuous linear function of temperature…and that’s it, at least in the simplest version of the model. I literally came up with it in an afternoon; actually less than an afternoon, it was more like 15 minutes from conception to implementation. It is very simple and I’ve worried that if I use it as an example of what I’ve done, people will think that level of sophistication is all I can do!

But maybe that’s a silly thing to worry about. People had been analyzing the same kinds of data for decades and in all that time nobody had come up with this model (as far as I know!), simple though it is. And jeez, I had no idea it was so widely used, it’s not just that one company. This is what a search turns up, for tools and methods that rely on the TOWT model:

RM&V (Reference Measurement & Verification Tool): An open-source M&V reference tool released by LBNL that implements piecewise linear regression using TOWT and machine learning methods for multi-site facility analysis.
NMECR: An R-based package by kW Engineering implementing peer-reviewed Normalized Metered Energy Consumption (NMEC) algorithms, featuring an expanded TOWT model that supports operational schedules.
OpenEEmeter / OpenDSM: An open-source analytics toolkit (formerly OpenEEmeter, now operating via OpenDSM) that incorporates the CalTRACK hourly methodology, which is directly adapted from LBNL’s TOWT model for demand response and energy efficiency baselines.
Energy Twin (ET Analytics / ET Core): A commercial energy analytics platform and Python/R library package (ecpEnergyTwinAnalytics) that uses custom TOWT-based implementations to derive operational KPIs and anomaly detection across commercial portfolios.
Power TakeOff (Verify™): An automated M&V 2.0 platform that uses TOWT-based modeling to isolate weather versus schedule loads for utility energy efficiency programs.
CalTRACK Hourly Methodology: A standardized protocol for calculating normalized metered energy savings using an adapted hourly TOWT architecture for commercial and residential buildings.
NMEC (Normalized Metered Energy Consumption): A protocol used by utilities (such as in California) where the TOWT model serves as a benchmark standard to quantify weather-normalized energy savings.
TOWST (Time-Of-Week, Solar and Temperature): An extension of the standard TOWT model that integrates physical constraints like solar gain and building orientation via iterative clustering to improve prediction accuracy.

So, OK, I’m in. That’s my model. I’ll try not to be sheepish about it anymore.

As far as the general lessons to be learned, I think there are a few:

Simplicity can be good. So what if the TOWT model was originally conceived and implemented in a few minutes. Indeed, that’s a good thing! If a simple model does the job, use the simple model.
Don’t hide your light under a bushel basket. If you’ve done something that people find useful, it doesn’t matter if you don’t think it’s your best work or if it doesn’t demonstrate the full range of your skills. This shamelessly self-promotional post was prompted by comments from people who looked at my LinkedIn profile and the letters I’ve sent to potential clients and told me I wasn’t nearly clear enough about what I’ve accomplished. In trying not to come across as bragging, I was selling myself short.
You often don’t know which of your work matters. I had no idea TOWT was in widespread use until a stranger told me at a memorial gathering. I would have guessed that some of my more elaborate work was my most influential. It wasn’t… as far as I know (but how would I know?). It’s worth writing up even the simple stuff (and publishing it, or anyway putting out there somehow).
“It only took 15 minutes” leaves out the work before it. The model took 15 minutes to conceive and implement because I’d spent a couple of months looking at building energy data and knew what mattered and what didn’t, and of course I was proficient in linear regression and other stats methods. Maybe I under-valued my work because I forgot about this fact.

This post is by Phil

Survey Statistics: Sampling from a Bayesian Menu

Posted on November 10, 2026 4:00 PM by shira

A week ago Andrew wrote about Alan Zaslavsky, who recently passed. Alan’s obituary has information about an upcoming celebration of life.

Alan Mark Zaslavsky Profile Photo

Andrew highlighted Alan’s work on the census. In Sampling from a Bayesian Menu, Alan wrote about objections not only to model-based adjustments for undercount, but the Supreme Court’s ruling in 1999 against even purely “design-based” estimates using sampling for nonresponse followup.

Alan was on the panel who wrote “Sampling for Nonresponse Follow-up: Achieving Adequate Precision at Acceptable Cost”. They explained that after initial response to the census, the bureau conducts nonresponse follow-up on remaining housing units. In the 1990 census, this was massively delayed, costing both money and accuracy. They recommended that the 2000 census should conduct nonresponse follow-up on a random sample of nonresponding housing units for unbiased estimation of the total number of people in all nonresponding housing units. This is what Alan described in Sampling from a Bayesian Menu as the “design-based” approach that was rejected by the Supreme Court in 1999.

Despite this bad (old) news, the paper is fun to read, and will make you hungry.

Constructing a latent-variable model to fit measurement data: A Stan story, part 2

Posted on November 10, 2026 11:11 AM by Andrew

Following up on our post from the other day, we were trying to come up with parameters for the following measurement-error model:

– People have latent a distribution of attractiveness:
latent_attractiveness ~ normal(0, sigma_latent)
– Attractiveness for each person is measured twice, with different mean levels and different levels of added variance:
At time t, observed_attractiveness ~ normal(mu_obs_t + latent_attractiveness, sigma_obs_t), for t = 1, 2
– Attractiveness at each time is then recorded as a discrete variable, rounded to the nearest integer and constrained to fall between 1 and 5: That is, if observed attractiveness is less than 1.5, it is recorded as 1, if observed attractiveness is between 1.5 and 2.5, it is recorded as 2, etc.
– The only information easily accessible is recorded attractiveness binned into three categories: 1 or 2, 3 or 4, and 5. This summary is available for each person at both time points.

Here are the data, from the first two waves of the Add Health survey:

So, here’s the deal. We want to choose five parameters (sigma_latent, mu_obs_1, mu_obs_2, sigma_obs_1, sigma_obs_2) to fit 9 data points. Actually 8 data points because what we’re fitting are the proportions that sum to 1.

This shouldn’t be so hard, right?

Actually, it is kind of hard! At least, it’s hard to get a fit that gets real close to the data table.

1. Trial and error to find parameter values that match the data

I started with trial and error, first computing the mean and sd of the attractiveness ratings for each wave of the data (approximating the three binned levels as 1.5, 3.5, and 5 on the 1-to-5 scale), then setting mu_obs_1 and mu_obs_2 to the mean values and partitioning the variance so that sd_latent and sd_obs are roughly equal. In R:

discretize <- function(a) {
  ifelse(a < 2.5, 1, ifelse(a < 4.5, 2, 3))
}
N_table <- rbind(c(208, 640, 40), c(501, 9879, 464), c(51, 1439, 1432))
N_data <- sum(N_table)
long_data <- array(NA, c(N_data, 2))
long_data[,1] <- rep(c(1,2,3, 1,2,3, 1,2,3), as.vector(N_table))
long_data[,2] <- rep(c(1,1,1, 2,2,2, 3,3,3), as.vector(N_table))
long_data_1to5 <- ifelse(long_data==1, 1.5, ifelse(long_data==2, 3.5, ifelse(long_data==3, 5, NA)))
means <- apply(long_data_1to5, 2, mean)
sds <- apply(long_data_1to5, 2, sd)
mu_obs <- means
sigma_latent <- sqrt(0.5*mean(sds^2))
sigma_obs <- sqrt(sds^2 - sigma_latent^2)

The result:

> print(round(mu_obs, 2))
[1] 3.68 3.59
> print(round(sigma_latent, 2))
[1] 0.54
> print(round(sigma_obs, 2))
[1] 0.61 0.45

I then estimate the probabilities in the 9 cells by simulating a million pairs of observations from this model, first simulating latent attractiveness, then the two observed attractiveness measures, then rounding:

N <- 1e6
y_latent <- rnorm(N, 0, sigma_latent)
y_obs <- array(NA, c(N, 2))
for (j in 1:2) {
  y_obs[,j] <- rnorm(N, mu_obs[j] + y_latent, sigma_obs[j])
}
y_discrete <- discretize(y_obs)

And here's the result:

> print(table(y_discrete[,1], y_discrete[,2]))
         1      2      3
  1  18715  54792    393
  2  41167 674479  55053
  3    653 110776  43972

OK, it's hard to compare this to the data table because the numbers are so big, so let's rescale so that it sums to the number of people in the data:

> print(table(y_discrete[,1], y_discrete[,2])*N_data/N)
              1           2           3
  1  274.249610  802.921968    5.759022
  2  603.261218 9883.815266  806.746662
  3    9.569062 1623.311504  644.365688

This is hard to read so I'll round it:

> print(round(table(y_discrete[,1], y_discrete[,2])*N_data/N))
       1    2    3
  1  274  803    6
  2  603 9884  807
  3   10 1623  644

And, for convenience, let's look again at the data:

> print(N_table)
     [,1] [,2] [,3]
[1,]  208  640   40
[2,]  501 9879  464
[3,]   51 1439 1432

Comparing the cells [1,3] and [3,1] in the simulated table to the observed data, we see that our simulations don't have enough variation between the two measurements. So let's rejigger by decreasing sigma_latent and increasing sigma_obs:

sigma_latent <- sqrt(0.25*mean(sds^2))
sigma_obs <- sqrt(sds^2 - sigma_latent^2)

Simulate again with these new parameter values:

y_latent <- rnorm(N, 0, sigma_latent)
y_obs <- array(NA, c(N, 2))
for (j in 1:2) {
  y_obs[,j] <- rnorm(N, mu_obs[j] + y_latent, sigma_obs[j])
}
y_discrete <- discretize(y_obs)

And here's what we get:

> print(round(table(y_discrete[,1], y_discrete[,2])*N_data/N))
   
       1    2    3
  1  145  896   37
  2  687 9592 1011
  3   54 1827  404

Now the [1,3] and [3,1] entries are basically OK, but [1,1] and [3,3] are too small.

I tried some more tinkering with the five parameters but couldn't get the table to line up. And then I thought: Hey, why not cut out the middleman and just solve the optimization problem directly?

2. Stan program: first try

It goes like this:

data: the observed table, N_table
parameters: the five parameters, sigma_latent, mu_obs, sigma_obs
model: compute table(y_discrete) and compare to N_table

In this model, table(y_discrete)*N_data/N represents E(N_table|parameters), and we were computing this using simulation . . . but you can't simulate random data in the Stan model block! Or we could compute the expectations directly, but these are two-dimensional normal integrals, and there's no formula for that.

So I did a hack. You can't do simulations in the Stan model block, but you can simulate random numbers offline and pass them into the Stan model. Here's our Stan program, table_match.stan:

functions {
  int discretize (real y) {
    if (y < 2.5)
      return 1;
    else if (y > 2.5 && y < 4.5)
      return 2;
    else
      return 3;
  }
}
data {
  array[3, 3] int N_table;
  int N;
}
transformed data {
  int N_data = 0;
  for (j1 in 1:3) {
    for (j2 in 1:3) {
      N_data += N_table[j1,j2];
    }
  }
  vector[N] z_latent;
  matrix[N, 2] z_obs;
  for (n in 1:N){
    z_latent[n] = std_normal_rng();
    for (j in 1:2) {
      z_obs[n, j] = std_normal_rng();
    }
  }
  real a = 1.*N_data/N;
}
parameters {
  vector[2] mu_obs;
  real s_latent;
  vector[2] s_obs;
}
transformed parameters {
  matrix[3,3] N_sim = rep_matrix(0, 3, 3);
  for (n in 1:N){
    real y_latent = s_latent * z_latent[n];
    array[2] int y_discrete;
    for (j in 1:2){
      real y_obs = y_latent + mu_obs[j] + s_obs[j] * z_obs[n,j];
      y_discrete[j] = discretize(y_obs);
    }
    N_sim[y_discrete[1], y_discrete[2]] += 1;
  }
  N_sim = a*N_sim;
}
model {
  mu_obs ~ normal(3,3);
  s_latent ~ exponential(1);
  s_obs ~ exponential(1);
  for (j1 in 1:3) {
    for (j2 in 1:3) {
      N_table[j1, j2] ~ poisson(N_sim[j1, j2]);
    }
  }
}

At the very end we're using the Poisson likelihood to compare the observed data to the model. And we include weak priors on the model parameters to keep everything stable.

We run the program from R:

library("cmdstanr")
table_match <- cmdstan_model("table_match.stan")
table_data <- list(N_table=N_table, N=1000)
fit <- table_match$sample(data=table_data, parallel_chains=4, refresh=0)
print(fit, max_rows=30)

Aaaand . . . it doesn't work:

   variable     mean   median       sd      mad       q5      q95  rhat ess_bulk ess_tail
 lp__       77677.57 77880.76 10570.22 13247.66 63077.85 91926.34  3.66        4        5
 mu_obs[1]      0.24     0.10     0.61     0.51    -0.42     1.22  3.94        4       11
 mu_obs[2]      0.47     0.40     1.64     2.24    -1.47     2.54  3.70        4       11
 s_latent       2.78     1.82     2.52     1.37     0.49     7.00  3.42        4       12
 s_obs[1]       1.84     1.60     1.18     1.08     0.46     3.68  3.47        4       13
 s_obs[2]       1.65     1.70     1.36     1.93     0.18     3.03  3.43        4       11
 N_sim[1,1]  8290.04  8865.67  1305.70   640.92  6096.06  9319.94  5.54        4       NA
 N_sim[2,1]   890.16   593.49   669.06   358.48   351.70  2007.60 32.33        4       NA
 N_sim[3,1]   468.94    80.60   711.48    65.18    14.65  1699.86 13.53        4        4
 N_sim[1,2]  1759.81  1472.73  1368.41  1292.70   219.81  3883.31  5.95        4        4
 N_sim[2,2]   852.42   622.80   655.06   369.34   219.81  1948.98  9.66        4       NA
 N_sim[3,2]   322.47   307.73   251.19   336.75    29.31   644.78  8.22        4        4
 N_sim[1,3]   626.65   256.44   779.88   260.71    14.65  1963.64 18.31        4       NA
 N_sim[2,3]   408.84   300.41   339.12   271.58    58.62   952.51 16.39        4       NA
 N_sim[3,3]  1034.67   381.00  1331.88   304.16    29.31  3326.46 16.98        4       NA

What about Stan in optimize setting?

fit_optimize <- table_match$optimize(data=table_data)

Nope:

Initial log joint probability = 69172.8 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
       1       69172.8             0       9.10807    0.001749       0.001       45    
Chain 1 Optimization terminated with error: 
Chain 1   Line search failed to achieve a sufficient decrease, no more progress can be made
Warning: Fitting finished unexpectedly! Use the $output() method for more information.

Damn.

3. Stan program: revision

What went wrong? The above code is mathematically correct, but the discreteness of all these simulations leads to difficulty with the Hamiltonian Monte Carlo. The problem is that, due to the use of fixed error terms and a finite N, continuous changes in the parameters result in jumps in the expected data and thus jumps in the likelihood.

But then I realized something: I can't compute two-dimensional normal integrals in closed form, but one-dimensional normal integrals are trivial. So, instead of passing in the vectors of random numbers z_latent, z_obs[,1], and z_obs[,2], I can just pass in z_latent, and then we can use the normal cdf function in Stan to integrate over z_obs[,1] and z_obs[,2]. This offers two advantages: first, by integrating over two of the three dimensions analytically rather than by simulation, we're computing the cell probabilities more accurately and we don't need such a large N; second, the target function will now be a smoothly continuous function of the parameters, so HMC should work better.

Here's our revised Stan program, table_match_2.stan:

data {
  array[3, 3] int N_table;
  int N;
}
transformed data {
  int N_data = 0;
  for (j1 in 1:3) {
    for (j2 in 1:3) {
      N_data += N_table[j1,j2];
    }
  }
  vector[N] z_latent;
  for (n in 1:N){
    z_latent[n] = std_normal_rng();
  }
  real a = 1.*N_data/N;
}
parameters {
  vector[2] mu_obs;
  real s_latent;
  vector[2] s_obs;
  real sigma;
}
transformed parameters {
  matrix[3,3] N_sim = rep_matrix(0, 3, 3);
  for (n in 1:N){
    real y_latent = s_latent * z_latent[n];
    matrix[3,2] p;
    for (j in 1:2) {
      p[1,j] = normal_cdf(2.5 | y_latent + mu_obs[j], s_obs[j]);
      p[2,j] = normal_cdf(4.5 | y_latent + mu_obs[j], s_obs[j]) - p[1,j];
      p[3,j] = 1 - (p[1,j] + p[2,j]);
    }
    N_sim += p[,1] * p[,2]';
  }
  N_sim = a*N_sim;
}
model {
  mu_obs ~ normal(3,3);
  s_latent ~ exponential(1);
  s_obs ~ exponential(1);
  sigma ~ exponential(1);
  for (j1 in 1:3) {
    for (j2 in 1:3) {
      sqrt(N_table[j1, j2]) ~ normal(sqrt(N_sim[j1, j2]), sigma);
    }
  }
}

You can see from the last line of code that I also changed the likelihood. The Poisson likelihood is appropriate if the model is correct, but it won't be correct; by fitting on the square root scale but including a free parameter, sigma, we're allowing for overdispersion.

We call it from R:

fit_2 <- table_match_2$sample(data=table_data, parallel_chains=4, refresh=0)
print(fit_2, max_rows=30)

Aaahhhh, much better:

   variable    mean  median     sd    mad      q5      q95 rhat ess_bulk ess_tail
 lp__        -25.68  -25.29   2.22   2.08  -29.93   -22.80 1.00     1009     1687
 mu_obs[1]     3.79    3.79   0.08   0.07    3.67     3.92 1.00     2090     2047
 mu_obs[2]     3.66    3.66   0.07   0.06    3.56     3.78 1.00     1762     2092
 s_latent      0.68    0.68   0.06   0.05    0.59     0.77 1.00     2562     1803
 s_obs[1]      0.51    0.51   0.08   0.07    0.39     0.64 1.00     2562     1894
 s_obs[2]      0.27    0.28   0.11   0.10    0.08     0.44 1.00     1791     1139
 sigma         4.22    4.07   1.01   0.93    2.84     6.14 1.00     1511     2172
 N_sim[1,1]  338.69  326.65 114.39 103.68  173.55   539.20 1.00     2433     2009
 N_sim[2,1]  434.37  418.15 164.48 152.14  199.93   722.66 1.00     2775     2271
 N_sim[3,1]    1.63    0.52   3.59   0.70    0.01     6.71 1.00     2207     2454
 N_sim[1,2]  588.17  570.14 193.09 180.59  304.30   926.40 1.00     2580     2098
 N_sim[2,2] 9734.95 9745.35 520.62 488.97 8858.11 10563.91 1.00     3395     2533
 N_sim[3,2] 1637.06 1621.71 322.01 300.94 1134.19  2184.29 1.00     3159     2436
 N_sim[1,3]    0.66    0.13   1.92   0.18    0.00     2.98 1.00     2309     2220
 N_sim[2,3]  595.11  575.31 193.51 175.52  319.79   954.80 1.00     2496     2634
 N_sim[3,3] 1323.36 1298.78 300.80 278.37  863.93  1847.16 1.00     1967     1755

Although, if you look at the expected values in the table (N_sim), it's not so great! The [1,3] and [3,1] values are 0.66 and 1.63, much lower than the data values of 40 and 51.

I think the problem is with the Poisson or square-root likelihood, that it allows relatively large error for small-count cells.

4. Stan program: another revision

So I'll try a new model, this time with a lognormal likelihood. This will go into the file table_match_3.stan. I'll just show you the one line in the model block because that's the only part of the program that's changed:

      log(N_table[j1, j2]) ~ normal(log(N_sim[j1, j2]), sigma);

We call it from R:

fit_3 <- table_match_3$sample(data=table_data, parallel_chains=4, refresh=0)
print(fit_3, max_rows=30)

Which yields:

   variable    mean  median      sd     mad      q5     q95 rhat ess_bulk ess_tail
 lp__         -6.98   -6.42    3.35    3.34  -13.32   -2.51 1.01      339      559
 mu_obs[1]     4.12    4.00    0.59    0.37    3.48    5.25 1.00      819      721
 mu_obs[2]     3.86    3.74    0.52    0.22    3.38    4.89 1.01      761      600
 s_latent      0.82    0.69    0.48    0.22    0.40    1.69 1.01      434      306
 s_obs[1]      0.98    0.95    0.37    0.28    0.46    1.58 1.00     1040      744
 s_obs[2]      0.55    0.51    0.34    0.23    0.09    1.16 1.01      777      833
 sigma         0.84    0.74    0.38    0.30    0.42    1.56 1.01      448      810
 N_sim[1,1]  488.28  286.72  716.36  190.52   91.35 1536.07 1.00     1071      888
 N_sim[2,1]  666.97  552.96  480.97  294.04  220.03 1474.32 1.00     1772     1300
 N_sim[3,1]  115.33   64.10  308.25   42.04   19.43  302.69 1.01     1794     1046
 N_sim[1,2] 1043.18  944.89  607.16  461.48  313.54 2133.92 1.00     1109     1153
 N_sim[2,2] 5820.30 6119.15 2055.18 2028.03 1832.99 8729.56 1.01      432      610
 N_sim[3,2] 2887.52 2773.01 1308.13 1259.89  941.96 5226.20 1.00     1142     1137
 N_sim[1,3]   69.16   36.91  145.99   24.02   12.50  203.99 1.00     1730     1098
 N_sim[2,3] 1139.58  969.93  709.41  479.43  390.12 2416.68 1.00     1478     1611
 N_sim[3,3] 2423.68 1726.31 2045.37 1221.81  495.16 7118.74 1.01      526      635

Or, in optimize setting:

fit_optimize <- table_match$optimize(data=table_data)
print(fit_3_optimize, max_rows=30)

   variable estimate
 lp__           1.62
 mu_obs[1]      3.86
 mu_obs[2]      3.68
 s_latent       0.54
 s_obs[1]       0.82
 s_obs[2]       0.52
 sigma          0.38
 N_sim[1,1]   209.05
 N_sim[2,1]   549.51
 N_sim[3,1]    49.04
 N_sim[1,2]   931.32
 N_sim[2,2]  8175.10
 N_sim[3,2]  2736.90
 N_sim[1,3]    30.44
 N_sim[2,3]   984.07
 N_sim[3,3]   988.56

Not bad! Not perfect---you can see that the [3,3] cell is much too big. But it'll do.

5. Conclusion

The model doesn't quite fit! Maybe it could make sense to try longer-tailed distributions such as the t with low degrees of freedom. But for our original purpose, which is to understand the impact of measurement error on a nonlinear relationship, this will be just fine. Also we've learned something about statistical computing.

Can peer review survive unrestricted AI use by reviewers?

Posted on November 9, 2026 2:07 PM by Jessica Hullman

This is Jessica. Reviews from ICLR (one of the three biggest ML conferences) came out last Friday, and a large number of them appear to have been primarily AI-generated, at least if we believe what authors are reporting on X. I believe it. On my own submissions, about half had Pangram scores over 50% (with some papers getting majority reviews with 90+% Pangram scores). There are other signs as well, like appendix mention reportedly being nearly twice as common in reviews compared to last year.

Was this rampant AI use in review surprising? The official review policy for ICLR 2027 permitted “limited and responsible use of AI tools to assist human writing of peer reviews and meta-reviews,” but said it was “mandatory to disclose use of AI tools,” with the review form soliciting this information and asking that if you use AI you also report your original self-written paper assessment and any LM interactions.

It’s a little mind-boggling to me that anyone would still think asking people to disclose AI use is a viable policy. By now it should be common knowledge that when you say something like “mandatory to disclose use of AI tools” that will be interpreted as “suggested to disclose,” and no one will do it.

It’s also a bit mind boggling that no one anticipated this mess. At ICML 2026, many reviewers who deliberately selected into a pool that would not be allowed to use AI were caught using AI by watermarking. Turns out that when you pump high numbers of papers at each reviewer their good intentions start to crumble. And ICLR received a record number of abstract submissions – over 60,000, suggesting total valid submissions probably surpassed last year’s total near 20,000.

When reviewers feel comfortable telling authors “I have used AI tool to do the review since I do not have enough time”* something is broken. Good for that reviewer for being honest, but who is sending the impression that an AI review you produced but can’t understand is a viable way to fulfill your reviewer agreement?

This situation makes a few things clear. One is that it’s irresponsible to give reviewers high review loads but not restrict how reviewers can use AI, at least in a field like ML where many reviewers are inexperienced and required to do multiple reciprocal reviews.

The second thing is that it’s irresponsible not to prioritize preserving review over hurting authors’ feelings. ICLR did introduce a cap of 20 submissions per author, plus a limit of one paper per author if none of the authors qualified for reviewing. But it’s too conservative. Compare to the policy adopted by STOC, a CS theory conference, which this year capped submissions at five per author, dropped anonymity, required authors to submit a video in which they explain their paper, and reserved the right to desk reject for poor exposition. We can debate whether these are the right steps, but at least someone on that organizing committee had some balls.

Why letting external reviewers use AI freely is a bad idea

At this point, every author is (or should be) using AI to check their own work if they care about quality. Frontier models can find errors that you and your co-workers miss. It can suggest limitations that reviewers might bring up. I don’t send out papers anymore if I haven’t run them through Refine.ink.

So what is the point of reviewers, with less expertise on what the authors did than they have (and in ML, we’re often talking junior PhD students), meat-proxying to send them more AI-generated feedback? Seemingly none.

Consider what we lose if we ban AI in peer review. The main risk would seem to be to peer review’s supposed function as a soundness filter. LLMs enhance a reviewer’s ability to comb through 50+ pages of paper and appendix to verify that what is presented is accurate. But even if the authors are lazy and miss some problems themselves, letting individual reviewers feed papers to LLMs using basic prompts to “write a review for [conference X]” is far from the best solution. Strong LLMs too can miss major errors when not directed what to look for (see e.g., examples from Nihar Shah’s work).

Error identification and basic soundness checks are much better addressed by the venue itself, using validated AI tools. Centralize checks like whether key results can be reproduced from the provided data/code, verify whether equations, proofs, or calculated statistics are correct, detect leakage or violated assumptions, make sure key claims can be traced to evidence, and so on. There’s precedent for all of this in existing research, see Table 1 here or arguments others have made for using AI in review for such purposes here and here.

Some of these automated checks can be used to determine if a paper should be desk rejected. Some will be more controversial (and they should be, as I discuss in more detail here, because there are no hard-and-fast filters to separate good from bad science). Anything contestable can be surfaced for reviewers.

Naturally this can shift what cues get attention during review. When we give reviewers a shared summary representation of the paper, we correlate errors relative to what we’d get if they relied on their private signals alone. But these are problems that machine learning researchers can define and study! See e.g., our work formalizing what it means to have a good LLM-implemented checklist (i.e., use review corpora to learn a set of checks that reduce the combination of bias and variance under LLM-assisted review relative to the variance you’d have from human-only review). We can also educate reviewers on what their role is–-figure out what matters that can’t be verified directly across the paper (significance, novelty, etc.) and decide how to weight the surfaced checks.

But can’t good reviewers write better reviews with AI?

One counterargument to banning reviewers from using AI at their discretion is that we are crippling review quality. It’s easier to get more information from a paper if you can query an LLM about it while you read. A reviewer may reach better insights if they can rely on LLMs to help them run simulations or produce counterexamples to understand the method or theory. It’s not that all AI use is bad. It depends on whether the critiques or suggestions that AI surfaces are within the convex hull of the reviewers’ expertise. If they are, the value of their review might be amplified by AI use. But when they’re not, we get reviewers confidently reporting claims from an LLM that they cannot themselves verify, which the authors may already have decided were not good advice in their own interactions with AI. This is where review becomes pointless.

To retain the value of AI for engaging with papers, venues may also want to provide access to LLMs for interacting with the paper. This makes it possible to preserve some provenance in where claims in reviews are coming from, which authors and ACs can use to appraise reviews as needed. It becomes possible for venues to run checks that reviewers are contributing some kind of unique information to the process, rather than parroting directly. If reviewers are adding minimal information, then those reviews should be removed.

Ultimately, we need reviewers to be using their own judgment to weigh in on the significance and importance and novelty of ideas. Otherwise fields may drift in directions that are misaligned with what the humans working in those fields value. The problem with unchecked AI use is that even when people might want to stay in control of their ultimate judgments, they are not always good at recognizing when they are taking on biases reflected by LLM output. There was a study out of Anthropic recently where they endowed LLMs with preferences that made them lean toward particular recommendations, then had people use them to do a task. The majority of the time, people adopted the preference that was endowed to their LLM, and many of them also attributed the recommendations to themselves.

Why are we still protecting authors’ rights to submit slop?

Reviewers turning to AI en masse is a symptom of conferences being hesitant about desk reject policies. For example, at the top ML conferences, one can still submit a mostly AI-produced paper and expect a team of human reviewers to provide feedback. I’ve heard many students who are reviewing for ICLR this year complain about how many papers in their review piles were hard-to-read AI slop.

Organizers add some minimal friction for authors, and expect the reviewers to suck it up and deal with the remaining onslaught. But when you don’t prioritize reviewers’ ability to actually do their job by requiring more authorial accountability, you hasten the “review death spiral”, where as review gets noisier, authors are further incentivized to submit low quality papers, in the hopes that they get lucky.

There are better ways to deal with this problem than the ultra-conservative don’t-hurt-any-feelings approaches we’ve been seeing. We should not be seeing teams of undergrads submitting 9 papers at a time only to withdraw them once they see the initial reviews, as reported by one AC. For example, start desk rejecting papers that are flagged by Pangram as mostly AI-written, or that are very hard to read. If an author can’t take the time to signal that there is a human on the other end who is accountable for the work, do we really owe them 3+ other humans’ time providing feedback?

I really hope that policies tighten. As imperfect as peer review has always been, removing people from it entirely is not going to lead us to better science.

*From a screenshot of a review shared on X before the author took it down.

Elizabeth Holmes, Jensen Huang, etc.: If all this had happened fifteen to thirty years ago, I’m guessing that all these people would’ve been Democrats. What changed?

Posted on November 9, 2026 10:24 AM by Andrew

We recently saw You Can See Everything—yes, I highly recommend it, as with The Odyssey, it’s as good as everybody says, it lives up to the hype—and then there was the recent news about tech executive Jensen Huang. Put these two stories together and you get the above charming bit of incoherence. “Bring back the gallows,” indeed. Billy Evans seems to be discriminating in his thoughts on the recommended severity of punishment for white-collar crime.

In the movie, Billy Evans, the husband of disgraced medtech fraudster Elizabeth Holmes, says that he’s a fan of the late conservative firebrand Charlie Kirk, and, as you can see below, the couple are all-in on Donald Trump. Huang too is closely linked to the president, who in turn is a promoter of Wild West ventures such as crypto and betting markets.

While watching You Can See Everything and reading the news, I had this thought: If all this had happened fifteen to thirty years ago, I’m guessing that all these people would’ve been Democrats. Big tech money is right-coded today but was left-coded in the (Bill) Clinton and Obama eras.

The question is, why has this happened? The trend is particularly notable when juxtaposed with the general shift of college-educated voters away from the Republican party during the past two decades.

One difference, perhaps, is that information technology is more closely tied to the military than it was twenty years ago. Back then, the defense-contractors and tech were connected, but I think less intertwined. So it could make sense for Apple, Google, Microsoft, etc., to side against George W. Bush and the Iraq War, while current AI companies don’t want to go against a militarily-active government. Taking over Greenland is more exciting to today’s tech investors; twenty years ago they were more interested in opening up consumer markets around the world, and an unpopular foreign invasions were not a useful step in that direction.

Another take on it is that business investors are inclined to go with the winner. In 2008, this was Obama; in 2024, it was Trump. And I do think this is part of the story but not the whole thing. In 2008, Obama’s team had the high-tech luster and the Republicans were the grumpy outsiders. Since 2016, it’s been the opposite. Even setting aside Jensen Huang and Elizabeth Holmes, not to mention Elon Musk and that space mirrors guy, there’s a general sense that conservatives are now the tech cheerleaders and liberals are the new Debbie Downers.

There’s a story to explain any given individual case, or, more generally, you can say that U.S.-style conservatives support something closer to an unregulated free market, while U.S.-style liberals are more concerned about the balance of power between business and labor, so in that sense this all makes sense. One theory that’s been floating around is that tech executives have moved to the right in reaction to efforts by the left-wingers in their technical workforce to reduce executives’ freedom of action. Go back twenty or thirty years, and the Democratic party was led by figures such as Bill Clinton, Al Gore, and Barack Obama, all of whom were friendly not just to big business but to the tech industry in particular. This brings us closer to political scientist Thomas Ferguson’s view of national politics as reflecting divisions within business executives in competing industries. Also, tech lords were rich a generation ago—remember all those jokes about Bill Gates being rich—but not nearly as rich in aggregate as they are now, so that has to be part of the story too.

It’s also interesting to compare how different people react to right-coded business fraud vs. left-coded academic fraud.

“Causal Effects Versus Causal Mechanisms: Two Traditions With Different Requirements and Contributions Towards Causal Understanding”

Posted on November 8, 2026 9:55 AM by Andrew

Following up on our post from last year on the importance of mechanistic models, Jim Grace points to his above-titled paper in Ecology Letters. From the abstract:

When methods are described as ‘causal,’ emphasis is increasingly placed on statistical techniques for isolating associations so as to quantify causal effects. In contrast, natural scientists have historically approached the pursuit of causal knowledge through the investigation of mechanisms that interconnect the components of systems. . . . the case is made that an explicit assessment of existing mechanistic knowledge should be an initial step in causal investigations.

Grace also has this webpage.

I’m posting this because I continue to be interested in the connections between mechanistic and statistical models of causality, and I continue to think that the fields of statistics, economics, and policy analysis took a wrong turn by associating “causal inference” with black-box statistical analysis. I understand the appeal of the idea–indeed, the very first applied statistics class I ever took, 40 years ago, was Rubin’s course on causal inference–but in many ways it’s been a very long detour that only in recent years have we been able to take in a more useful direction.

2026 compared to other dramatic midterm election reversals of the past half-century

Posted on November 7, 2026 9:53 AM by Andrew

As I said a few days ago, you don’t need me anymore to do your fancy post-election data analyses. So instead I can offer speculation (Could the 2026 election results lead to a decrease in political polarization?) and meta-commentary (NYT columnist offers epically stupid take on 2026 midterm elections).

But we can also learn from some really simple analyses.

Comparing 2026 to previous midterm swings

How big of a deal was the 2026 national election?

To look at consequences, you can count the number of seats won by the two parties. To look at the aggregate opinions of voters, it makes sense to look at the national popular vote for the House of Representatives, where the Democrats received ~~56.2% of the two-party vote, a swing of 7.5 percentage points from their House vote share in 2024~~. Update: Checking the latest, the Democrats received 57.1% of the House votes, so a swing of 8.4% from their 2024 share.

Is ~~56.2~~ 57.1% with a swing of ~~7.5~~ 8.4pp a little or a lot? One way to judge this is to compare to the midterm stunners of recent decades:

2010 is the election that put the brakes on the Obama government: the Tea Party movement convinced the voters that government expansion was going too far. In that election, the Republicans won 53.5% of the two-party vote, a 9.0 percentage point swing from 2008.

2006 was the year that George W. Bush’s post-9/11 high finally wore off. The Democrats’ share of the two-party vote in that election: 54.1%, a 5.5pp swing from the previous election.

1994 was the year that Newt Gingrich and the Republicans offered a clear alternative to the liberal policies of the Clinton administration, taking control of Congress with 53.5% of the two-party vote, a 6.1pp swing in their favor. 53.5% is lower than the other numbers I’ll be showing today, but it arguably understates the Republicans’ accomplishment in 1994, in that this was their first stint in the majority in forty years.

1982 isn’t remembered much now, because the Republicans roared back two years later, but, with the worst economic crisis since the Great Depression, and President Reagan having very low job approval, the Democrats did very well at the ballot box, winning 56.0% of the two-party House vote, a 4.1pp swing from the previous election.

1974 was the post-Watergate election, bringing many fresh faces into Congress. The Democrats capitalized on economic stagnation and discontent with the criminal doings of the Nixon administration to gather 56.9% of the two-party vote, a 4.2pp swing.

In summary, the Democrats’ performance in 2026—a 7.5 percentage point swing yielding 56.2% of the two-party vote—was indeed exceptional, at or above the level of previous repudiations of the incumbent party during the midterm elections.

This can be understood in various ways.

One way that 2026 was different is that this time the incumbent party was really the incumbent party. Unlike in those previous campaigns, this time the Republicans had unified control of the federal government—the legislative, executive, and judicial branches—and the incumbent party was itself highly unified (as can be seen, for example, in the Senate confirmation of some extreme cabinet officers). So, voters who wanted party balancing had a clear choice in a way that wasn’t so obvious in previous years where there was divided government.

The other way is to compare to historical baselines. The Democrats’ victories in 1974 and 1982 cemented their status as the default majority party, but that’s not as impressive as the Republicans’ feat in 1994, and then in 2010, of moving from the minority to gain firm control of Congress. The Democrats achieved this in 2026.

What comes next?

In most of the above cases, the midterm winner did not follow up with a decisive win at the next presidential election. 2006, followed by Obama’s victory two years later, was the only exception. Consider the other cases:
After 1974, the Republican party was left for dead: it had been the majority party for awhile and then was wiped out in the congressional elections. But, two years later, the not-particularly-popular Gerald Ford came within 2 percentage points of being elected president.
After the repudiation of Republicans in 1982, the economy came back to life—as the president quipped, now that things are going better, “they don’t call it Reaganomics anymore”—and the presidential election was a landslide.
And the Republicans’ congressional victories in 1994 and 2010 were followed by Clinton and Obama doing just fine in their subsequent reelection campaigns.

That said, the Republicans’ successes in 1994 and 2010 represent one possible template: the triumphant return of a party that had been given up for dead. In both cases, the victorious congressional party did not follow up with a presidential election win, but the party remained strong enough to determine much of the political agenda going forward.

The take-home message based on this history: the Democrats should value their congressional win for itself and not assume they’re on a glide path for 2028. Conversely, the Republicans should figure out how best to work with the Democratic majority without thinking of themselves as being about to lose all power.

Nationalization of politics

Along with partisan polarization has come a nationalization of the vote, which has led to two things. First, the vote swings in presidential elections have been much smaller. After their losses in 1982, the Republicans rebounded a bit in 1984, upping their House vote share to 47.4%, a gain of 3.4 percentage points. In that same year, Reagan was reelected with a two-party vote share of 59.2%!

That sort of thing doesn’t happen anymore, in two ways. First, with partisan polarization, you wouldn’t expect either party’s presidential candidate to come close to 60% of the vote. Second, congressional elections are nationalized, and you wouldn’t expect to see a stunning 12 percentage point gap between a party’s congressional and presidential vote share.

Constructing a latent-variable model to fit measurement data: A Stan story, part 1

Posted on November 6, 2026 9:39 AM by Andrew

This simulation is a little different than what came before. For my earlier post I simulated hypothetical future data from the beauty-and-sex-ratio study. For my new post I’m simulating data from a hypothetical strong relationship between two variables. Here I’m entirely focused on the question of what happens to a nonlinear (but monotonic) pattern when independent measurement error is added to the predictor. To make the patterns clear, I will assume a large effect; indeed, I’m only concerned here with the effect size, not its variation.

The one thing I’ll preserve from the previous simulation is the distribution of attractiveness measurements in the sample: 2%, 5%, 45%, 37%, and 11% in the five attractiveness categories.

I want to set up a model where there’s an assumed true nonlinear relationship between some outcome of interest and a continuous underlying latent attractiveness variable.

The nonlinear model isn’t hard to come up with. The hard thing is creating a distribution for the latent attractiveness, under the constraint that it be roughly consistent with the data.

I’ll assume a symmetric model, a normal distribution on the measurement scale. This might seem like a strong assumption but it seems reasonable enough given the data, which have an average that’s between 3 and 4. The tails are asymmetric but that’s just because the scale has 5 as an upper bound.

My first step is to compute the mean and sd of the observed data:

n <- round(2792*c(.02,.05,.45,.37,.11))
mean_ratings <- sum(n*(1:5))/sum(n)
sd_ratings <- sqrt(sum(n*((1:5) - mean_ratings)^2)/sum(n))

These come to 3.5 and 0.83. Assuming random zero-mean measurement error, the latent attractiveness distribution should have the same mean but a lower variance. We should subtract 1/12 to account for the rounding (as this is approximately equivalent to adding an independent uniform(-0.5,0.5) random variable, which has a variance of 1/12) and then more to account for measurement error.

But how much to subtract? For this, we need to make some assumption about variation in the attractiveness measurement. Attractiveness seems pretty subjective to me (also it can vary a lot over time, as we know from those high-school yearbook photos of people who later became movie stars). I gotta assume something here, so let's assume that the standard deviation of the measurements is 0.5. I'm picking this because 0.5 is a large value, but it's not as large as the sd of 0.83 in the data. Indeed, 0.5^2 + 1/12 = 0.33 and 0.83^2 = 0.69, so I'm assuming that about half of the variance in the data is coming from measurement error.

But is 0.5 a reasonable value? To answer that question, it helps to have data.

The data at hand come from wave 3 of the Add Health study. I already downloaded those data awhile ago in preparing my earlier analysis, but I've had difficulty working with the numbers. I can't remember why, but there's no single dataset and no single codebook. It's kind of a mess. I spent a few minutes going through to see if there were any repeated measurements of attractiveness there, but I wasn't sure where to look.

So I googled *reliability of add health attractiveness measure*. I was hoping to find the names of the relevant variables in the survey, and that didn't come up in the search. What did appear was a bunch of research papers on the topic. I looked at a few of them and then came across this one: "Ugly criminals," by Naci Mocan and Erdal Tekin. I'm not so interested in the topic of the paper---I'll leave it to others to decide if its conclusions are reasonably supported by its data---; what's relevant here is that they discuss the Add Health attractiveness measure. Which is indeed how it came up in my search.

They looked at data from respondents whose attractiveness had been recorded in two waves. That's exactly what we need for a reliability analysis:

There's only one mild annoyance here: they pooled categories 1 and 2, and they pooled categories 3 and 4. So we don't have the full 5 x 5 table, only a 3 x 3 collapsing of it.

But we can deal with that. We'll just simulate and see what happens.

First we write a couple of simple functions:

round_and_bound <- function (a, lower=1, upper=5){
  pmin(upper, pmax(lower, round(a)))
}

clumped <- function(a) {
  ifelse(a==1 | a==2, 1, ifelse(a==3 | a==4, 2, ifelse(a==5, 3, NA)))
}

We set up the parameters for the latent data:

mu_obs <- c(3.87, 3.69)
s_latent <- 0.53
s_obs <- c(0.82, 0.52)

Then we simulate two independent measurements on each of 14654 people---that's the number of cases in the above table:

y_latent <- rnorm(N, 0, s_latent)
y_obs <- array(NA, c(N, 2))
for (j in 1:2) {
  y_obs[,j] <- round_and_bound(rnorm(N, mu_obs[j] + y_latent, s_obs[j]))
}

And then we look at the two margins and the table showing the joint distribution of the two measurements:

print(table(clumped(observed_1)))
print(table(clumped(observed_2)))
print(table(clumped(observed_1), clumped(observed_2)))

I also looked at the unclumped data but I won't bother showing you that here.

Here are the results:

> print(table(clumped(observed_1)))
    1     2     3 
  668 11186  2800 

> print(table(clumped(observed_2)))
    1     2     3 
  723 11221  2710 

> print(table(clumped(observed_1), clumped(observed_2)))
       1    2    3
  1  227  437    4
  2  490 9292 1404
  3    6 1492 1302

These don't quite match the table above. The simulations have not enough "unattractive" people and too many "very attractive." Also they two waves in the simulations are essentially the same (by construction, as we simulated them from the same distribution), but the real data show a few more unattractive people and a lot more very attractive people in wave 1.

There would be a way to set up my simulation to alter its parameters and get the closest possible match to the observed two-way table---I guess it wouldn't be hard to program this in Stan, actually! But here I took the lazy way out and just tinkered with the numbers to get things to work. It took about five tries of playing with the inputs, and eventually I converged on this:

sd_measurement <- 0.5
mu_latent <- 3.8
sd_latent <- 0.67
latent_beauty <- rnorm(n_sim, mu_latent, sd_latent)
observed_1 <- round_and_bound(rnorm(n_sim, latent_beauty, sd_measurement))
observed_2 <- round_and_bound(rnorm(n_sim, latent_beauty - 0.1, sd_measurement*0.6))

As you can see, I hacked the simulations to give the wave 2 measurement a lower mean and lower variance.

Here's the result:

> print(table(clumped(observed_1)))
    1     2     3 
  845 10839  2970 

> print(table(clumped(observed_2)))
    1     2     3 
  731 11852  2071

> print(table(clumped(observed_1), clumped(observed_2)))
       1    2    3
  1  330  514    1
  2  400 9727  712
  3    1 1611 1358

OK, now I have the margins pretty good, but the joint table doesn't have enough people whose measurements shifted between the lowest and highest categories.

I played around some more with the parameters of the model but I couldn't get the counts to line up with the data.

So I bit the bullet and did what I should've done from the start, which was to write a Stan program to find a best fit. That was its own struggle---as is often the case, I learn a lot from these little examples---so I'll devote a separate post to it.

Here's what I ended up using:

mu_obs <- c(3.87, 3.69)
s_latent <- 0.53
s_obs <- c(0.82, 0.52) 
y_latent <- rnorm(N, 0, s_latent)
y_obs <- array(NA, c(N, 2))
for (j in 1:2) {
  y_obs[,j] <- round_and_bound(rnorm(N, mu_obs[j] + y_latent, s_obs[j]))
}

And here's the result, which is still not perfect (as will be discussed in my followup post) but is enough for my goals here:

> for (j in 1:2) print(table(clumped(y_obs[,j])))
   1     2    3 
 1225 9684 3745  

   1     2    3 
 815 11857 1982 

> print(table(clumped(y_obs[,1]), clumped(y_obs[,2])))
       1    2    3
  1  208  990   27
  2  563 8141  980
  3   44 2726  975

Which is far from a perfect match to the actual table of responses from the two waves:

> print(N_table)
     [,1] [,2] [,3]
[1,]  208  640   40
[2,]  501 9879  464
[3,]   51 1439 1432

But it's close enough, as we're just trying to make a general point about attenuation of effects under measurement error.

P.S. For the next story, see Constructing a latent-variable model to fit measurement data: A Stan story, part 2.

Could the 2026 election results lead to a decrease in political polarization?

Posted on November 5, 2026 9:41 AM by Andrew

I have some hope, for the following reasons:

1. The incoming Democrats in Congress include some centrists and some on the left. Mix them in with the ones who are already there, and we can expect divisions on general ideological grounds and on specific issues such as Israel/Palestine, data centers, etc.

2. There have always been congressional Republicans who detest Trump, but they’ve shown close to zero defections from the president’s policies. But now that they’re the minority, the forces keeping them in line will be weaker. Also they are getting the point that general elections as well as primaries are threats to them.

My hope is that there will be some important issues that split within each party, leading to cross-party coalitions. I’m not saying this would be a cure for polarization—back in 1993, the vote on NAFTA split the two parties, and that did nothing to stop the Newt Gingrich-led polarization that took over soon after—but maybe it would be a step forward.

To put it another way, I think of there being two components to partisan polarization. One component is a natural consequence of connections between issues: there’s a logic to the left and right representing different coalitions and being represented by two political parties. The other component is not so good, it’s the lockstep conformity that leads to candidates taking extreme positions to appeal to primary election voters, presidents crafting policies to satisfy moderates just enough to get the threshold of required votes, and party-line votes to confirm unqualified cabinet officials. Maybe some open disputes on high-profile issues will make pure partisan voting less of an expected thing, so that legislators can vote against their party leadership without it being such a big deal, which in turn could motivate moderation among their leaders. Which does not mean that non-centrist policies on the left or right could not be passed, just that they’d require coalition building of some sort..

Political polarization has been steadily increasing for the past thirty-five years, so it seems like it would be a mistake to bet against more of the same. But trends can occasionally reverse; nothing lasts forever. Back when my colleagues and I were starting to do research in political science, the incumbency advantage had been steadily increasing for forty years . . . but then the trend reversed.

Also, partisan polarization is a funny thing: some people are strongly motivated to vote for Republicans because they hate the Democrats and what they stand for; others are strongly motivated to vote for Democrats because they hate the Republicans etc. There aren’t too many voters who actually like either party, at least not on a national level. To borrow a physical analogy, you could say that political polarization is maintained by tension, not compression.

P.S. Somewhat relatedly, NYT columnist offers epically stupid take on 2026 midterm elections.

P.P.S. I think I blew it by not using a clickbait title such as “Polarization is so over.”

“So in that sense the silly statistics in the paper is no big deal, it serves as a sort of rhetorical flourish?”

Posted on November 4, 2026 9:36 AM by Andrew

Richard Feldman writes:

I’m a long-time reader and sometimes commenter on the blog. I’ve learned a lot – and been entertained – by many of the posts.

A paper was just published in Science on a theme that is potentially of interest to your blog readership, and that has the potential to enter the news, given its political nature. It’s not my research area, but the statistical practice to which I object crops up in my discipline – ecology – quite often. In fact, I’m reviewing a paper right now that does something similar.

The paper is “Partisan disparities in the funding of science in the United States.”

What has me question my sanity as someone trying to engage in good statistical practice is that the authors ran a regression with the covariates of interest (party control of the presidency, house, and senate) and many “control” variables and then a simple model with the covariates of interest and just one “control” variable. (A whole other issue is using the term control to describe the additional covariates and to refer to party control…in the same figure!). In the main figure in the text, they present both results – from the simple and full model – but clearly emphasize the results from the simple model. They justify their emphasis on the simpler model with what sounds ludicrous to me: “Collinear controls may absorb relevant variance, and the descriptive meaning of a regression coefficient becomes increasingly opaque as more variables are included in the model.”

I call the above the “have your cake and eat it too” approach to interpreting regressions. You propose that there are a series of factors to consider when measuring the effect of a variable of interest, but then interpret the effect as if those other factors weren’t important. I think there tends to be a default sentiment that coefficients in multiple regression are difficult to interpret – in the paper I’m reviewing, this is what the authors say – because of a lack of understanding that it is the effect when the other variables are at a pre-specified value (usually zero or their mean, but now with R packages like marginaleffects it is easy to specify those values). Am I right to have an unease about the approach of running multiple and simple regressions at the same time, but interpreting just the simple regression?

In addition to the above, the paper analyses an untold number of models – “alternative specifications” – and uses an NHST approach without considering the inflation of Type I error.

I agree that this is a pretty empty paper, published in a prestigious journal just because it’s on a topic of interest to the editors and also it’s wrapped in a clean-looking package. In addition to the flaws noted above, there’s also the problem that they don’t start their analysis earlier, a period during which the Democrats were in continuous control of Congress for several decades. On the other hand, their general point seems correct, that until this year science funding has had steady bipartisan support. So in that sense the silly statistics in the paper is no big deal, it serves as a sort of rhetorical flourish?