Replacing IMSL call ?
Author Message Replacing IMSL call ?

Hi - I am supporting a Win32 Console Application written with Microsoft
fortran Powerstation.  It makes 1 IMSL call to "eqtil".

It reads numbers from a file, crunches them furiously, then writes out
several files as an analysis of the data read, including one file with
any error messages for the user.

This compiler was some-what buggy, and sooo much is available now that
wasn't when this program was written.

I am wondering if anyone knows if there is free-ware fortran code out
there now that replaces eqtil ?  (see manual on eqtil below)

Thanks - Gary

According to the manual:

EQTIL/DEQTIL (Single/Double precision)
Compute empirical quantiles.
Usage
CALL EQTIL (NOBS, X, NQPROP, QPROP, Q, XLO, XHI, NMISS)
Arguments
NOBS Number of observations. (Input)
NOBS must be greater than or equal to one.
X Vector of length NOBS containing the data. (Input)
NQPROP Number of quantiles. (Input)
NQPROP must be greater than or equal to one.
QPROP Vector of length NQPROP containing the quantile proportions.
(Input)
The elements of QPROP must lie in the interval (0, 1).
Q Vector of length NQPROP containing the empirical quantiles. (Output)
Q(i) corresponds to the empirical quantile at proportion QPROP(i). The
quantiles
are determined by linear interpolation between adjacent ordered sample
values.
XLO Vector of length NQPROP containing the largest element of X less
than or
equal to the desired quantile. (Output)
XHI Vector of length NQPROP containing the smallest element of X greater
than or equal to the desired quantile. (Output)
NMISS Number of missing values. (Output)
1. Automatic workspace is allocated only if X is not sorted on input. The
amount allocated is
EQTIL NOBS units, or
DEQTIL 2 * NOBS units.
Workspace may be explicitly provided, if desired, by use of
E2TIL/DE2TIL. The reference is
CALL E2TIL (NOBS, X, NQPROP, QPROP, Q, XLO, XHI,
NMISS, WK)
36 Chapter 1: Basic Statistics IMSL STAT/LIBRARY
WK Workspace of length NOBS containing the sorted data. (Output)
If X is sorted in ascending order with all missing values at the end of X,
then X and WK may share the same storage location.
2. Informational error
Type Code
3 1 All of the observations are missing values. The
elements of Q, XLO, and XHI have been set to NaN
(not a number).
3. Missing values (NaN) are excluded from the analysis. Empirical
quantiles are based on the NOBS - NMISS nonmissing elements of X.
Algorithm
The routine EQTIL determines the empirical quantiles, as indicated in
the vector
QPROP, from the data in X. The routine EQTIL first checks to see if X is
sorted; if
X is not sorted, the routine does either a complete or partial sort,
depending on
how many order statistics are required to compute the quantiles requested.
The routine EQTIL returns the empirical quantiles and, for each
quantile, the two
order statistics from the sample that are at least as large and at least
as small as
the quantile. For a sample of size n, the quantile corresponding to the
proportion
p is defined as
Q(p) = (1 - f )x(M) + f x(M + 1)
where j = ?p(n + 1)?, f = p(n + 1) - j, and x(M) is the j-th order
statistic, if
1 j < n; otherwise, the empirical quantile is the smallest or largest
order statistic.
Example
In this example, five empirical quantiles from a sample of size 30 are
obtained.
Notice that the 0.5 quantile corresponds to the sample median. The data
are from
Hinkley (1977) and Velleman and Hoaglin (1981). They are the
measurements (in
inches) of precipitation in Minneapolis/St. Paul during the month of
March for 30
consecutive years.
INTEGER NOBS, NQPROP
PARAMETER (NOBS=30, NQPROP=5)
C
INTEGER I, NMISS, NOUT
REAL QPROP(NQPROP), X(NOBS), XEMP(NQPROP), XHI(NQPROP),
& XLO(NQPROP)
EXTERNAL EQTIL, UMACH
C
DATA X/0.77, 1.74, 0.81, 1.20, 1.95, 1.20, 0.47, 1.43, 3.37,
& 2.20, 3.00, 3.09, 1.51, 2.10, 0.52, 1.62, 1.31, 0.32, 0.59,
& 0.81, 2.81, 1.87, 1.18, 1.35, 4.75, 2.48, 0.96, 1.89, 0.90,
& 2.05/
DATA QPROP/0.01, 0.50, 0.90, 0.95, 0.99/
IMSL STAT/LIBRARY Chapter 1: Basic Statistics 37
C
CALL UMACH (2, NOUT)
CALL EQTIL (NOBS, X, NQPROP, QPROP, XEMP, XLO, XHI, NMISS)
WRITE (NOUT,99997)
99997 FORMAT ( Smaller Empirical Larger, /,
& Quantile Datum Quantile Datum)
DO 10 I=1, NQPROP
WRITE (NOUT,99998) QPROP(I), XLO(I), XEMP(I), XHI(I)
10 CONTINUE
99998 FORMAT (4X, F4.2, 8X, F4.2, 8X, F4.2, 8X, F4.2)
WRITE (NOUT,99999) NMISS
99999 FORMAT (/, There are , I2, missing values.)
END
Output
Smaller Empirical Larger
Quantile Datum Quantile Datum
0.01 0.32 0.32 0.32
0.50 1.43 1.47 1.51
0.90 3.00 3.08 3.09
0.95 3.37 3.99 4.75
0.99 4.75 4.75 4.75
There are 0 missing values.

Mon, 21 Jul 2008 16:20:33 GMT  Replacing IMSL call ?
I see no one else has responded to this, so I will.  It may be
difficult to find a replacement for this routine because the
problem is so simple.  Therefore, no one has bothered to publish
a routine for finding percentiles.

You find percentiles by sorting the data and then counting up
(or down) the desired fraction of the ordered values. Some
linear interpolation may be required.

So in summary you can write your own replacement routine with
the aid of a sort routine (you will find plenty on the Web) and
a basic statistics book (for defining percentiles where there
are ties or when interpolation is needed).  You might not even
need the stats book, just type "define percentile" into Google.

I hope that helps.

...Mike Prager

Quote:

>Hi - I am supporting a Win32 Console Application written with Microsoft
>Fortran Powerstation.  It makes 1 IMSL call to "eqtil".

>It reads numbers from a file, crunches them furiously, then writes out
>several files as an analysis of the data read, including one file with
>any error messages for the user.

>This compiler was some-what buggy, and sooo much is available now that
>wasn't when this program was written.

>I am wondering if anyone knows if there is free-ware fortran code out
>there now that replaces eqtil ?  (see manual on eqtil below)

>Thanks - Gary

>According to the manual:

>EQTIL/DEQTIL (Single/Double precision)
>Compute empirical quantiles.
>Usage
>CALL EQTIL (NOBS, X, NQPROP, QPROP, Q, XLO, XHI, NMISS)
>Arguments
>NOBS Number of observations. (Input)
>NOBS must be greater than or equal to one.
>X Vector of length NOBS containing the data. (Input)
>NQPROP Number of quantiles. (Input)
>NQPROP must be greater than or equal to one.
>QPROP Vector of length NQPROP containing the quantile proportions.
>(Input)
>The elements of QPROP must lie in the interval (0, 1).
>Q Vector of length NQPROP containing the empirical quantiles. (Output)
>Q(i) corresponds to the empirical quantile at proportion QPROP(i). The
>quantiles
>are determined by linear interpolation between adjacent ordered sample
>values.
>XLO Vector of length NQPROP containing the largest element of X less
>than or
>equal to the desired quantile. (Output)
>XHI Vector of length NQPROP containing the smallest element of X greater
>than or equal to the desired quantile. (Output)
>NMISS Number of missing values. (Output)
>1. Automatic workspace is allocated only if X is not sorted on input. The
>amount allocated is
>EQTIL NOBS units, or
>DEQTIL 2 * NOBS units.
>Workspace may be explicitly provided, if desired, by use of
>E2TIL/DE2TIL. The reference is
>CALL E2TIL (NOBS, X, NQPROP, QPROP, Q, XLO, XHI,
>NMISS, WK)
>36 Chapter 1: Basic Statistics IMSL STAT/LIBRARY
>WK Workspace of length NOBS containing the sorted data. (Output)
>If X is sorted in ascending order with all missing values at the end of X,
>then X and WK may share the same storage location.
>2. Informational error
>Type Code
>3 1 All of the observations are missing values. The
>elements of Q, XLO, and XHI have been set to NaN
>(not a number).
>3. Missing values (NaN) are excluded from the analysis. Empirical
>quantiles are based on the NOBS - NMISS nonmissing elements of X.
>Algorithm
>The routine EQTIL determines the empirical quantiles, as indicated in
>the vector
>QPROP, from the data in X. The routine EQTIL first checks to see if X is
>sorted; if
>X is not sorted, the routine does either a complete or partial sort,
>depending on
>how many order statistics are required to compute the quantiles requested.
>The routine EQTIL returns the empirical quantiles and, for each
>quantile, the two
>order statistics from the sample that are at least as large and at least
>as small as
>the quantile. For a sample of size n, the quantile corresponding to the
>proportion
>p is defined as
>Q(p) = (1 - f )x(M) + f x(M

+
1)

- Show quoted text -

Quote:
>where j = ?p(n + 1)?, f = p(n + 1) - j, and x(M) is the j-th order
>statistic, if
>1 j < n; otherwise, the empirical quantile is the smallest or largest
>order statistic.
>Example
>In this example, five empirical quantiles from a sample of size 30 are
>obtained.
>Notice that the 0.5 quantile corresponds to the sample median. The data
>are from
>Hinkley (1977) and Velleman and Hoaglin (1981). They are the
>measurements (in
>inches) of precipitation in Minneapolis/St. Paul during the month of
>March for 30
>consecutive years.
>INTEGER NOBS, NQPROP
>PARAMETER (NOBS=30, NQPROP=5)
>C
>INTEGER I, NMISS, NOUT
>REAL QPROP(NQPROP), X(NOBS), XEMP(NQPROP), XHI(NQPROP),
>& XLO(NQPROP)
>EXTERNAL EQTIL, UMACH
>C
>DATA X/0.77, 1.74, 0.81, 1.20, 1.95, 1.20, 0.47, 1.43, 3.37,
>& 2.20, 3.00, 3.09, 1.51, 2.10, 0.52, 1.62, 1.31, 0.32, 0.59,
>& 0.81, 2.81, 1.87, 1.18, 1.35, 4.75, 2.48, 0.96, 1.89, 0.90,
>& 2.05/
>DATA QPROP/0.01, 0.50, 0.90, 0.95, 0.99/
>IMSL STAT/LIBRARY Chapter 1: Basic Statistics 37
>C
>CALL UMACH (2, NOUT)
>CALL EQTIL (NOBS, X, NQPROP, QPROP, XEMP, XLO, XHI, NMISS)
>WRITE (NOUT,99997)
>99997 FORMAT ( Smaller Empirical Larger, /,
>& Quantile Datum Quantile Datum)
>DO 10 I=1, NQPROP
>WRITE (NOUT,99998) QPROP(I), XLO(I), XEMP(I), XHI(I)
>10 CONTINUE
>99998 FORMAT (4X, F4.2, 8X, F4.2, 8X, F4.2, 8X, F4.2)
>WRITE (NOUT,99999) NMISS
>99999 FORMAT (/, There are , I2, missing values.)
>END
>Output
>Smaller Empirical Larger
>Quantile Datum Quantile Datum
>0.01 0.32 0.32 0.32
>0.50 1.43 1.47 1.51
>0.90 3.00 3.08 3.09
>0.95 3.37 3.99 4.75
>0.99 4.75 4.75 4.75
>There are 0 missing values.

--
Mike Prager, NOAA, Beaufort, NC
* Opinions expressed are personal and not represented otherwise.
* Any use of tradenames does not constitute a NOAA endor{*filter*}t.

Tue, 22 Jul 2008 02:16:58 GMT  Replacing IMSL call ?

Quote:
> You find percentiles by sorting the data and then counting up
> (or down) the desired fraction of the ordered values. Some
> linear interpolation may be required.

Surprisingly, you can find percentiles in O(N) time, whereas
it takes O(N*log(N)) time to sort an array of N elements.
This may not be significant for a given application, however.

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end

Tue, 22 Jul 2008 03:18:29 GMT  Replacing IMSL call ?

Quote:
>Surprisingly, you can find percentiles in O(N) time, whereas
>it takes O(N*log(N)) time to sort an array of N elements.
>This may not be significant for a given application, however.

Yes, that is interesting.  Can you supply an algorithm or
reference to one?

--
Mike Prager, NOAA, Beaufort, NC
* Opinions expressed are personal and not represented otherwise.
* Any use of tradenames does not constitute a NOAA endor{*filter*}t.

Tue, 22 Jul 2008 06:09:36 GMT  Replacing IMSL call ?

Quote:

> >Surprisingly, you can find percentiles in O(N) time, whereas
> >it takes O(N*log(N)) time to sort an array of N elements.
> >This may not be significant for a given application, however.
> Yes, that is interesting.  Can you supply an algorithm or
> reference to one?

My recollection could be wrong on this one.  I thought I saw
the algorithm in a CS algorithms book; look through a couple
of likely sources -- I couldn't google up anything useful
before my next appointment.

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end

Tue, 22 Jul 2008 09:01:14 GMT  Replacing IMSL call ?
On 2006-02-02 18:09:36 -0400, Michael Prager

Quote:

>> Surprisingly, you can find percentiles in O(N) time, whereas
>> it takes O(N*log(N)) time to sort an array of N elements.
>> This may not be significant for a given application, however.

> Yes, that is interesting.  Can you supply an algorithm or
> reference to one?

The gotchas are:

1: fine as long as you want fewer than log(n) of them
2: the coefficient is larger than for sorting so fewer than log(n)

so think finding quartiles of a million values.

The original observation was that it is possible to obtain
the median in linear time. The basis was to think of the
data as a two way layout. Find the medians of the columns
and then the median of the column medians. Toss the elements
that are greater than the column median for those columns
with larger medians or less than the column median for those
columns with smaller medians. Repeat. Enough stuff gets tossed
to make it a linear algorithm but the coefficient is noticable.
Many many optimizations are possible. Generalizations for other
than the median.

There are special versions of quicksort in the Applied Statisitcs
algorithms section circa 1970. This is differeent algorithm than
the linear median algorithm but it works nicely with the usual
caveats on quicksort. Good on average but occasional high cost.

Tue, 22 Jul 2008 21:55:38 GMT  Replacing IMSL call ?

Quote:
> My recollection could be wrong on this one.  I thought I saw
> the algorithm in a CS algorithms book; look through a couple
> of likely sources -- I couldn't google up anything useful
> before my next appointment.

OK, this is as far as I have been able to get by my next
appointment:

! File: order_stat.f90
module mykinds
implicit none
integer, parameter :: sp = selected_real_kind(6,30)
integer, parameter :: dp = selected_real_kind(15,300)
end module mykinds

module order_stat
use mykinds, only: wp => sp
implicit none
contains
recursive subroutine select(array, N, k)
integer, intent(in) :: N
real(wp) array(N)
integer, intent(in) :: k
real(wp) med_loc((N+4)/5)
integer med_last
integer i
integer N_last
real(wp) mark
integer mark_count
integer lo_insert, hi_insert
real(wp) lo_hold(10), hi_hold(10)
integer lo_temp, hi_temp
integer i_hi, i_lo
integer j
real(wp), allocatable :: merge_temp(:)

! Handle wierd cases
if(k <= 0) then
write(*,'(a,i0,a)') "Invalid k = ", k, &
" detected in subroutine select"
stop
else if(k == 1) then
array((/1,minloc(array)/)) = array((/minloc(array),1/))
return
else if(k < N) then
! Normal case follows IF block
else if(k == N) then
array((/N,maxloc(array)/)) = array((/maxloc(array),N/))
return
else
write(*,'(a,i0,a)') "Invalid k = ", k, &
" detected in subroutine select"
stop
end if

! Normal case
! Handle small arrays via sorting
if(N <= 20) then
allocate(merge_temp((N+1)/2))
call merge_sort(array, N, merge_temp)
return
end if
! Find medians of subsets of 5
do i = 1, N/5
call median_5(array(5*i-4:5*i))
med_loc(i) = array(5*i-2)
end do
N_last = mod(N,5)
med_last = N+1-(N_last+1)/2
if(N_last > 0) then
call sort_short(array(N-N_last+1:N),N_last)
med_loc(i) = array(med_last)
end if
! Find median of medians
call select(med_loc,size(med_loc),(size(med_loc)+1)/2)
! Partition array around mark = median of medians
mark = med_loc((size(med_loc)+1)/2)
! Initialize partition code variables
mark_count = 0
lo_insert = 1
hi_insert = N
lo_temp = 1
hi_temp = 1
! Handle possible deficient subset
if(N_last > 0) then
if(mark <= array(med_last)) then
if(array(med_last) <= mark) then
mark_count = mark_count+1
hi_insert = med_last
else
hi_insert = med_last-1
end if
do j = med_last-1, N-N_last+1, -1
if(array(j) <= mark) then
lo_hold(lo_temp) = array(j)
lo_temp = lo_temp+1
else
array(hi_insert) = array(j)
hi_insert = hi_insert-1
end if
end do
else
do j = N, med_last+1, -1
if(array(j) <= mark) then
lo_hold(lo_temp) = array(j)
lo_temp = lo_temp+1
else
array(hi_insert) = array(j)
hi_insert = hi_insert-1
end if
end do
do j = med_last, N-N_last+1, -1
lo_hold(lo_temp) = array(j)
lo_temp = lo_temp+1
end do
end if
end if
i_hi = i-1
i_lo = 1
! Loop for general subsets
outer_partition: do
inner_partition_lo: do
if(i_lo > i_hi) exit outer_partition
if(array(5*i_lo-2) <= mark) then
if(mark <= array(5*i_lo-2)) then
mark_count = mark_count+1
do j = 5*i_lo-4, 5*i_lo-3
array(lo_insert) = array(j)
lo_insert = lo_insert+1
end do
else
do j = 5*i_lo-4, 5*i_lo-2
array(lo_insert) = array(j)
lo_insert = lo_insert+1
end do
end if
do j = 5*i_lo-1, 5*i_lo
if(array(j) <= mark) then
array(lo_insert) = array(j)
lo_insert = lo_insert+1
else
hi_hold(hi_temp) = array(j)
hi_temp = hi_temp+1
end if
end do
else
do j = 5*i_lo-4, 5*i_lo-3
if(array(j) <= mark) then
array(lo_insert) = array(j)
lo_insert = lo_insert+1
else
hi_hold(hi_temp) = array(j)
hi_temp = hi_temp+1
end if
end do
do j = 5*i_lo-2, 5*i_lo
hi_hold(hi_temp) = array(j)
hi_temp = hi_temp+1
end do
end if
i_lo = i_lo+1
if(hi_temp > 5) exit inner_partition_lo
end do inner_partition_lo
inner_partition_hi: do
if(i_lo > i_hi) exit outer_partition
if(mark <= array(5*i_hi-2)) then
if(array(5*i_hi-2) <= mark) then
mark_count = mark_count+1
do j = 5*i_hi, 5*i_hi-1, -1
array(hi_insert) = array(j)
hi_insert = hi_insert-1
end do
else
do j = 5*i_hi, 5*i_hi-2, -1
array(hi_insert) = array(j)
hi_insert = hi_insert-1
end do
end if
do j = 5*i_hi-3, 5*i_hi-4, -1
if(mark <= array(j)) then
array(hi_insert) = array(j)
hi_insert = hi_insert-1
else
lo_hold(lo_temp) = array(j)
lo_temp = lo_temp+1
end if
end do
else
do j = 5*i_hi, 5*i_hi-1, -1
if(mark <= array(j)) then
array(hi_insert) = array(j)
hi_insert = hi_insert-1
else
lo_hold(lo_temp) = array(j)
lo_temp = lo_temp+1
end if
end do
do j = 5*i_hi-2, 5*i_hi-4, -1
lo_hold(lo_temp) = array(j)
lo_temp = lo_temp+1
end do
end if
i_hi = i_hi-1
if(lo_temp > 5) exit inner_partition_hi
end do inner_partition_hi
do j = 1, 5
lo_temp = lo_temp-1
array(lo_insert) = lo_hold(lo_temp)
lo_insert = lo_insert+1
hi_temp = hi_temp-1
array(hi_insert) = hi_hold(hi_temp)
hi_insert = hi_insert-1
end do
end do outer_partition
do lo_temp = lo_temp-1, 1, -1
array(lo_insert) = lo_hold(lo_temp)
lo_insert = lo_insert+1
end do
do hi_temp = hi_temp-1, 1, -1
array(hi_insert) = hi_hold(hi_temp)
hi_insert = hi_insert-1
end do
array(lo_insert:hi_insert) = mark
if(k < lo_insert) then
call select(array(1:lo_insert-1), lo_insert-1, k)
else if(k > hi_insert) then
call select(array(hi_insert+1:N), N-hi_insert, k-hi_insert)
end if

contains
subroutine median_5(A5)
real(wp) A5(5)

if(.NOT. A5(1) <= A5(2)) then
A5(1:2) = A5((/2,1/))
end if
if(.NOT. A5(2) <= A5(3)) then
if(.NOT. A5(1) <= A5(3)) then
A5(1:3) = A5((/3,1,2/))
else
A5(2:3) = A5((/3,2/))
end if
end if
if(.NOT. A5(4) <= A5(5)) then
A5(4:5) = A5((/5,4/))
end if
if(A5(4) <= A5(2)) then
if(A5(2) <= A5(5)) then
A5(2:4) = A5((/4,2,3/))
else
A5(2:5) = A5((/4,5,2,3/))
end if
else
if(.NOT. A5(3) <= A5(4)) then
A5(3:4) = A5((/4,3/))
end if
end if
end subroutine median_5

subroutine sort_short(A4,N)
integer, intent(in) :: N
real(wp) :: A4(N)
integer i, j

do j = N, 2, -1
do i = 1, j-1
if(.NOT. A4(i) <= A4(i+1)) then
A4(i:i+1) = A4((/i+1,i/))
end if
end do
end do
end subroutine sort_short
end subroutine select

recursive subroutine merge_sort(array, N, temp)
integer, intent(in) :: N
real(wp) array(N)
real(wp) temp(*)
integer lo_remove, hi_remove
integer j

! Handle small arrays
if(N < 0) then
write(*,'(a,i0,a)') "Invalid array size = ", N, &
"detected in subroutine merge_sort"
...