PROGRAM se_logistic ! Generate artificial data to see how much bias there is in ! the estimated standard errors. USE Logistic_Regression IMPLICIT NONE REAL (dp) :: chisq, devnce, beta(0:1), se_beta(0:1), p, z(10), dev, & aver(0:1) = 0.0_dp, sumbb(0:1) = 0.0_dp, aver_se(0:1) = 0.0_dp INTEGER :: k, i, case, n(5), s(5), ngroups = 5, ndf, ier INTEGER, ALLOCATABLE :: seed(:) REAL (dp), PARAMETER :: x(5,1) = RESHAPE( & (/ -2.1972_dp, -0.8473_dp, 0.0_dp, & 0.8473_dp, 2.1972_dp /), (/ 5, 1 /) ) CALL RANDOM_SEED(SIZE=k) ALLOCATE( seed(k) ) WRITE(*, '(a, i2, a)')' Enter ', k, ' integers for random no. seeds: ' READ(*, *) seed WRITE(*, '(a, (7i10))') ' Random no. seeds: ', seed CALL RANDOM_SEED(PUT=seed) n = 10 DO case = 1, 10000 ! Generate binomial samples from samples of size 10, ! with p = 0.1, 0.3, 0.5, 0.7 & 0.9 DO i = 1, 5 p = 0.1_dp * (2*i - 1) CALL RANDOM_NUMBER( z ) s(i) = COUNT( z < p ) END DO CALL logistic(ngroups, x, 1, s, n, chisq, devnce, ndf, beta, se_beta, ier) ! Update the average dev = beta(0) - aver(0) aver(0) = aver(0) + dev/case sumbb(0) = sumbb(0) + dev*(beta(0) - aver(0)) dev = beta(1) - aver(1) aver(1) = aver(1) + dev/case sumbb(1) = sumbb(1) + dev*(beta(1) - aver(1)) aver_se = aver_se + (se_beta - aver_se) / case END DO sumbb = SQRT( sumbb / 9999. ) WRITE(*, '(a, 2f9.5)') ' Average estimates of beta = ', aver WRITE(*, '(a, 2f9.5)') ' Std. devn. of sample beta = ', sumbb WRITE(*, '(a, 2f9.5)') ' Average se_beta = ', aver_se STOP END PROGRAM se_logistic