289 lines
12 KiB
FortranFixed
289 lines
12 KiB
FortranFixed
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
subroutine jbtest(nsamp,x,alpha,h,p,jbstat,cv,hb)
|
|
implicit none
|
|
|
|
double precision fn9999
|
|
parameter(fn9999=-9999.0d0)
|
|
integer nsiglev
|
|
parameter(nsiglev=17)
|
|
|
|
integer nsamp,h,hb
|
|
!h=1, reject the hypothesis of normal distribution
|
|
!h=0, fail to reject the nul hypothesis
|
|
|
|
integer i,k
|
|
integer iorder(nsiglev)
|
|
|
|
double precision x(nsamp),alpha,p,jbstat,cv,
|
|
&zscores(nsamp),skew,kurt,fstd,fmean,fmin,fmax
|
|
|
|
double precision cvs(nsiglev),alphas(nsiglev),tmpy2(nsiglev),
|
|
&alphasort(nsiglev)
|
|
|
|
data alphas/0.001, 0.0016681,0.0027826,0.0046416,0.0077426,
|
|
& 0.01, 0.012915, 0.021544, 0.025, 0.035938,
|
|
& 0.05, 0.059948, 0.1, 0.15, 0.2,
|
|
& 0.25, 0.50/
|
|
|
|
if (nsamp .lt. 2) then
|
|
write(*,*) 'stats:jbtest:NotEnoughData,
|
|
&Sample vector X must have at
|
|
&least 2 valid observations.'
|
|
h = fn9999
|
|
p = fn9999
|
|
jbstat = fn9999
|
|
cv = fn9999
|
|
hb=fn9999
|
|
|
|
elseif (nsamp .eq. 2) then
|
|
write(*,*) 'The J-B stat is a constant
|
|
&when the sample size is 2'
|
|
h = 0
|
|
p = 1.0
|
|
jbstat = 1.0/3.0
|
|
cv = fn9999
|
|
hb = 1
|
|
else
|
|
|
|
skew = 0.0
|
|
kurt = 0.0
|
|
|
|
call stdmaxmeanmin(nsamp,x,fstd,fmean,fmin,fmax)
|
|
|
|
! readjust sample size from (n-1) to n according to matlab
|
|
|
|
fstd = dsqrt(fstd*fstd*dble(nsamp-1)/dble(nsamp))
|
|
do i=1,nsamp
|
|
zscores(i) = (x(i)-fmean)/fstd
|
|
enddo
|
|
|
|
do i=1,nsamp
|
|
skew = skew + zscores(i)*zscores(i)*zscores(i)
|
|
kurt = kurt + zscores(i)*zscores(i)*zscores(i)*zscores(i)
|
|
enddo
|
|
|
|
skew = skew/dble(nsamp)
|
|
kurt = kurt/dble(nsamp) - 3
|
|
jbstat = nsamp*(skew*skew/6 + kurt*kurt/24)
|
|
|
|
|
|
! Get a row of the critical value table for the current sample size
|
|
call cvtbl(nsamp,cvs)
|
|
|
|
! 1-D interpolation into the tabulated quantiles.
|
|
call SPLINE(alphas,cvs,nsiglev,1.0d+31,1.0d+31,tmpy2)
|
|
call SPLINT(alphas,cvs,tmpy2,nsiglev,alpha,cv)
|
|
|
|
! Compute the P-value. Warn if the P-value is not found within the
|
|
! available 'alphas' of the table and return one of the extremes.
|
|
|
|
if (jbstat < cvs(nsiglev)) then ! smallest critval at end
|
|
write(*,*) 'P is greater than the largest tabulated value'
|
|
p = alphas(nsiglev)
|
|
hb=1
|
|
elseif (cvs(1) < jbstat) then ! largest critval at beginning
|
|
write(*,*) 'P is less than the smallest tabulated value'
|
|
p = alphas(1)
|
|
hb=1
|
|
else
|
|
|
|
call sort_shell(nsiglev,cvs,iorder)
|
|
do k=1,nsiglev
|
|
alphasort(k) = alphas(iorder(k))
|
|
enddo
|
|
call SPLINE(cvs,alphasort,nsiglev,1.0d+31,1.0d+31,tmpy2)
|
|
call SPLINT(cvs,alphasort,tmpy2,nsiglev,jbstat,p)
|
|
|
|
hb=0
|
|
if(jbstat .gt. cv) then
|
|
h=1
|
|
else
|
|
h=0
|
|
endif
|
|
|
|
endif
|
|
endif
|
|
|
|
return
|
|
end
|
|
|
|
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
subroutine cvtbl(n,cvs)
|
|
implicit none
|
|
|
|
integer n,nsiglev,nsamples,i,j,k
|
|
parameter(nsiglev=17,nsamples=32)
|
|
integer sampleSizes(nsamples) !Sample sizes
|
|
double precision criticalValues(nsamples,nsiglev) !Critical value
|
|
double precision cvs(nsiglev)
|
|
|
|
double precision x(nsamples),y(nsamples),y2(nsamples),
|
|
& ysort(nsamples),tmp(nsamples)
|
|
integer iorder(nsamples)
|
|
|
|
data sampleSizes/3, 4, 5, 10, 15, 20, 25, 30, 35,
|
|
& 40, 45, 50, 60, 70, 80, 90, 100, 125,
|
|
& 150, 175, 200, 250, 300, 400, 500, 800, 1000,
|
|
& 1200,1400,1600,1800,2000/
|
|
|
|
data criticalValues/ 0.5312, 0.9606, 1.8289,10.9719,
|
|
& 19.5425,25.0722,28.4885,30.6106,
|
|
& 31.9343,32.7514,33.2273,33.4825,
|
|
& 33.5009,33.2610,32.8742,32.3418,
|
|
& 31.8884,30.6089,29.4290,28.4225,
|
|
& 27.5140,26.0239,24.8353,23.0771,
|
|
& 21.8170,19.5539,18.6707,18.0126,
|
|
& 17.5190,17.1305,16.8158,16.5585,
|
|
& 0.5312, 0.9590, 1.8053, 9.8430,
|
|
& 16.7207,21.0001,23.6332,25.2501,
|
|
& 26.2695,26.9025,27.2686,27.4765,
|
|
& 27.5458,27.3738,27.0975,26.7278,
|
|
& 26.3990,25.4748,24.6043,23.8582,
|
|
& 23.1914,22.0764,21.1758,19.8542,
|
|
& 18.9068,17.1748,16.4879,15.9830,
|
|
& 15.6036,15.3029,15.0600,14.8633,
|
|
& 0.5312, 0.9564, 1.7727, 8.6751,
|
|
& 14.0454,17.2981,19.2818,20.5242,
|
|
& 21.3105,21.8016,22.1030,22.2863,
|
|
& 22.3790,22.2848,22.1135,21.8770,
|
|
& 21.6407,21.0100,20.3979,19.8658,
|
|
& 19.3922,18.5877,17.9384,16.9669,
|
|
& 16.2666,14.9865,14.4745,14.0990,
|
|
& 13.8149,13.5883,13.4102,13.2657,
|
|
& 0.5312, 0.9520, 1.7276, 7.4815,
|
|
& 11.5527,13.9762,15.4635,16.4088,
|
|
& 17.0176,17.4145,17.6594,17.8172,
|
|
& 17.9431,17.9209,17.8305,17.6977,
|
|
& 17.5510,17.1460,16.7425,16.3846,
|
|
& 16.0635,15.5123,15.0633,14.3843,
|
|
& 13.8903,12.9765,12.6143,12.3493,
|
|
& 12.1501,11.9883,11.8682,11.7666,
|
|
& 0.5312, 0.9448, 1.6661, 6.2927,
|
|
& 9.2814,11.0602,12.1658,12.8805,
|
|
& 13.3572,13.6749,13.8848,14.0345,
|
|
& 14.1788,14.2092,14.1911,14.1283,
|
|
& 14.0491,13.8271,13.5878,13.3664,
|
|
& 13.1643,12.8124,12.5206,12.0743,
|
|
& 11.7469,11.1441,10.9072,10.7362,
|
|
& 10.6072,10.5059,10.4296,10.3636,
|
|
& 0.5312, 0.9396, 1.6275, 5.7077,
|
|
& 8.2365, 9.7531,10.7058,11.3263,
|
|
& 11.7488,12.0355,12.2302,12.3739,
|
|
& 12.5255,12.5801,12.5841,12.5542,
|
|
& 12.5067,12.3565,12.1804,12.0139,
|
|
& 11.8602,11.5923,11.3653,11.0157,
|
|
& 10.7604,10.2914,10.1106, 9.9803,
|
|
& 9.8811, 9.8072, 9.7478, 9.6981,
|
|
& 0.5311, 0.9329, 1.5828, 5.1350,
|
|
& 7.2613, 8.5489, 9.3658, 9.9063,
|
|
& 10.2818,10.5399,10.7201,10.8586,
|
|
& 11.0182,11.0884,11.1155,11.1112,
|
|
& 11.0882,10.9976,10.8779,10.7589,
|
|
& 10.6494,10.4542,10.2839,10.0226,
|
|
& 9.8318, 9.4852, 9.3521, 9.2587,
|
|
& 9.1869, 9.1333, 9.0911, 9.0549,
|
|
& 0.5310, 0.9133, 1.4723, 4.0456,
|
|
& 5.5238, 6.4401, 7.0389, 7.4472,
|
|
& 7.7395, 7.9520, 8.1098, 8.2340,
|
|
& 8.3968, 8.4921, 8.5520, 8.5850,
|
|
& 8.6004, 8.6020, 8.5719, 8.5307,
|
|
& 8.4904, 8.4114, 8.3384, 8.2266,
|
|
& 8.1454, 8.0027, 7.9501, 7.9127,
|
|
& 7.8828, 7.8611, 7.8446, 7.8296,
|
|
& 0.5309,0.9056,1.4343,3.7474,
|
|
& 5.0729,5.8998,6.4463,6.8222,
|
|
& 7.0942,7.2949,7.4469,7.5661,
|
|
& 7.7287,7.8286,7.8938,7.9361,
|
|
& 7.9589,7.9816,7.9718,7.9496,
|
|
& 7.9259,7.8751,7.8276,7.7526,
|
|
& 7.6984,7.6040,7.5691,7.5441,
|
|
& 7.5239,7.5089,7.4979,7.4872,
|
|
& 0.5305,0.8817,1.3294,3.0691,
|
|
& 4.0773,4.7176,5.1501,5.4578,
|
|
& 5.6856,5.8598,5.9945,6.1045,
|
|
& 6.2620,6.3696,6.4462,6.5028,
|
|
& 6.5429,6.6071,6.6387,6.6571,
|
|
& 6.6688,6.6803,6.6842,6.6877,
|
|
& 6.6882,6.6858,6.6847,6.6832,
|
|
& 6.6800,6.6789,6.6777,6.6763,
|
|
& 0.5297,0.8519,1.2185,2.5239,
|
|
& 3.2985,3.8011,4.1494,4.4039,
|
|
& 4.5973,4.7481,4.8689,4.9697,
|
|
& 5.1203,5.2305,5.3134,5.3796,
|
|
& 5.4314,5.5277,5.5919,5.6408,
|
|
& 5.6783,5.7343,5.7728,5.8248,
|
|
& 5.8581,5.9096,5.9282,5.9408,
|
|
& 5.9482,5.9546,5.9600,5.9635,
|
|
& 0.5290,0.8316,1.1516,2.2555,
|
|
& 2.9215,3.3596,3.6684,3.8968,
|
|
& 4.0734,4.2126,4.3258,4.4214,
|
|
& 4.5688,4.6803,4.7662,4.8378,
|
|
& 4.8957,5.0071,5.0863,5.1474,
|
|
& 5.1957,5.2681,5.3194,5.3889,
|
|
& 5.4336,5.5043,5.5296,5.5471,
|
|
& 5.5587,5.5677,5.5755,5.5806,
|
|
& 0.5251,0.7553,0.9442,1.6231,
|
|
& 2.0533,2.3505,2.5707,2.7431,
|
|
& 2.8827,2.9987,3.0973,3.1834,
|
|
& 3.3246,3.4374,3.5292,3.6071,
|
|
& 3.6730,3.8025,3.8987,3.9732,
|
|
& 4.0327,4.1224,4.1873,4.2748,
|
|
& 4.3320,4.4246,4.4580,4.4814,
|
|
& 4.4979,4.5105,4.5214,4.5289,
|
|
& 0.5176,0.6721,0.7945,1.2821,
|
|
& 1.5965,1.8239,1.9986,2.1390,
|
|
& 2.2547,2.3524,2.4361,2.5097,
|
|
& 2.6316,2.7298,2.8104,2.8788,
|
|
& 2.9372,3.0523,3.1384,3.2050,
|
|
& 3.2584,3.3402,3.3988,3.4790,
|
|
& 3.5318,3.6180,3.6500,3.6718,
|
|
& 3.6876,3.6997,3.7101,3.7173,
|
|
& 0.5074,0.6303,0.7302,1.1235,
|
|
& 1.3779,1.5631,1.7063,1.8216,
|
|
& 1.9172,1.9980,2.0674,2.1283,
|
|
& 2.2297,2.3115,2.3789,2.4361,
|
|
& 2.4851,2.5819,2.6543,2.7106,
|
|
& 2.7559,2.8250,2.8749,2.9434,
|
|
& 2.9886,3.0631,3.0907,3.1099,
|
|
& 3.1237,3.1345,3.1433,3.1499,
|
|
& 0.4946,0.5947,0.6878,1.0198,
|
|
& 1.2336,1.3885,1.5079,1.6040,
|
|
& 1.6835,1.7508,1.8085,1.8592,
|
|
& 1.9434,2.0114,2.0674,2.1151,
|
|
& 2.1557,2.2363,2.2964,2.3436,
|
|
& 2.3812,2.4388,2.4808,2.5382,
|
|
& 2.5760,2.6388,2.6624,2.6787,
|
|
& 2.6904,2.6996,2.7072,2.7129,
|
|
& 0.4063,0.4739,0.5285,0.6951,
|
|
& 0.7916,0.8577,0.9071,0.9457,
|
|
& 0.9771,1.0033,1.0256,1.0449,
|
|
& 1.0768,1.1023,1.1231,1.1408,
|
|
& 1.1557,1.1852,1.2072,1.2243,
|
|
& 1.2382,1.2591,1.2744,1.2956,
|
|
& 1.3097,1.3334,1.3421,1.3484,
|
|
& 1.3530,1.3564,1.3595,1.3619/
|
|
|
|
|
|
! Interpolate a row of critical values for the given sample size.
|
|
do i=1,nsiglev
|
|
do j=1,nsamples
|
|
x(j) = 1.0/dble(sampleSizes(j))
|
|
y(j) = criticalValues(j,i)
|
|
enddo
|
|
|
|
call sort_shell(nsamples,x,iorder)
|
|
do k=1,nsamples
|
|
ysort(k) = y(iorder(k))
|
|
enddo
|
|
|
|
call SPLINE(x,ysort,nsamples,1.0d+31,1.0d+31,y2)
|
|
call SPLINT(x,ysort,y2,nsamples,1.0/dble(n),cvs(i))
|
|
|
|
enddo
|
|
|
|
return
|
|
end
|
|
|
|
|
|
|