Fit isn't restricted to polynomials. You should choose your basis functions to match the nature of the problem. For example, suppose I believe my data has a jump (of unknown height) at time 0, added to an unknown linear trend with an unknown offset and measurement noise. Something like this:
data = Table[{x,
3 UnitStep[x] - 0.3 x + 2 +
RandomVariate[NormalDistribution[]]}, {x, -10, 10}]
{{-10, 4.5653}, {-9, 4.72512}, {-8, 4.44479}, {-7, 4.18572}, {-6,
2.58606}, {-5, 3.00913}, {-4, 2.09279}, {-3, 3.71923}, {-2,
2.69165}, {-1, 1.55887}, {0, 4.69223}, {1, 3.07118}, {2,
5.43961}, {3, 3.50293}, {4, 2.9189}, {5, 3.70878}, {6, 3.57311}, {7,
3.41627}, {8, 2.79404}, {9, 3.04386}, {10, 2.47882}}
You may ListPlot this if you like: my eyes can't easily see the underlying regularity in all of the noise. Nevertheless, Fit figures it out pretty well:
Fit[data, {1, UnitStep[x], x}, x]
2.10015 - 0.228675 x + 2.55593 UnitStep[x]
It works well here because my data represents a sum of known functions with unknown coefficients, and I fed those known functions to Fit as a basis. More phenomenological uses, like fitting a polynomial or Fourier series to data whose genesis is different, are trickier. High-order polynomials represented as power series (1,x^1,x^2,x^3,...) are especially tricky, as the results are generally numerically unstable.