# Maple functions for the Cramer-von Mises Statistic # # cvm(x) - distribution function for the asymptotic # distribution - returns the probability for the CVM # statistic to be less than or equal to x # # cdf(x,n) - approximate distribution function for sample # size n - returns the probability for the CVM # statistic to be less than or equal to x # # icdf(p,n) - inverse distribution function for sample # size n - returns the quantile for probability p. # # Reference: "The exact and asymptotic distributions of # Cramer-von Mises statistics" by Sandor Csorgo and Julian # Faraway appearing in the Journal of the Royal Statistical # Society, Series B (1996) vol 58 pages 221-234. # # This code written by Julian Faraway (faraway@umich.edu) D2 := proc(y) (y/2)^(3/2)*(BesselK(1/4,y^2/4)+BesselK(3/4,y^2/4))/(Pi^(1/2)); end: D3 := proc(y) local q: q := y^2/4: (y/2)^(5/2)*(2*BesselK(1/4,q)+3*BesselK(3/4,q)-BesselK(5/4,q))/(Pi^(1/2)); end: ED2 := proc(y) exp(-y^2/4)*D2(y); end: ED3 := proc(y) exp(-y^2/4)*D3(y); end: Ak := proc(k,x) (2*k+1)*GAMMA(k+1/2)*ED2((4*k+3)/(2*sqrt(x)))/(9*x^(3/4)) + GAMMA(k+1/2)*ED3((4*k+1)/(2*sqrt(x)))/(72*x^(5/4)) + 2*(2*k+3)*GAMMA(k+3/2)*ED3((4*k+5)/(2*sqrt(x)))/(12*x^(5/4)) + 7*(2*k+1)*GAMMA(k+1/2)*ED2((4*k+1)/(2*sqrt(x)))/(144*x^(3/4)) + 7*(2*k+1)*GAMMA(k+1/2)*ED2((4*k+5)/(2*sqrt(x)))/(144*x^(3/4)): end: psi1 := proc(x) local k,tot,z; k := -1: tot := 0: while k < 20 do k := k+1: z := evalf(-Ak(k,x)/(Pi*k!)); tot := tot+z; if(evalf(abs(z)) < 10.0^(-7) ) then break fi: od: RETURN(cvm(x)/12+tot); end: cvm := proc(x) local k,tot,z,q; k := -1; tot := 0; while k < 10 do k := k + 1: q := (4*k+1)^2/(16*x): z := evalf((-1)^k*binomial(-1/2,k)*sqrt(4*k+1)* exp(-q)*BesselK(1/4,q)/sqrt(x)): tot := tot + z; if(evalf(abs(z)) < 10.0^(-7) ) then break fi: od: RETURN(evalf(tot/Pi)); end: cdf := proc (x,n) cvm(x)+psi1(x)/n: end: icdf := proc(p,n) local x1,x2,fx1,fx2,xn; x1 := 1/(12*n)+0.01; x2 := x1*2; fx1 := cdf(x1,n)-p; fx2 := cdf(x2,n)-p; while abs(x1-x2) > 10.0^(-7) do xn := x2 - fx2*(x2-x1)/(fx2-fx1); if ( xn < 0) then xn := 1/(12*n) fi; x1 := x2; x2 := xn; fx1 := fx2; fx2 := cdf(x2,n)-p; od: RETURN(x2); end: