GPT4 is even better at Stan than GPT3 and even better at C++. It’s sooo much better at longer chains of reasoning, writing style, and being able to maintain coherence. Check this two turn dialogue to get it to write a non-centered nonlinear mixed effects model. I followed it up with a one-prompt C++ question to get it to write a full N-body simulation, which it nailed. I’m sure it could’ve used other integration method like Euler’s method or Runge-Kutta.
It’s easy at this point to have it start riffing on the model by adding effects, removing effects, switching to varying slope/varying intercept, etc. It can store up to about 5000 words of context when generating the next 1000 words (and there’s a version you can pay even more for with a roughly 25K word context).
The fact that it can keep all the indexes coherent and line everything up is still mind-boggling to me. GPT4 is qualitatively better at this than GPT3 was. This thing was never taught to code Stan per se. It was just taught which kinds of randomly generated code look like real code! Nevertheless, it picked up things that look like the weakly informative priors people use in practice for logistic regression.
Of course, it has trained on everything available about Stan online, which isn’t much. It’s way better at C++, Python, and R. I included a C++ example using physics below that you should appreciate.
I had GPT3 coach me through writing a speech input, speech output interface to GPT. I built it in Python using the Qt GUI toolkit, the Elevenlabs text to speech API, and the Whisper speech to text API. I now have a (push to talk) speech interface to GPT.
BC: I want to write a Stan program to define a nonlinear mixed effects model. I want to use a logistic regression to predict response to an intervention in hospital trauma centers. Each binomial measurement is a number of successes out of number of trials. Each hospital has a state in which it is located (50 states), and a trauma center level (1, 2, or 3), both of which get random effects with independent hierarchical priors in the model. Then there is are fixed effect for the average income of the hospital patients and the average severity of the cases it sees. And there should be an overall intercept. The fixed effects get fixed priors.
GPT4: To write a Stan program that defines a nonlinear mixed-effects model with the described specifications, you can follow the below structure:
data {
int<lower=0> N; // Number of observations
int<lower=0> J; // Number of hospitals
int<lower=0> K; // Number of states
int<lower=1, upper=3> level[J]; // Trauma center level
int<lower=1, upper=K> state[J]; // State in which each hospital is located
int<lower=0> successes[N]; // Number of successes
int<lower=0> trials[N]; // Number of trials
real<lower=0> avg_income[N]; // Average income of hospital patients
real<lower=0> avg_severity[N]; // Average severity of the cases
}
parameters {
real<lower=0> intercept; // Overall intercept
real<lower=0> beta_income; // Fixed effect of income
real<lower=0> beta_severity; // Fixed effect of severity
real<lower=0> sigma_state; // State random effect standard deviation
real<lower=0> sigma_level; // Trauma center level random effect
standard deviation
vector[K] state_effect_raw; // State random effects (non-centered)
vector[3] level_effect_raw; // Trauma center level random effects
(non-centered)
}
transformed parameters {
vector[K] state_effect = sigma_state * state_effect_raw; // State random effects
vector[3] level_effect = sigma_level * level_effect_raw; // Trauma center level random effects
}
model {
// Priors for fixed effects
intercept ~ normal(0, 10);
beta_income ~ normal(0, 10);
beta_severity ~ normal(0, 10);
// Priors for random effects standard deviations
sigma_state ~ cauchy(0, 2.5);
sigma_level ~ cauchy(0, 2.5);
// Priors for non-centered random effects
state_effect_raw ~ normal(0, 1);
level_effect_raw ~ normal(0, 1);
// Model
for (i in 1:N) {
int hospital_idx = ...; // Index of the hospital corresponding to
the i-th observation
successes[i] ~ binomial_logit(trials[i], intercept + beta_income *
avg_income[i] + beta_severity * avg_severity[i] +
state_effect[state[hospital_idx]] +
level_effect[level[hospital_idx]]);
}
}
This Stan program defines a logistic regression model with fixed effects for average income and average severity of cases, and random effects for the state and trauma center level. The program takes in the number of observations, number of hospitals, number of states, trauma center level, and state of each hospital, as well as the number of successes, number of trials, average income, and average severity of the cases. The model parameters include the overall intercept, fixed effects for income and severity, standard deviations for the random effects, and non-centered random effects for state and trauma center level.
The transformed parameters block calculates the actual random effects for the state and trauma center level by scaling the non-centered random effects with their respective standard deviations. In the model block, priors are set for the fixed effects and the standard deviations of the random effects. The non-centered random effects are given standard normal priors.
Finally, the model block iterates through each observation, and for each observation, it identifies the corresponding hospital index. The binomial_logit distribution is used to model the number of successes out of the number of trials, incorporating the fixed effects, the state random effect, and the trauma center level random effect.
BC: That’s really good, but there are still two errors. (1) There should not be lower bounds on the random effects because we are going to put the linear predictor through an inverse logit transform (implicit in the binomial_logit distribution you used). (2) The indexing is wrong in the model block. The number of observations and the number of hospitals is the same and the index can just represent the hospital (e.g., in order to index the random effects). Please fix these errors and give me just the corrected Stan program.
GPT4:
My apologies for the errors in the previous response. I have corrected the errors in the Stan program below:
data {
int<lower=0> N; // Number of observations (equals to the number of hospitals)
int<lower=0> K; // Number of states
int<lower=1, upper=3> level[N]; // Trauma center level
int<lower=1, upper=K> state[N]; // State in which each hospital is located
int<lower=0> successes[N]; // Number of successes
int<lower=0> trials[N]; // Number of trials
real<lower=0> avg_income[N]; // Average income of hospital patients
real<lower=0> avg_severity[N]; // Average severity of the cases
}
parameters {
real intercept; // Overall intercept
real beta_income; // Fixed effect of income
real beta_severity; // Fixed effect of severity
real<lower=0> sigma_state; // State random effect standard deviation
real<lower=0> sigma_level; // Trauma center level random effect standard deviation
vector[K] state_effect_raw; // State random effects (non-centered)
vector[3] level_effect_raw; // Trauma center level random effects
(non-centered)
}
transformed parameters {
vector[K] state_effect = sigma_state * state_effect_raw; // State random effects
vector[3] level_effect = sigma_level * level_effect_raw; // Trauma center level random effects
}
model {
// Priors for fixed effects
intercept ~ normal(0, 10);
beta_income ~ normal(0, 10);
beta_severity ~ normal(0, 10);
// Priors for random effects standard deviations
sigma_state ~ cauchy(0, 2.5);
sigma_level ~ cauchy(0, 2.5);
// Priors for non-centered random effects
state_effect_raw ~ normal(0, 1);
level_effect_raw ~ normal(0, 1);
// Model
for (i in 1:N) {
successes[i] ~ binomial_logit(trials[i], intercept + beta_income * avg_income[i] + beta_severity * avg_severity[i] + state_effect[state[i]] + level_effect[level[i]]);
}
}
I have removed the lower bounds from the random effects in the parameters block and corrected the indexing in the model block, as the number of observations and the number of hospitals are the same.
==================physics time============================
BC: Can you write a C++ program to carry out an N-body problem simulation using the leapfrog integrator for the Hamiltonian dynamics derived from Newtonian mechanics. It should randomly initialize positions and masses and velocities.
GPT4: Here’s a simple C++ program for an N-body simulation using the leapfrog integrator. This program initializes positions, masses, and velocities randomly and then simulates the system using the leapfrog method.
cpp
#include <iostream>
#include <vector>
#include <cmath>
#include <random>
const double G = 6.67430e-11; // Gravitational constant
struct Particle {
double m;
double x, y, z;
double vx, vy, vz;
};
void compute_forces(const std::vector<Particle>& particles,
std::vector<Particle>& forces) {
for (size_t i = 0; i < particles.size(); ++i) {
double fx = 0, fy = 0, fz = 0;
for (size_t j = 0; j < particles.size(); ++j) {
if (i == j) continue;
double dx = particles[j].x - particles[i].x;
double dy = particles[j].y - particles[i].y;
double dz = particles[j].z - particles[i].z;
double d = std::sqrt(dx*dx + dy*dy + dz*dz);
double f = G * particles[i].m * particles[j].m / (d*d);
fx += f * dx / d;
fy += f * dy / d;
fz += f * dz / d;
}
forces[i].x = fx;
forces[i].y = fy;
forces[i].z = fz;
}
}
void leapfrog_step(std::vector<Particle>& particles, double dt) {
std::vector<Particle> forces(particles.size());
compute_forces(particles, forces);
for (size_t i = 0; i < particles.size(); ++i) {
particles[i].vx += forces[i].x / particles[i].m * dt / 2;
particles[i].vy += forces[i].y / particles[i].m * dt / 2;
particles[i].vz += forces[i].z / particles[i].m * dt / 2;
particles[i].x += particles[i].vx * dt;
particles[i].y += particles[i].vy * dt;
particles[i].z += particles[i].vz * dt;
}
compute_forces(particles, forces);
for (size_t i = 0; i < particles.size(); ++i) {
particles[i].vx += forces[i].x / particles[i].m * dt / 2;
particles[i].vy += forces[i].y / particles[i].m * dt / 2;
particles[i].vz += forces[i].z / particles[i].m * dt / 2;
}
}
int main() {
const int n_particles = 100;
const double t_total = 1000.0;
const double dt = 0.01;
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<> mass_dist(1e20, 1e22);
std::uniform_real_distribution<> pos_dist(-1e11, 1e11);
std::uniform_real_distribution<> vel_dist(-1e4, 1e4);
std::vector<Particle> particles(n_particles);
for (int i = 0; i < n_particles; ++i) {
particles[i].m = mass_dist(gen);
particles[i].x = pos_dist(gen);
particles[i].y = pos_dist(gen);
particles[i].z = pos_dist(gen);
particles[i].vx = vel_dist(gen);
particles[i].vy = vel_dist(gen);
particles[i].vz = vel_dist(gen);
}
int n_steps = static_cast<int>(t_total / dt);
for (int step = 0; step < n_steps; ++step) {
leapfrog_step(particles, dt);
// Print the state of the system at each step (optional)
/*
for (const auto& p : particles) {
std::cout << "Mass: " << p.m << ", ";
std::cout << "Position: (" << p.x << ", " << p.y << ", " << p.z << "), ";
std::cout << "Velocity: (" << p.vx << ", " << p.vy << ", " << p.vz << ")\n";
}
std::cout << "-----\n";
*/
}
// Print final state of the system
for (const auto& p : particles) {
std::cout << "Mass: " << p.m << ", ";
std::cout << "Position: (" << p.x << ", " << p.y << ", " << p.z << "), ";
std::cout << "Velocity: (" << p.vx << ", " << p.vy << ", " << p.vz << ")\n";
}
return 0;
}
This program initializes the positions, masses, and velocities randomly, and then simulates the system using the leapfrog method for a given total time `t_total` and time step `dt`. The state of the system can be printed at each step by uncommenting the optional print statements in the loop. The final state of the system is printed after the simulation is complete.
BC: What is the best integration method to use for an N-body problem of the scale you just coded?
GPT4: For N-body problems, especially when simulating gravitational interactions, the best integration methods are those that are symplectic, time-reversible, and have good long-term energy conservation properties.
One of the popular choices for such problems is the leapfrog integrator, which was used in the code I provided earlier. Leapfrog is a second-order symplectic integrator with good energy conservation properties and is relatively simple to implement. It is suitable for many N-body simulations and provides accurate results for a wide range of problems.
However, when dealing with highly eccentric orbits or close encounters between bodies, higher-order symplectic integrators like the Verlet, Forest-Ruth, or Yoshida integrators might be more appropriate. These methods provide better accuracy and stability for such challenging scenarios, albeit at the cost of increased complexity and computational effort.
I asked Bob if he’d like to blog this, and he said, “Overall, everyone’s fatigued from everyone saying how great GPT is, so I’m not so keen to blog it. But hey, feel free if you’d like to.”
So I decided to post it on lag. By the time the post is appeared, perhaps people won’t be talking about GPT anymore. The only thing I added above were the html tags to make the code more readable on the blog.
My overall reaction is the same as Bob’s: I’m impressed. As we see more and more of these chatbot dialogues, it becomes more clear how they work, in general terms. As Bob and others have explained, they’ve pulled tons of text off the internet and are reconstructing it, monkey-on-typewriter style. Oddly enough, the regurgitating nature of the algorithms becomes clearer when the chatbots work well, as above, than when they fail, as in Gary Smith’s examples we discussed awhile ago. So, yeah, it doesn’t understand Stan; it can just find various Stan codes on the web that are connected to various statistical questions, and put these together. Kinda like the way that, if I’m trying to get a widget working in Latex and I’m too lazy to try to figure out what’s going on, I’ll search the web for some relevant code. But . . . damn, yeah, it’s pretty impressive how well this works, how close it comes to a reasonable answer, even while not following the logic that a human would use to answer these questions. Or, maybe it’s not so far from the logic that a human would use to try to get in the right ballpark.
Also, I don’t like the blank lines between paragraphs of the code. When I’m coding, a key resource is real estate on the screen, and I don’t like to waste space in that way. Every blank line represents a line of code somewhere else that I don’t get to see. But we can’t blame the chatbot—it’s just copying all the code-with-blank-lines that’s already on the web!
