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

 (*    Complete and Hermite spline interpolation of Shadrin' 
        s null spline phi              *) 
 (*    on [0, 1] with 2  k - 
        4 interior knots                                                   *) 
 (*    Monte - 
        carlo program without summary plots                                           *) 
 (*    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


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)
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}]]
BasisNull[x_, vec_List, k_] := 
      Join[{1}, Table[B1[x, j, 0], {j, 1, 2*k - 2}], 
        Table[B2[x, vec[[j]], k, 1]/(2*k - 1), {j, 1, 
            2*k - 4}]]
BasisNullDer[x_, l_, vec_List, k_] := 
      Join[{l!}, Table[B1[x, j, l], {j, l + 1, 2*k - 2}], 
        Table[B2[x, vec[[j]], k, l + 1]/(2*k - 1), {j, 1, 
            2*k - 4}]]

CompInterp[x_, coeff_List, vec_List, k_] := 
      coeff . Basis[x, vec, k]
NullSpline[x_, coeff_List, vec_List, k_] := 
      coeff . BasisNull[x, vec, k]
NullSplDer[x_, l_, k_, vec_List, coeff_List] := 
      coeff[[Table[j, {j, l + 1, 4*k - 5}]]] . 
        BasisNullDer[x, l, vec, k]
HermiteInterp[x_, coeff_List, vec_List, k_] := 
      coeff . Basis[x, vec, k]
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}]]
VectNull[k_, l_, vec_List] := 
      Join[Table[0, {m, 1, l}], Table[B1[1, i, l], {i, l, 2*k - 2}], 
        Table[B2[1, vec[[i]], k, l + 1]/(2*k - 1), {i, 1, 
            2*k - 4}]]

mxlistComp[coeff1_List, coeff2_List, vec_List, k_] := 
      Table[Abs[
          CompInterp[i/100, coeff1, vec, k] - 
            NullSplDer[i/100, k - 1, k, vec, coeff2]], {i, 1, 
          100}]
mxlistHermite[coeff1_List, coeff2_List, vec_List, k_] := 
      Table[Abs[
          HermiteInterp[i/100, coeff1, vec, k] - 
            NullSplDer[i/100, k - 1, k, vec, coeff2]], {i, 1, 
          100}]

Error[
        k_] := {
VectauUnif = 
          Sort[
            Table[Random[Real, {0, 1}, d], {j, 1, 
                2*k - 4}]]; 
Vectau = 
          Append[Prepend[VectauUnif, 0], 1]; 
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}]]; 
 
        MM =  Join[MM1, MM2, MM3]; 
MM1Null =  
          AppendRows[DiagonalMatrix[Table[j!, {j, 0, (k - 2)}]], 
            Table[0, {i, 1, k - 1}, {j, 1, 
                3*k - 4}]]; 
MM2Null =  
          Table[Join[{1}, 
              Table[B1[VectauUnif[[j]], i, 0], {i, 1, 2*k - 2}], 
              Table[B2[VectauUnif[[j]], VectauUnif[[i]], k, 
                    1]/(2*k - 1), {i, 1, 2*k - 4}]], {j, 1, 
              2*k - 4}]; 
Vec3Null =  
          Join[Table[1, {m, 1, 2*k - 1}],  
            Table[B2[1, VectauUnif[[i]], k, 1]/(2*k - 1), {i, 1, 
                2*k - 4}]]; 
MM3Null =  
          Join[{Vec3Null}, 
            Table[VectNull[k, l, VectauUnif], {l, 1, 
                k - 1}]]; 
MMNull =  
          Join[MM1Null, MM2Null, MM3Null]; 
 
        VecCNull = 
          Join[Table[
              0, {l, 1, 
                4*k - 6}], {(k - 
                    1)!}]; 
NullSplCoeff = 
          LinearSolve[MMNull, VecCNull]; 
VecPhi = 
          Join[Table[
              NullSplDer[0, l, k, VectauUnif, NullSplCoeff], {l, k - 1, 
                2*k - 3}],  {NullSplCoeff[[2*k - 1]]*(2*k - 
                      2)!}, 
            Table[NullSplDer[VectauUnif[[j]],  k - 1, k, VectauUnif,  
                NullSplCoeff], {j, 1, 2*k - 4}],  
            Table[NullSplDer[1, l, k, VectauUnif, NullSplCoeff], {l, k - 1, 
                2*k - 2}]]; 
 
        CpltSplineCoeff = 
          LinearSolve[MM, 
            VecPhi]; 

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[
              NullSplDer[Vectau[[j]], k - 1, k, VectauUnif,  
                NullSplCoeff], {j, 1, 2*k - 2}], 
            Table[NullSplDer[Vectau[[j]],  k, k, VectauUnif,  
                NullSplCoeff], {j, 1, 
                2*k - 2}]]; 
HermiteCoeff = 
          LinearSolve[MMHermite, 
            VecHermite]; 
maxComp  =  
          Max[mxlistComp[CpltSplineCoeff, NullSplCoeff, VectauUnif, k]]/
            Sum[Abs[NullSplCoeff[[j]]]*(2*k - 2)!, {j, 2*k, 
                4*k - 5}]; 
maxHerm =  
          Max[mxlistHermite[HermiteCoeff, NullSplCoeff, VectauUnif, k]]/
            Sum[Abs[NullSplCoeff[[j]]]*(2*k - 2)!, {j, 2*k, 
                4*k - 5}]; 
maxCompHerm = {maxHerm, 
            maxComp}
}

Ert = 
        Table[Flatten[Error[k]], {j, 1, nn}];
ErtH = 
        Table[Ert[[i, 1]], {i, 1, nn}];
ErtC = 
        Table[Ert[[i, 2]], {i, 1, nn}];
Print["Statistics of Hermite Error and Complete Error"]
Print["Mean Hermite, Mean Complete, Ratio"]
{Mean[ErtH], Mean[ErtC],  
      Mean[ErtH]/Mean[ErtC ]}
Print["Median Hermite, Median Complete"]
{Median[ErtH], Median[ErtC],  
      Median[ErtH]/Median[ErtC]}
Print["StandardDev Hermite, StandardDev Complete, Ratio"]
{StandardDeviation[ErtH], StandardDeviation[ErtC], 
      Standard Deviation[ErtH]/
          StandardDeviation[ErtC]}
Print["0.95-q Hermite, 0.95-q Complete, ratio"]
{Quantile[ErtH, 0.95], Quantile[ErtC, 0.95], 
      Quantile[ErtH, 0.95]/Quantile[ErtC, 0.95]}
Print["0.99-q Hermite, 0.99-q Complete, ratio"]
{Quantile[ErtH, 0.99], Quantile[ErtC, 0.99], 
      Quantile[ErtH, 0.99]/Quantile[ErtC, 0.99]}
Print["Max Hermite, Max complete, ratio"]
{Max[ErtH], Max[ErtC], 
      Max[ErtH]/Max[ErtC]}