In an attempt to get a formula, we have to consider some approximations.
Into the allowed cooling air ambient temperature (Ta), the total operating temperature is the addition of Ta plus the machine temperature rise (Tr). See NEMA MG1 – 20.40.
T = Ta + Tr
To calculate Tr based on your statistical data, (at least two sets of readings) we will assume:
Tr1 = K1*Ii^2 + K2
Tr2 = K1*Ij^2 + K2
-------------------------
Were K1*I^2 = Temp rise due to Stator I^2R+RotorI^2R+Stray Load Losses
K2 = Temp rise due to Core loss + Friction and windage loss (constant if voltage and frequency are held constant)
Tr1, Tr2 = recorded temperature rises.
From two sets of winding temperature records, subtract ambient (Ta) and get Tr.
Resolving two simultaneous equations we find K1 and K2.
K1= (Tr1-Tr2)/(Í^2-Ij^2)
K2= Tr1 – I1^2(Tr1-Tr2)/(Í^2-Ij^2)
For a new load current I, calculate Tr using the values calculated for K1 and K2
Tr = k1*I^2 + K2
Adding the air ambient temperature T = Ta + Tr
You could check how close this approach is by comparing against the recorded results of a third set of readings.
I see that this is almost an identical proposition to that by Electicpete.