I think that any non-linear fitting scheme would do, as. e.g. the Newton-Raphson scheme. You can use even numeriacl derivatives (finite difference) instead of analytical ones. You can estimate the static and high frequency limiting permiitivities (or moduli) and the approximate relaxation time form the data. You may start with the Debye approximation (exponential decay) and run the Newton-Raphson procedure. It can be programmed relatively easily even for complex numbers. You will get a good estimate in less than 10 cycles. The improvement of the parameter estimation can be judged by comparing the residuals in subsequent cycles.

http://www.mathworks.com/help/optim/nonlinear-least-squares-curve-fitting.html

http://mathfaculty.fullerton.edu/mathews/n2003/Web/Newton%27sMethodMod/Newton%27sMethodMod.html

Again, let us assume that you are sure you are measuring bulk effect and not electrode-sample interface as well.

The quickest way I found to fit the non-Debye response is to display data as complex capacitance and fit Havrilliak general phenomenological function (see A. Jonscher book on dielectric relaxation in solids - KWW function is just one of these universal responses in real time). NO involved numerical fitting procedures are needed, manual fit (5 to 10 loops) is usually sufficient to capture the essential parts of the model when compared with experiment.

Best regards

Petr