Normalise max value of probability function for all frames

I have a working script that plots an animated probability density function (PDF). I have normalised the contour so values fall within -1 and 1. However, I'm calling a separate distribution function for each frame. So the maximum values change and hence the probability will be transformed differently for each frame. I'm hoping to first calculate the PDFs for all frames, compute the max values for all frames and then normalise all frames at once with this 'global' maximum value. Below is the function I use to normalise the PDF for a single frame.

normPDF = (PDFs[0]-PDFs[1])/max(PDFs[0].max(),PDFs[1].max())

The problem can be visualised using the colorbar in the attached animation.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as sts
import matplotlib.animation as animation
from mpl_toolkits.axes_grid1 import make_axes_locatable

DATA_LIMITS = [-30, 30]

def datalimits(*data):
    return DATA_LIMITS  

def rot(theta):
    theta = np.deg2rad(theta)
    return np.array([
        [np.cos(theta), -np.sin(theta)],
        [np.sin(theta), np.cos(theta)]
    ])

def getcov(radius=1, scale=1, theta=0):
    cov = np.array([
        [radius*(scale + 1), 0],
        [0, radius/(scale + 1)]
    ])
    r = rot(theta)
    return r @ cov @ r.T

def mvpdf(x, y, xlim, ylim, radius=1, velocity=0, scale=0, theta=0):
    X,Y = np.meshgrid(np.linspace(*xlim), np.linspace(*ylim))
    XY = np.stack([X, Y], 2)
    x,y = rot(theta) @ (velocity/2, 0) + (x, y)
    cov = getcov(radius=radius, scale=scale, theta=theta)
    PDF = sts.multivariate_normal([x, y], cov).pdf(XY)
    return X, Y, PDF

def mvpdfs(xs, ys, xlim, ylim, radius=None, velocity=None, scale=None, theta=None):
    PDFs = []
    for i,(x,y) in enumerate(zip(xs,ys)):
        kwargs = {
            'radius': radius[i] if radius is not None else 0.5,
            'velocity': velocity[i] if velocity is not None else 0,
            'scale': scale[i] if scale is not None else 0,
            'theta': theta[i] if theta is not None else 0,
            'xlim': xlim,
            'ylim': ylim
        }
        X, Y, PDF = mvpdf(x, y,**kwargs)
        PDFs.append(PDF)

    return X, Y, np.sum(PDFs, axis=0)

fig, ax = plt.subplots(figsize = (10,4))

ax.set_xlim(DATA_LIMITS)
ax.set_ylim(DATA_LIMITS)

line_a, = ax.plot([], [], '.', c='red', alpha = 0.5, markersize=5, animated=True)
line_b, = ax.plot([], [], '.', c='blue', alpha = 0.5, markersize=5, animated=True)
lines=[line_a,line_b] 

cfs = None

def plotmvs(tdf, xlim=None, ylim=None, fig=fig, ax=ax):
    global cfs  
    if cfs:
        for tp in cfs.collections:
            tp.remove()

    df = tdf[1]

    if xlim is None: xlim = datalimits(df['X'])
    if ylim is None: ylim = datalimits(df['Y'])

    PDFs = []

    for (group, gdf), group_line in zip(df.groupby('group'), lines):
        if group in ['A','B']:
            group_line.set_data(*gdf[['X','Y']].values.T)
            kwargs = {
            'radius': gdf['Radius'].values if 'Radius' in gdf else None,
            'velocity': gdf['Velocity'].values if 'Velocity' in gdf else None,
            'scale': gdf['Scaling'].values if 'Scaling' in gdf else None,
            'theta': gdf['Rotation'].values if 'Rotation' in gdf else None,
            'xlim': xlim,
            'ylim': ylim
            }
            X, Y, PDF = mvpdfs(gdf['X'].values, gdf['Y'].values, **kwargs)
            PDFs.append(PDF)

    normPDF = (PDFs[0]-PDFs[1])/max(PDFs[0].max(),PDFs[1].max())

    cfs = ax.contourf(X, Y, normPDF, cmap='viridis', alpha = 0.8, levels=10)

    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.1)
    cbar = fig.colorbar(cfs, ax=ax, cax=cax)
    cbar.set_ticks([-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1])

    return  cfs.collections + [line_a,line_b] 

n = 9
time = range(n)  

d = ({
     'A1_X' :    [13,14,12,13,11,12,13,12,11,10],
     'A1_Y' :    [6,6,7,7,7,8,8,8,9,10],
     'A2_X' :    [7,6,5,7,6,3,4,5,6,6],
     'A2_Y' :    [11,12,11,10,11,12,10,11,10,9],
     'B1_X' :    [8,9,8,7,6,7,5,6,7,6],
     'B1_Y' :    [3,4,3,2,3,4,2,1,2,3],
     'B2_X' :    [13,14,14,14,13,13,13,12,12,12],
     'B2_Y' :    [5,4,3,2,4,5,4,6,3,3],               
    })

tuples = [((t, k.split('_')[0][0], int(k.split('_')[0][1:]), k.split('_')[1]), v[i])
    for k,v in d.items() for i,t in enumerate(time) ]

df = pd.Series(dict(tuples)).unstack(-1)
df.index.names = ['time', 'group', 'id']

interval_ms = 1000
delay_ms = 2000
ani = animation.FuncAnimation(fig, plotmvs, frames=df.groupby('time'), interval=interval_ms, repeat_delay=delay_ms,)

plt.show()

Animation

from Normalise max value of probability function for all frames