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

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


k :  (2k-1) degree spline interpolation
3
d = number of digits of accuracy
30
nn = number of monte-carlo replications
100
Statistics of Hermite Error and Complete Error
Mean Hermite, Mean Complete, Ratio
{0.016227552869699047,0.002877360604248092,5.639735542962856}
Median Hermite, Median Complete
{0.0055986147927716,0.00115900616426514,4.8305306437446}
StandardDev Hermite, StandardDev Complete, Ratio
{0.056642969867939067,0.003748218814617313,
  266.7933889292155 Standard Deviation[{0.00206393521073861344654927,
        0.001777299979702536305830197,0.000686932893214305461290085,
        0.005701292648788075136095062,0.00217266751712148997002897,
        0.016200107679030361590,0.0035114658514587047473847,
        0.0011348299790105589401923,0.003001755585967085132961729,
        0.0013195835425556599691931,0.00478379464402114936265868,
        0.00712055762229376894091427,0.0013813310991816504355405,
        0.0006936336825029853607511518,0.0018758086864107822423574,
        0.0077737941771966339041494,0.013060595595921401928348,
        0.0162565093870817770207,0.0097452246537160717459221,
        0.01229398160293263908749,0.076300001533357216839,
        0.01044093066041340013830,0.0012452607815226737053568,
        0.0093352694786162783503517,0.00157205997134402122864241,
        0.00364039848305144302068452,0.013813026607105307443366,
        0.0116243697143571646441,0.0380302816296489061471,
        0.004435224368945632083373,0.0052602618595376270557321,
        0.0166373649428508046,0.016578927142007473204578,
        0.007296287747973397813181,0.0072586479830753863986641,
        0.00769121231378186179289920,0.00256462923207779471351798,
        0.01265920122703047449736323,0.004086967404529379762592,
        0.003732447298018800447365953,0.0281092774896864021484,
        0.00413308508567238047571332,0.002071767261652271758594590,
        0.01155057752747021326010201,0.006612624658455724336654875,
        0.001290039614038396539013904,0.019734477595886350552555,
        0.0011848127337561058143681,0.011117102473918737890700,
        0.00464329996937150641484972,0.00532435213775391346211782,
        0.0170640570367780934952642,0.0022719953902965001246505289,
        0.00260539620775349833695502,0.00184843864358034472615211,
        0.00529014208671024888933789,0.00493858434953057517781975,
        0.011909457515307503520612,0.09698650302871357164262912,
        0.0087238155517991292439426,0.001703642815944583871830,
        0.0018155335437054989837286161,0.012735892458489067958990,
        0.012637029036384726805212,0.0165573074368540,
        0.0028210736389263544146364,0.0026227440415397354728665,
        0.009781867301444768256617795,0.0022704799409354191471415,
        0.124641070503621674313,0.001245627628989707802142,
        0.0116745265609084534068829,0.00388641855740506523657998,
        0.03205287208132855206547231,0.003000180567346966480704651,
        0.0015562309993455613974808,0.0011728758590366784316406,
        0.0077751147616020847984201,0.013758753869573731773990,
        0.0007206583924112144726003124,0.548444345027578617743191,
        0.004369205940134128865765097,0.04262687219983609134294936,
        0.005495936936755168684284,0.002085924861655358306740803,
        0.00177262235018472548288451,0.02580248585181266975609,
        0.002333866435821479136631450,0.006426312418202877504733631,
        0.004586416872941605303311153,0.0062842672724473220842035,
        0.006247117520394510647226778,0.02096696936682296879098,
        0.020420728495428141751184470,0.011949779647573637486701778,
        0.013915490000461421521136,0.0087524486160709349924252,
        0.00213382176290624816323354,0.002096833362056181956217511,
        0.00744626325880407835501182}]}
0.95-q Hermite, 0.95-q Complete, ratio
{0.038030281629649,0.0122760590775260,3.0979226630859}
0.99-q Hermite, 0.99-q Complete, ratio
{0.124641070503622,0.0164603286270675,7.5722103323418}
Max Hermite, Max complete, ratio
{0.548444345027578617743191,0.0167213038035724,32.7991376432266}
Statistics of their ratio
Mean ratio
31.4214557685426342
Median ratio
2.9763122885411
StandardDev ratio
207.867721666748374
0.95-q ratio
48.531577919400
0.99-q ratio
279.06249532212
Max ratio
2067.87135368833107