<< Statistics`DescriptiveStatistics`
<< 
      LinearAlgebra`MatrixManipulation`


 (*    Hermite spline and complete spline interpolation of f (x)  =  
        x^(2  k)        *) 
 (*        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                                                      *) 
 (*    with isolation of  the knot configurations corresponding to the           *) 
 (*    10 largest Hermite interpolation errors                                   *) 
 (*    in the form                                                               *) 
 (*    {hermite error,  
        complete error,  {knot_  1 ,  knot_  2 ,   ...  ,  
          knot_ (2  k - 
                4)}}   *) 
 (*        Restriction :   k  ≥  
          3                                          nn       *) 

    
Print["k :  (2k-1) degree spline interpolation"]
k = 7
Print["d = number of digits of accuracy"]
d = 100
Print["nn = number of monte-carlo replications"]
nn = 100


    
B1[x_, j_, l_] := ((j!)/(j - l)!)*
        x^(j - l)
B2[x_, y_, k_, 
        l_] := ((2*k - 1)!/(2*k - 1 - l)!)*((x - y)*
              UnitStep[x - y])^(2*k - 1 - l)
hd[x_, k_, l_] := ((2*k)!/(2*k - l)!)*
        x^(2*k - l)

    
Basis[x_, vec_List, k_] := 
      Join[{1}, Table[B1[x, j, 0], {j, 1, 2*k - 1}], 
        Table[B2[x, vec[[j]], k, 0], {j, 1, 
            2*k - 4}]]

    
Vect[k_, l_, vec_List] :=  
      Join[Table[0, {m, 1, l}], Table[B1[1, i, l], {i, l, 2*k - 1}], 
        Table[B2[1, vec[[i]], k, l], {i, 1, 
            2*k - 4}]]

    
HermiteInterp[x_, coeff_List, vec_List, k_] := 
      coeff . Basis[x, vec, k]
CompInterp[x_, coeff_List, vec_List, k_] := 
      coeff . Basis[x, vec, k]

    
mxlistH[coeff_List, vec_List,  k_] := 
      Table[Abs[
          HermiteInterp[i/1000, coeff, vec, k] - hd[i/1000,  k, 0]], {i, 1, 
          1000}]
mxlistC[coeff_List, vec_List, k_] := 
      Table[Abs[CompInterp[i/1000, coeff, vec, k] - hd[i/1000,  k, 0]], {i, 
          1, 1000}]


Error[
        k_] := {
VectauUnif = 
          Sort[Table[
              Random[Real, {0, 1}, d], {j, 1, 
                2*k - 4}]]; 
Vectau = 
          Append[Prepend[VectauUnif, 0], 1]; 
MMHermite = 
          Join[
Table[
              Join[{1}, 
Table[
                  B1[Vectau[[j]], i, 0], {i, 1, 2*k - 1}], 
                Table[B2[Vectau[[j]], Vectau[[i]], k, 0], {i, 2, 
                    2*k - 3}]], {j, 1, 2*k - 2}], 
Table[
              Join[{0, 1}, 
Table[
                  B1[Vectau[[j]], i, 1], {i, 2, 2*k - 1}], 
                Table[B2[Vectau[[j]], Vectau[[i]], k, 1], {i, 2, 
                    2*k - 3}]], {j, 1, 
                2*k - 2}]]; 
VecHermite = 
          Join[Table[hd[Vectau[[j]], k, 0], {j, 1, 2*k - 2}], 
            Table[hd[Vectau[[j]], k, 1], {j, 1, 
                2*k - 2}]]; 
HermiteCoeff = 
          LinearSolve[MMHermite, VecHermite]; 
Hermitemxl = 
          mxlistH[HermiteCoeff, VectauUnif,  k]; 
Hmx = 
          Max[Hermitemxl]; 
MM1 =  
          AppendRows[DiagonalMatrix[Table[j!, {j, 0, (k - 1)}]], 
            Table[0, {i, 1, k}, {j, 1, 3*k - 4}]]; 
MM2 =  
          Table[Join[{1}, 
              Table[B1[VectauUnif[[j]], i, 0], {i, 1, 2*k - 1}], 
              Table[B2[VectauUnif[[j]], VectauUnif[[i]], k, 0], {i, 
                  1, 2*k - 4}]], {j, 1, 
              2*k - 4}]; 
Vec3 =  
          Join[Table[1, {m, 1, 2*k}],  
            Table[B2[1, VectauUnif[[i]], k, 0], {i, 1, 
                2*k - 4}]]; 
MM3 =  
          Join[{Vec3}, 
            Table[Vect[k, l, VectauUnif], {l, 1, 
                k - 1}]]; 
MMComp =  
          Join[MM1, MM2, MM3]; 
VecC = 
          Join[Table[hd[0, k, l], {l, 0, k - 1}],  
            Table[hd[VectauUnif[[j]], k, 0], {j, 1, 2*k - 4}],  
            Table[hd[1, k, l], {l, 0, 
                k - 1}]]; 
CompCoeff = 
          LinearSolve[MMComp, VecC]; 
Compmxl = 
          mxlistC[CompCoeff, VectauUnif, k]; 
Compmx =  
          Max[Compmxl]; 
mx = {Hmx, Compmx, 
            VectauUnif}

}

Ert = 
        Table[Error[k][[1]], {j, 1, 
            nn}];
Print["Block of hermite error, complete error, knots, 1st MC \
rep"]
Ert1 = Ert[[1]];
Print["Transpose of Ert"]
ErtT  =  Transpose[Ert];
Print["Hermite errors"]
ErtT1 = ErtT[[1]];
Print["sorted Hermite errors"]
Sort[ErtT1];
Print["Ranks of 10 largest Hermite errors"]
BigIO = Ordering[ErtT1, -10]
Print["block of hermite error, complete error, knots of largest 10 \
hermite errors"]
Ert[[BigIO]]

Print["all the hermite errors"]
ErtH = 
        Table[Ert[[i, 1]], {i, 1, nn}];
Print["all the complete spline errors"]
ErtC = 
        Table[Ert[[i, 2]], {i, 1, nn}];
Print["Statistics of Hermite Error and Complete Error"]
Print["Mean Hermite Error, Mean Complete Error"]
{Mean[ErtH], Mean[ErtC]}
Print["Median Hermite Error, Median Complete Error"]
{Median[ErtH], Median[ErtC]}
Print["StandardDev Hermite Error, StandardDev Complete Error"]
{StandardDeviation[ErtH], 
      StandardDeviation[ErtC]}
Print["0.95-q Hermite error, 0.95-q Complete error"]
{Quantile[ErtH, 0.95], Quantile[ErtC, 0.95]}
Print["0.99-q Hermite error,  0.99-q Complete error"]
{Quantile[ErtH, 0.99], Quantile[ErtC, 0.99]}
Print["Max Hermite Error, Max Complete Error"]
{Max[ErtH], 
      Max[ErtC]}



zzH = Sort[ErtH];
zzC = Sort[ErtC];
hh = Table[k, {k, 1, nn}]*(1/nn); // 
      N
JH = 
        Table[{zzH[[i]], hh[[i]]}, {i, 1, 
            nn}];
JC = 
        Table[{zzC[[i]], hh[[i]]}, {i, 1, 
            nn}];
Print["listplot, empirical distribution of hermite errors"]
ListPlot[JH, PlotRange -> All];
Print["listplot, empirical distribution of complete errors"]
ListPlot[JC, PlotRange -> All];
Delta[x_] := 0 /; x < 0
Delta[x_] := 1 /; x ≥ 0
DEH[x_] := 
      Table[{Delta[zzH[[i]] - x]}, {i, 1, nn}]
DEC[x_] := 
      Table[{Delta[zzC[[i]] - x]}, {i, 1, nn}]
EDFH[nn_, x_] := 
      1.  - Apply[Plus, DEH[x]]*(1/nn) // N
EDFC[nn_, x_] := 
      1.  - Apply[Plus, DEC[x]]*(1/nn) // 
        N

Print["empirical distribution function of Hermite Errors"]
P1 = 
      Plot[EDFH[nn, x], {x, 0, Quantile[ErtH, 0.99]}, 
        PlotStyle -> {Thickness[1/100], RGBColor[0.000, 1.000, 0.196], 
            Dashing[{}]}, 
        PlotRange  ->  All ]
Print["empirical distribution function of Complete Errors"]
P2 = 
      Plot[EDFC[nn, x], {x, 0, Quantile[ErtH, 0.99]}, 
        PlotStyle -> {Thickness[1/100], RGBColor[0.50, 0.000, 0.50], 
            Dashing[{}]}, 
        PlotRange  ->  All ]
Print["empirical distribution functions of both errors \
compared"]
Show[P1, P2]