Find a best fitting curve for some data with no regular pattern

May be something here of use: https://mathematica.stackexchange.com/a/14232/363

data=
{{1, 0.08888888888888872},{2, 0.13842641081508972},{3, 0.17153280014540232},
 {4, 0.19586894934728782},{5, 0.21482500502852542},{6, 0.23017213935351592},
 {7, 0.24294479920054762},{8, 0.25379708095478332},{9, 0.26316771779066082},
 {10, 0.27136463872018062},{11, 0.27861184976714762},{12, 0.28507709754351612},
 {13, 0.29088903838678592},{14, 0.29614834727916052},{15, 0.30093516063430942},
 {16, 0.30531421254734282},{17, 0.30933847047870212},{18, 0.31305176609215322},
 {19, 0.3164907360488469}};

nlm = NonlinearModelFit[data, {a + b Log[c x], c > 0}, {a, b, c}, x];
Print@Normal@nlm;
Print@FindFit[data, {a + b Log[c x], c > 0}, {a, b, c}, x];
Show[Plot[nlm[x], {x, 1, Length@data}, PlotRange -> All],
 ListPlot[data, PlotStyle -> Directive[Red, PointSize[0.02]]]]

See also, http://reference.wolfram.com/mathematica/tutorial/CurveFitting.html - particularly the section: 'Searching for general fits to data'.

In version 10.2 there is a new experimental function which might be what you are looking for: FindFormula.

I suspect that a genetic programming algorithm (symbolic regression) is behind this new feature.

See also my question here: What is behind experimental function: FindFormula?

Functions are like vectors. Actually you can define a vector space over functions. We can describe any vector in terms of 3 independent vector which may not be orthogonal to each other. Similarly, you can fit a function in terms of other linearly independent functions. The easiest case is fitting to a polynomial of order n. Depending on how well your data can be fit you can set n. Another example is fitting to a series of Sine and Cose functions(Discrete Fourier transform). One other example is fitting to Gaussian functions with different mean and standard deviation.

data = {{1, 0.08888888888888872}, {2, 0.13842641081508972}, {3, 
0.17153280014540232}, {4, 0.19586894934728782}, {5, 
0.21482500502852542}, {6, 0.23017213935351592}, {7, 
0.24294479920054762}, {8, 0.25379708095478332}, {9, 
0.26316771779066082}, {10, 0.27136463872018062}, {11, 
0.27861184976714762}, {12, 0.28507709754351612}, {13, 
0.29088903838678592}, {14, 0.29614834727916052}, {15, 
0.30093516063430942}, {16, 0.30531421254734282}, {17, 
0.30933847047870212}, {18, 0.31305176609215322}, {19, 
0.3164907360488469}};

Fit to Cos functions:

model2 = Sum[Subscript[a, n]*Cos[(2*Pi)/T*n*x], {n, 0, 20}];
sol = FindFit[data , {model2, T > 100},  Flatten[{Table[{Subscript[a, i]}, {i, 0, 20}], T}], x]

Fit to polynomial:

model1 = Sum[Subscript[a, i]*x^i, {i, 0, 5}];
sol = FindFit[data , {model1},  Flatten[{Table[{Subscript[a, i]}, {i, 0, 5}]}], x]

Fit to Gaussian functions:

Here I took only 3 Gaussian function, for better fit either you need more Gaussian functions or a better initial guess and constraints.