\(\(\(
\)\(<< 
      Statistics`DescriptiveStatistics`

 (*\ \ \ \ Hermite\ spline\ interpolation\ of\ f_t\ \((x)\)\  = \ \((x - 
                t)\) _ + \(\(^\)\((k - 
                1)\)\)\ \ \ \ \ \ \ *) 
 (*\ \ \ \ on\ [0, 
            1]\ with\ 2  k - 
        4\ interior\ knots\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
 (*\ \ \ \ generated\ as\ order\ statistics\ of\ 2  
          k - 4\ i . i . d . \ uniform \((0, 
            1)\)\ \ \ \ *) 
 (*\ \ \ \ random\ variables\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
 (*\ \ \ \ Restriction : \ \ k\  ≥ \ 
          3\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 


    Print["k :  (2k-1) degree spline interpolation"]

    k = 3

    Print["d = number of digits of accuracy"]

    d = 30

    Print["nn = number of monte-carlo replications"]

    nn = 100


    f[x_, j_] := x^j

    g[x_, y_, 
        k_] := \((\((x - y)\)*UnitStep[x - y])\)^\((2*k - 
            1)\)

    fp[x_, j_] := j*x^\((j - 1)\)

    gp[x_, y_, 
        k_] := \((2*k - 
            1)\)*\((\((x - y)\)*UnitStep[x - y])\)^\((2*k - 
              2)\)

    h[x_, t_, 
        k_] := \((\((x - t)\)*UnitStep[x - t])\)^\((k - 1)\)/\(\((k - 
              1)\)!\)

    hp[x_, t_, 
        k_] := \((\((x - t)\)*UnitStep[x - t])\)^\((k - 2)\)/\(\((k - 
              2)\)!\)


    \(Basis[x_, vec_List, k_] := 
        Join[{1}, Table[f[x, j], {j, 1, 2*k - 1}], 
          Table[g[x, vec[\([j]\)], k], {j, 2, 
              2*k - 3}]];\)


    \(HermiteInterp[x_, coeff_List, vec_List, k_] := 
        coeff . Basis[x, vec, k];\)

    mxlist[coeff_List, vec_List, t_, \ k_] := 
      Table[Abs[
          HermiteInterp[i/1000, coeff, vec, k] - h[i/1000, \ t, \ k]], {i, 1, 
          1000}]



    Error[
        k_] := {
tt = 
          Random[Real, {0, 1}, d]; 
VectauUnif = 
          Sort[Table[
              Random[Real, {0, 1}, d], {j, 1, 
                2*k - 4}]]; 
Vectau = 
          Append[Prepend[VectauUnif, 0], 1]; 
MM = 
          Join[
Table[
              Join[{1}, 
Table[
                  f[Vectau[\([j]\)], i], {i, 1, 2*k - 1}], 
                Table[g[Vectau[\([j]\)], Vectau[\([i]\)], k], {i, 2, 
                    2*k - 3}]], {j, 1, 2*k - 2}], 
Table[
              Join[{0, 1}, 
Table[
                  fp[Vectau[\([j]\)], i], {i, 2, 2*k - 1}], 
                Table[gp[Vectau[\([j]\)], Vectau[\([i]\)], k], {i, 2, 
                    2*k - 3}]], {j, 1, 2*k - 2}]]; 
VecH = 
          Join[Table[h[Vectau[\([j]\)], tt, k], {j, 1, 2*k - 2}], 
            Table[hp[Vectau[\([j]\)], tt, k], {j, 1, 
                2*k - 2}]]; 
HermiteCoeff = 
          LinearSolve[MM, VecH]\ ; 
mtl = 
          mxlist[HermiteCoeff, Vectau, tt, \ k]; 
mx = 
          Max[mtl]
}



    \(Ert = 
        Flatten[Table[
            Error[k], {j, 1, nn}]];\)


    \(zz = Sort[Ert];\)

    Print["mean, median, and harmonic mean of the maximum error (on grid of 1000 pts)"]

    LocationReport[Ert]

    Print["sample variance and standard deviation of the maximum error"]

    DispersionReport[Ert]

    Print["95% quantile of maximum error"]

    qu = Quantile[zz,  .95]

    Print["99% quantile of maximum error"]

    quu\  = Quantile[zz,  .99]

    Print["maximum (across monte-carlo replications) of maximum error"]

    ME = Max[Ert]



    \(hh = Table[k, {k, 1, nn}]*\((1/nn)\);\) // N

    \(J1 = 
        Table[{zz[\([i]\)], hh[\([i]\)]}, {i, 1, nn}];\)

    Print["listplot of empirical distribution of maximal error"]

    ListPlot[J1, PlotRange -> All]

    Delta[x_] := 0 /; x < 0

    Delta[x_] := 1 /; x ≥ 0

    DE[x_] := Table[{Delta[zz[\([i]\)] - x]}, {i, 1, nn}]

    EDF[nn_, x_] := 
      1.  - Apply[Plus, DE[x]]*\((1/nn)\) // N

    Print["empirical distribution of errors, up to 1.1 times 99% quantile"]

    P1 = Plot[EDF[nn, x], {x, 0, quu*1.1}, PlotRange -> All, 
        PlotStyle -> {Thickness[1/100], RGBColor[0.000, 1.000, 0.196], 
            Dashing[{}]}\ ]

    Print["empirical distribution of errors, up to 1.1 times max"]

    P2 = Plot[EDF[nn, x], {x, 0, ME*1.1}, PlotRange -> All, 
        PlotStyle -> {Thickness[1/100], RGBColor[0.000, 1.000, 0.196], 
            Dashing[{}]}\ ]

    \)\)\)