Prerequisites

True data generating process p(y=1|x)

g[x_]=0.5 1/(1+Exp[-100 x])+0.25+0.24Sin[20 x];
Plot[{0.5,g[x]},{x,-1,1},AxesLabel->{x,p(y==1|x)}]

To make optimization tractable, compute loss by summing losses over the following points

points=Table[i,{i,-1,1,0.01}];
ListPlot[Thread[{points,g[#]&/@points}]]

Define log - loss and truncated log loss. Truncated log loss is used to avoid numerical problems (getting Log[1/0] during optimization). Using truncated loss has an effect that once the function gets close to 1 or 0, getting it closer no longer decreases the loss.

logloss=-#1Log[((#2)/(1-#2))]-Log[1-#2]&;
loglosst=-#1Log[((Min[#2,0.999])/(1-Min[#2,0.999]))]-Log[1-Min[#2,.999]]&;
Plot[{logloss[0.6,x],loglosst[0.6,x]},{x,0,1}]

Linear/hinge loss. In this case we are considering probability models that never go outside of (0, 1) Since hinge loss and linear loss are identical over that domain, use linear loss

linloss=#2(1-2 #1)&;
Plot[linloss[0.6,x],{x,0,1}]

sqloss=(((#1-#2))^2)&;
Plot[sqloss[0.6,x],{x,0,1}]

Fitter

Construct a function of the form (1/(1+Exp[(Subscript[f, n])(x)])) where n is the n'th order polynomial. Represent (Subscript[f, n])(x) using an orthonormal basis for numerical reasons. (it has an effect of decorellating the parameters, which decreases the condition number of the Hessian). Find best fit by initializing the parameters to correspond to Bayes Optimal classifier and doing gradient descent.

myfito[loss_,dims_,points_,g_]:=Module[{params,polys,f,timing,sol,pt},
params=Thread[(Subscript[p, #])&/@Table[i,{i,0,dims-1}]];
(* orthonormal polynomials*)
polys=LegendreP[#,x]*(2#+1)&/@Table[i,{i,0,dims-1}]; f=1/(1+Exp[-Plus@@Thread[polys*params]]);
realloss=Plus@@(loss[g[#],f/.{x->#}]&/@points) ;
p0=Table[0,{i,0,dims-1}];
p0[[2]]=10; (* set starting point to the optimal classifier *)
timing=Timing[pt=FindMinimum[realloss,Thread[{params,p0}]]];
sol=pt[[2]];
{timing[[1]],f/.sol,params/.sol}
];

Best fit using log - loss and 2 degrees of freedom corresponds to Bayes optimal classifier

f[x_]=myfito[loglosst,2,points,g][[2]];
Plot[{f[x],0.5},{x,-1,1}]

Because some points are on the boundary (g (x) = 1), and our function class does not include those points, there may not be a minimum, however, because we initialize the starting value to function that produces the optimal classifier, we can be sure that the result of minimization gives better loss than the starting point (it might be equal in the case of the Hessian being 0, which is not the case here). You can see that with higher degree of freedom you can decrease log loss, but that increases 0-1 loss.

f[x_]=myfito[loglosst,17,points,g][[2]];
Plot[{f[x],0.5},{x,-1,1}]

Find minimizers values for all degrees of from 2 to 30. There are warnings, likely due to the fact that minimum is outside of the achievable region (step function is not achievable exactly)

fs1=myfito[loglosst,#,points,g][[2]]&/@Table[i,{i,1,30}];
fs2=myfito[linloss,#,points,g][[2]]&/@Table[i,{i,1,30}];

Plot[Evaluate[fs1],{x,-1,1}]

Plot[Evaluate[fs2],{x,-1,1}]

Evaluate approximation error

geterror[g_,f_]:=Module[{cf,n=100000,e},
e=Compile[x,Abs[UnitStep[g[x]-0.5]-UnitStep[f[x]-0.5]]];
Plus@@(e[#]&/@Table[i,{i,-1,1,2/n}]) *2./n
]

geterror[g,ffs1[[29]]]

errors=geterror[g,#]&/@ffs1

ListPlot[errors,Joined->True]