Some SPICEs have optimization utilities (SpiceOPUS, for one)
but any SPICE can be used "in the feedback position" to get
at a set of fitting params less intelligently (that is, using some
of your own). Nested .modif loops, and pick the best fit by
eye, like. This may take more than one pass to "find the range"
but it's doable.
You might attack it in layers, like getting the Vth ought to be
either a "pick off from model files" or a simple fit to the simple
transistor measurement testbench for VT.
But of course you want to fit against the proper equation, the
independent variable seems to be missing entirely? Or is vth
supposed to be the independent variable, the form perhaps
=(a0+a1*vth+a2*vth^2) - normal power series? But that
would seem to be leaving a lot of leakage influences off the
table... or, you choose this formula to represent only one
specific PVTXYZ condition and embed all of that to the
coefficients.