ChatGPT writes a (quite good) Matlab function.

On Quora a question was asked about individual football match results indicating the actual relative strengths of the teams. I replied with an explanation of the Poisson model of match scores. This model assumes that on the day each team have an expected number of goals (at full time) under the conditions pertaining, and that the actual number of goals scored by each team is a Poisson distributed random variable (RV) with the expected mean.

As an example I ran some Monte-Carlo runs (in Gnu-Octave) of games with the first team having an expected number of goals of 1.2, and the second team having an expected number of 1.5. This gave a probability of team 1 winning of $\approx 0.3$ and of team 2 winning of $\approx 0.45$ and of a draw of $\approx 0.25$. I also ran this again with means of 1.2 and 3 when the results were $\approx 0,12,\ 0.74,\ 0.14$.

In conversation with a fellow dog-walker in Jubilee Park, he explained the success he was having using LLMs as code generators from text prompts. So I thought "how well could one of these LLMs do with this model?"

Here are the prompt I used, and ChatGPT's (free version) response:

----------------------------------------------------------------------------------

Prompt:

would you write a Matlab function for me that computes the probability of a random variable with a Poisson distributes with mean L1 exceeding the value of a random variable with Poisson distribution with mean L1

Response:

 I assume you meant two independent Poisson random variables with means $L_1$ and $L_2$. (Your message says both means are $L_1$, but if they are the same then the probability is simply $(1-P(X=Y))/2)$ by symmetry.)

The probability you want is

$P(X>Y),\qquad X\sim\mathrm{Poisson}(L_1)$  $Y\sim\mathrm{Poisson}(L_2)$.

A numerically stable MATLAB implementation uses the identity

$P(X>Y)=\sum_{y=0}^{\infty}P(Y=y),P(X>y)$,

and truncates the infinite sum once the remaining probability mass is negligible.

function p = poissonGreater(L1, L2, tol)
%POISSONGREATER  Compute P(X > Y) for independent Poisson variables.
%
%   p = poissonGreater(L1,L2)
%   p = poissonGreater(L1,L2,tol)
%
%   Inputs:
%       L1  Mean of X (>=0)
%       L2  Mean of Y (>=0)
%       tol Optional truncation tolerance (default 1e-12)
%
%   Output:
%       p   Probability that X > Y

    if nargin < 3
        tol = 1e-12;
    end

    if L1 < 0 || L2 < 0
        error('Means must be nonnegative.');
    end

    % Upper limit chosen so omitted tail probability is below tol
    ymax = ceil(max(poissinv(1-tol/2,L1), poissinv(1-tol/2,L2)));

    y = (0:ymax)';

    % PMF of Y
    py = poisspdf(y,L2);

    % Tail probabilities P(X > y)
    pxTail = 1 - poisscdf(y,L1);

    % Sum
    p = sum(py .* pxTail);
end

Example

>> poissonGreater(5,3)

ans =

    0.7336

Special case: equal means

If $L_1=L_2=\lambda$, then

$P(X>Y)=\frac{1-P(X=Y)}{2}$,

where

$P(X=Y)=\sum_{k=0}^{\infty}\frac{\lambda^{2k}}{(k!)^2}e^{-2\lambda}$

The function above works for this case as well. For example:

>> poissonGreater(5,5)

ans =

    0.4361

If you intended both means to be (L_1), I can also provide a much faster implementation that exploits the symmetry above and computes (P(X>Y)) directly from the equal-mean formula (or via the Skellam distribution).

--------------------------------------------------------------------------------------

This code as requested by the prompt only gives the probability that the first team wins, but a simple wrapper will give equivalent results to the Monte-Carlo simulation:

function [p1,p2,pE]=FB_2(L1,L2)
%
% my function to fill in the cases where RV2 exceeds RV1
% and that they are equal.
%
  p1=poissonGreater(L1, L2);
  p2=poissonGreater(L2, L1);
  pE=1-p1-p2;
end

Then running this gives:

>> [p1,p2,pE]=FB_2(1.2,1.5)
p1 = 0.3037
p2 = 0.4415
pE = 0.2548
>> [p1,p2,pE]=FB_2(1.2,3.0)
p1 = 0.1212
p2 = 0.7367
pE = 0.1421
>

Which are close enough to Monte-Carlo results for government purposes.

Comments on the ChatGPT code: In some respects it is better than mine in that it does input data validation, also it is better commented (I do not comment experimental code of this sort). Also I just used a script virtually as a calculator:

disp("football_1");  

ngames=10000; 

mu1=1.2; mu2=3;  
r1=rand(1,ngames);  
r2=rand(1,ngames); 

s1=poissinv(r1,mu1);
s2=poissinv(r2,mu2);  

sum(s1>s2)/ngames
sum(s1==s2)/ngames
sum(s1<s2)/ngames  

On the whole I would say I'm impressed by the performance of ChatGPT as a code generator in this instance, it even corrected the typo in the prompt!