SUBROUTINE pearsn(xbar, sd, itype, rb1, b2, bndry, sigpt, ifault) ! Code converted using TO_F90 by Alan Miller ! Date: 1999-09-25 Time: 16:22:51 ! N.B. Array a has been removed from the list of arguments, ! and replaced with a PARAMETER statement in this routine. ! ALGORITHM AS 192 APPL. STATIST. (1983) VOL.32, NO.3 ! Computes approximate significance points of a Pearson curve with given ! first four moments, or first three moments and left or right boundary. ! Error codes ! IFAULT = 0 successful completion ! = 1 SD < 0.0 ! = 2 ITYPE < 1 or > 3 ! = 3 | RB1 | > 2.0 ! = 4 XBAR value impossible with the value entered for BNDRY ! = 5 B2 cannot be computed for the first 3 moments + left boundary ! = 6 B2 out of range, that is: ! If | RB1 | <= 1.0 then ! B2 < 1.5 * | RB1 | + 1.5 or B2 > 0.2 * | RB1 | + 10.8 ! If | RB1 | > 1.0 then ! B2 < 3.9 * | RB1 | - 0.9 or B2 > 4.8 * | RB1 | + 6.2 IMPLICIT NONE REAL, INTENT(IN) :: xbar REAL, INTENT(IN) :: sd INTEGER, INTENT(IN) :: itype REAL, INTENT(IN) :: rb1 REAL, INTENT(IN OUT) :: b2 REAL, INTENT(IN) :: bndry REAL, INTENT(OUT) :: sigpt(:) INTEGER, INTENT(OUT) :: ifault REAL :: b(9), rbeta1, ca, ca2, carb1, denom, absca, bnd1, bnd2, sgn, & pi1, pi2 INTEGER :: i, j, k, l REAL, PARAMETER :: zero = 0.0, point2 = 0.2, rthalf = 0.70710678, & point9 = 0.9, one = 1.0, onept5 = 1.5, two = 2.0, & three = 3.0, thrpt9 = 3.9, four = 4.0, forpt8 = 4.8, & sixpt2 = 6.2, tenpt8 = 10.8 REAL, PARAMETER :: a(418) = (/ & -1.9336, -1.6061, 2.6955, -1.7036, 1.7236, -1.3209, -2.5464, 1.6812, -0.11812, 0.050875, 0.16042, -1.0616, 0.43929, & -0.35799, 0.51040, 1.0958, -0.63453, 0.053175, -0.019945, 0.25597, 0.32520, -0.094571, 2.0418, -2.5786, 0.61148, -0.64242, & 1.1968, -0.33145, 0.015508, -4.1444, 0.13884, 4.6625, 0.37661, -0.20447, -1.9508, -0.20214, 0.11921, -0.0085268, & -1.6453, -1.8494, 2.3915, -1.5844, 1.8290, -1.1167, -3.6091, 2.7150, -0.45857, 0.048102, 0.89117, -1.2700, 0.37653, & -0.81542, 0.56815, 1.6562, -1.1908, 0.24574, -0.024375, -1.5063, 4.4876, -0.60765, -6.5584, 2.8944, -0.42381, 2.2664, & -1.3425, 0.25806, -0.010421, -1.9526, 0.21332, 2.1317, -1.1996, 0.21033, -0.52154, 0.53597, -0.12355, 0.0060658, & -1.1044, -1.1300, 1.7681, 0.10947, 0.30566, -0.90598, -0.98814, 0.49355, 0.30625, -0.018707, 0.74941, -1.3174, 0.23974, & -0.12433, 0.59105, 0.85714, -0.47292, -0.16028, 0.012684, -1.4645, 4.5349, -0.71053, -7.9213, 4.0018, -0.61932, 2.7320, & -1.7752, 0.35919, -0.014944, -2.0999, 0.30754, 3.6865, -2.3140, 0.39961, -1.1530, 1.0257, -0.22981, 0.011032, & -0.85842, 0.78929, 1.2196, -0.55088, -0.96385, -0.57312, -0.076974, 0.34882, 0.47311, -0.056120, -0.63153, -1.2697, 0.64649, & 0.92712, 0.49088, 0.48118, -0.46407, -0.39753, 0.052413, -0.92873, 2.2039, 0.25828, 0.74312, -2.4307, 0.49934, 0.029833, & 0.68249, -0.19687, 0.0068136, -2.9758, -0.050738, 0.067069, 2.0835, -0.45211, -0.28149, -0.58012, 0.17621, -0.0067195, & -2.4031, -1.4891, 0.75590, -2.3309, 9.2314, -3.8930, -7.0192, 2.3648, 0.44308, -0.029437, 16.616, -4.9826, 1.0626, & -18.888, 7.2830, 14.642, -2.6394, -1.2739, 0.060862, -1.6092, 0.81385, -0.059847, -0.18180, 1.2679, -0.36205, 0.29683, & -0.83310, 0.26630, -0.015081, 0.43536, 0.27965, 1.2516, -3.0536, 0.65936, -1.0579, 1.7886, -0.47962, 0.025311, & -2.0711E-4, 0.0068277, -7.5482E-6, 0.47086, -0.21213, 1.1003E-4, -0.22462, 0.11733, -0.079682, -1.0565E-5, 2.5087, -2.6169, & 0.94103, -2.3261, 1.6482, 2.3024, -1.4821, -0.020294, 9.9168E-4, -0.097078, 0.58192, -0.019936, -0.7025, 0.11526, 0.0013959, & 0.1595, -0.062382, 0.01338, 1.8557E-4, -1.8383, -0.54124, 0.68041, 1.7654, -0.37606, -1.113, 0.044447, 0.042456, -0.0017448, & 0.39672, -0.17140, -0.74246, 0.48553, 0.19650, 0.58942, -1.5768, 0.37082, -0.32481, 0.025296, -0.60818, -1.9005, 0.45451, & 0.79612, 1.0544, -1.9010, 0.29849, -0.51728, 0.050926, 0.45223, -2.4004, 0.38372, 3.6817, -1.0796, 0.016619, -1.6829, & 0.81147, -0.10376, 0.0052116, -4.3837, 0.67456, 6.3029, -1.8042, 0.023516, -2.9320, 1.4537, -0.20680, 0.010634, & 0.77212, 0.028272, -1.3932, 0.95628, -0.37918, 0.85487, 0.38225, -1.3009, 0.42400, -0.051933, 0.17133, -1.5330, 1.2721, & -0.58446, 0.78993, 0.33412, -1.2988, 0.39155, -0.046683, 2.2548, -1.7005, -0.65377, -2.9256, 2.1397, 0.093244, 3.2094, & -2.2896, 0.24580, -0.0020536, 0.30478, -0.58544, -4.0524, 1.8239, 0.083952, 3.4072, -2.0163, 0.20345, -7.5748E-4, & 1.1504, 0.22674, -1.7407, 1.1356, -0.58596, 0.88839, 0.34743, -1.2531, 0.42319, -0.014362, 0.17074, -1.2892, 1.1653, & -0.46699, 0.58479, 0.33975, -1.0005, 0.25271, -0.0091547, 2.1265, -8.7599, 2.3560, 8.0491, -3.3844, -0.028416, -4.3522, & 3.9741, -0.92840, 0.055908, -5.5206, 1.6024, 6.1077, -2.4508, -0.012514, -3.7719, 2.8877, -0.59134, 0.034299, & 1.5989, 1.5353, -2.2124, 0.83145, -1.2543, 0.92247, 2.2076, -1.4147, 0.053879, 0.058513, 0.42581, -1.1636, 0.67165, & -0.47849, 0.46159, 0.97601, -0.61165, -0.054400, 0.029952, 10.803, -5.6739, -8.3559, -12.404, 29.846, -7.7439, -6.2673, & -1.3612, 1.3966, -0.22122, 2.8678, -3.0580, -14.864, 16.489, -4.0334, -0.49139, -1.7927, 0.98545, -0.10118, & 2.4201, -1.9281, -3.3357, 2.3720, 1.7318, 1.4149, -0.64616, -2.7558, 0.28474, -0.034471, -0.078628, -1.2092, 0.43924, & 0.43093, 0.53223, 0.34235, -1.0741, 0.080494, -0.013004, -15.787, -3.9798, 23.933, 24.332, -46.762, 6.0862, 15.874, & 5.2360, -2.4644, 0.28404, -14.830, 8.5161, 23.701, -19.419, 2.4239, -1.8457, 4.8007, -1.2525, 0.099997 /) ! Check for invalid input arguments ifault = 1 IF (sd <= zero) RETURN ifault = 2 IF (itype < 1 .OR. itype > 3) RETURN rbeta1 = ABS(rb1) ifault = 3 IF (rbeta1 > two) RETURN IF (itype == 1) GO TO 10 ! Compute B2 using XBAR, SD, RB1, and the known boundary ifault = 4 IF ((itype == 2 .AND. xbar <= bndry) .OR. (itype == 3 .AND. xbar >= bndry)) & RETURN ifault = 5 ca = sd / (xbar - bndry) ca2 = ca * ca carb1 = ca * rb1 denom = four * ca2 - carb1 + two IF (denom <= zero) RETURN absca = ABS(ca) bnd1 = (ca2 - one) / ca IF (absca <= rthalf) THEN bnd2 = four * ca / (one - ca2) ELSE bnd2 = four * ca + two / ca END IF IF ((itype == 2 .AND. (rb1 < bnd1 .OR. rb1 > bnd2)) .OR. & (itype == 3 .AND. (rb1 < bnd2 .OR. rb1 > bnd1))) RETURN b2 = three * (carb1 * (carb1 + one) + rb1 * rb1 + two) / denom 10 ifault = 6 IF ((rbeta1 <= one .AND. (b2 < onept5 * rbeta1 + onept5 .OR. & b2 > point2 * rbeta1 + tenpt8)) .OR. (rbeta1 > one .AND. & (b2 < thrpt9 * rbeta1 - point9 .OR. b2 > forpt8 * rbeta1 + sixpt2))) RETURN ifault = 0 ! Initialization b(1) = rbeta1 b(2) = b2 b(3) = rbeta1 * rbeta1 b(4) = rbeta1 * b2 b(5) = b2 * b2 b(6) = b(3) * rbeta1 b(7) = b(3) * b2 b(8) = b(5) * rbeta1 b(9) = b(5) * b2 ! Significance point computation k = -19 IF (rbeta1 > one) k = 0 sgn = one IF (rb1 < zero) sgn = -one DO i = 1, 11 l = i IF (rb1 < zero) l = 12 - i k = k + 20 pi1 = a(k) DO j = 1, 9 k = k + 1 pi1 = pi1 + a(k) * b(j) END DO pi2 = one DO j = 1, 9 k = k + 1 pi2 = pi2 + a(k) * b(j) END DO ! Compute significance point sigpt(l) = xbar + sd * sgn * pi1 / pi2 END DO RETURN END SUBROUTINE pearsn PROGRAM test_as192 IMPLICIT NONE INTERFACE SUBROUTINE pearsn(xbar, sd, itype, rb1, b2, bndry, sigpt, ifault) IMPLICIT NONE REAL, INTENT(IN) :: xbar REAL, INTENT(IN) :: sd INTEGER, INTENT(IN) :: itype REAL, INTENT(IN) :: rb1 REAL, INTENT(IN OUT) :: b2 REAL, INTENT(IN) :: bndry REAL, INTENT(OUT) :: sigpt(:) INTEGER, INTENT(OUT) :: ifault END SUBROUTINE pearsn END INTERFACE INTEGER :: itype, ifault REAL :: xbar, sd, rb1, b2, bndry, sigpt(11) DO WRITE(*, '(a)', ADVANCE='NO') ' Enter mean & st.devn.: ' READ(*, *) xbar, sd WRITE(*, *) 'ITYPE = 1 4 moments used' WRITE(*, *) 'ITYPE = 2 3 moments + left boundary used' WRITE(*, *) 'ITYPE = 3 3 moments + right boundary used' WRITE(*, '(a)', ADVANCE='NO') ' Enter ITYPE: ' READ(*, *) itype WRITE(*, '(a)', ADVANCE='NO') ' Enter root(beta1) with sign of beta1: ' READ(*, *) rb1 IF (itype == 1) THEN WRITE(*, '(a)', ADVANCE='NO') ' Enter beta2: ' READ(*, *) b2 END IF IF (itype > 1) THEN WRITE(*, '(a)', ADVANCE='NO') ' Enter boundary: ' READ(*, *) bndry END IF CALL pearsn(xbar, sd, itype, rb1, b2, bndry, sigpt, ifault) IF (ifault /= 0) THEN WRITE(*, *) 'IFAULT = ', ifault ELSE WRITE(*, '(" Percentage points"/ " ", 11f7.2)') sigpt END IF END DO STOP END PROGRAM test_as192