Using Claude Code and MATLAB to Find a New Way to Estimate Pi
Bio: Today's guest blogger is Charlie Saunders. Charlie is a Senior Research Scientist in the MathWorks Advanced Research and Technology Office.
Recently, I came across a cool result about calculating π via coin flips. You may have seen other physical processes used to calculate π, too, such as throwing darts at a square and seeing what ratio land within a circle inside it. This one is a fun classroom demo, but the convergence is terrible... and you might not want to hand darts to a classroom of young kids! But this got me wondering whether there other simple random processes that produce πand are suitable to try out in a classroom. And, could I use an AI agent to help me search for one?
I spent an afternoon with Claude Code connected to a live MATLAB session using the MATLAB Agentic Toolkit and MCP server. My goal was to use Claude as a fast experimental partner by having it help brainstorm candidate processes, write and run MATLAB simulations, and then help check the algebra when something looked promising.
The Workflow: Claude Code + MATLAB MCP
Claude Code is a command-line AI agent that can read and write files, run commands, and interact with tools from the terminal. For this project, the important part was that it could write MATLAB code, run it through the MCP connection, inspect the output, and then revise the next attempt without me having to manually intervene too much.
The loop went like this:
- I described the main goal and a few starting approaches to consider: stopping rules, Poisson arrivals, and variants of coin-flip or birthday-problem processes. Claude then suggested some additional families of random processes to try.
- For each candidate, Claude wrote a MATLAB simulation, ran it through the MCP server, and analyzed the results.
- We dropped ideas that did not converge to π, converged too slowly, or appeared to have systematic bias.
- For the methods that survived numerically, I asked Claude to use Symbolic Math Toolbox™ to check whether there was a simple closed-form proof.
This workflow was great, as I could test a half-formed "what if..." idea almost immediately, instead of spending an hour doing algebra or hand-writing experimental code to check an approach only to find it wasn't working.
The Discovery: A Marble Bag Encodes Pi
We started by testing a whole bunch of ideas. I found a birthday-problem approach (how many people until two share a birthday?) that gives π only approximately, with bias that depends on the number of possible birthdays. Counting ties in coin flips was far too noisy to be practical. And an early version of the marble bag where you add one blue marble per round instead of two produced an answer close to π, yet upon further investigation resulted in $ (8 log(2) + 4)/3 $ instead.
The idea we ended with came from exploring hazard-rate processes: random stopping games where the probability of stopping changes each round. After testing several stopping rules in MATLAB, we came across one where the probability of stopping on round k, assuming the game has reached that round, is 1/(2k). At first, this may not sound especially connected to π. But this particular stopping rule has a useful hidden structure. The probability of surviving the first n rounds becomes
$ S(n)=\frac{\binom{2n}{n}}{4^n} $,
which involves the central binomial coefficients. Those same coefficients appear in the Taylor series for $ \arcsin(x) $, and when that series is evaluated at $ x=1 $, $ \arcsin(1)=\pi/2 $. So this was the first real clue that the stopping rule might be more than a numerical accident.
The remaining question was whether we could turn this abstract stopping rule into something simple enough to do physically.
The Marble Bag
Start with 1 red marble and 1 blue marble in the bag. Each round, draw a marble at random. If it is red, the game ends. If it is blue, put it back, add two more blue marbles, and draw again.
Here is what happens:
- Round 1: Bag has 1 red + 1 blue (2 marbles). Draw. P(red) = 1/2.
- Round 2: Bag now has 1 red + 3 blue (4 marbles). Draw. P(red) = 1/4.
- Round 3: Bag now has 1 red + 5 blue (6 marbles). Draw. P(red) = 1/6.
- Round k: Bag has 1 red + (2k-1) blue = 2k marbles. P(red) = 1/(2k).
Record the stopping round τ. Compute the score: $ \tau / (2\tau - 1) $. Repeat many times, and... the average score converges to $ \pi/4 $! (a slice of Pi, if you will).
This sounds reasonable, but does it actually work? I asked Claude to write a quick simulation and run it using the MATLAB MCP server to test this out. Here's the MATLAB code it ran:
N = 10000;
tau_vals = zeros(N, 1);
for t = 1:N
k = 1;
while rand() >= 1/(2*k)
k = k + 1;
end
tau_vals(t) = k;
end
scores = tau_vals ./ (2*tau_vals - 1);
running_est = 4 * cumsum(scores) ./ (1:N)';
fprintf('After 100 games: pi ~ %.4f\n', running_est(100));
The estimate is pretty good! With 100 games, the mean absolute error was only about 0.06, or roughly 2% of π. We also see that the estimator is unbiased in expectation (and doesn't depend on anyone having good aim with a dart!).
Proving It: The Symbolic Math Toolbox
Without leaving the conversation, I then asked Claude to switch from performing numerical analysis with MATLAB to using the Symbolic Math Toolbox to prove the result in closed form. I had not tried this before so I was unsure how well it would work. But to my delight, it wrote the symbolic code, executed it through MCP, and returned an exact answer on the first try.
The One-Line Proof
The stopping probability 1/(2k) at each round produces a PMF:
$ P(\tau = k) = \frac{\binom{2k-2}{k-1}}{k \cdot 2^{2k-1}} $
Claude defined this symbolically and asked MATLAB for the expected score:
syms k
pmf = nchoosek(sym(2*k-2), k-1) / (k * 2^(2*k-1));
score = k / (2*k - 1);
% The one-line proof:
E_score = symsum(score * pmf, k, 1, inf)
% Returns: E_score = pi/4
So with one call to symsum, the Symbolic Math Toolbox returns $ \pi/4 $.
Where Does Pi Come From?
The identity $ \tau/(2\tau-1) = 1/2 + 1/(2(2\tau-1)) $ means:
$ E[\tau/(2\tau-1)] = \frac{1}{2} + \frac{1}{2} \cdot E[1/(2\tau-1)] $.
The Symbolic Toolbox evaluates that inner expectation:
E_inv = symsum(pmf / (2*k - 1), k, 1, inf)
% ans = pi/2 - 1
This traces back to the arcsin series, one of the classical representations of π. Combining:
$ E[\tau/(2\tau-1)] = 1/2 + (1/2)(\pi/2 - 1) = \pi/4 $
If you haven't seen these before: the first series is the Taylor expansion of arcsin(x) at x = 1 (see here for more detail). The second sums the Catalan numbers with geometric weights (see here). The partial fraction step shows how our score function decomposes into one piece from each series, and that's where π sneaks in!
How long will a game last?
We saw from our Monte Carlo simulations that the stopping round τ is heavily right-skewed. You draw red on round 1 about half the time, but occasionally the game runs for 10+ rounds.
We asked Claude to verify that the PMF sums to 1 and explore the moments:
symsum(pmf, k, 1, inf) % ans = 1 (valid PMF)
symsum(k * pmf, k, 1, inf) % ans = Inf (infinite expected game length!)
That second result is a lovely surprise. The game has infinite expected length, even though the score $ \tau/(2\tau-1) $ has a perfectly finite mean.
This does have an implication for the physical version of the game. Most games finish quickly: 50% stop on round 1, and about 75% stop within 5 rounds. But occasionally, a game can keep going for a very long time. That is fine for MATLAB, but tedious if you are drawing marbles by hand. In fact, whilst discussing this post, Mike Croucher had a game that ran for 3,997,135,355 draws: to do that with marbles, you'd need 40,000 metric tonnes of marbles and ~125 years!
I posed this issue to Claude, asking whether there's an unbiased way to cap the game at a set number of rounds. It turns out, there is an elegant solution. If red has not appeared by round M, stop and record a fixed replacement score equal to the average score you would expect from the unfinished tail:
$ c_M = E\left[\frac{\tau}{2\tau-1} \mid \tau > M\right] $.
I felt that 5 rounds seems reasonable. For M=5, this replacement score is approximately $ c_5 = 0.5162 $.
So the practical classroom version works like this: draw at most 5 times. If you hit red, record $ \tau/(2\tau-1) $. If you do not hit red in 5 draws, record 0.5162. About 25% of games are truncated, but the estimator remains unbiased because those unfinished games are replaced by their exact conditional average.
What Made This Work
Looking back, three things came together to make this feasible in an afternoon:
- Claude Code for reasoning and brainstorming. Claude drew on its knowledge and reasoning to propose a non-obvious stopping rule. I probably wouldn't have thought to explore hazard-rate processes as I was not particularly familiar with them before embarking on this project.
- MATLAB via MCP for instant feedback. Every idea was tested in seconds, so bad ideas could be discarded quickly, and good ones refined and verified. Claude could write a simulation, run it, look at the plot, and decide whether to continue autonomously.
- The Symbolic Math Toolbox for proof. Moving from simulation to proof without switching tools turned the empirical hunch into a theorem. The symsum one-liner that returned $ \pi/4 $ was more convincing than any Monte Carlo run.
The project went from my question "could we find a classroom-friendly random process that estimates π?" to a working marble-bag game with simulations, figures, and a closed-form proof in a single afternoon. And I think that's pretty cool!
If you want to try this kind of AI-assisted research yourself, you'll need a coding agent and the MATLAB MCP server. You can find the companion MATLAB scripts for this post here.
댓글
댓글을 남기려면 링크 를 클릭하여 MathWorks 계정에 로그인하거나 계정을 새로 만드십시오.