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)
Comments
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)
The additional argument is
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)
>Comments
>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)
>The additional argument is
>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
Address spam-trapped; remove color to reply.
* 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
Address spam-trapped; remove color to reply.
* 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"
...

read more »



Thu, 24 Jul 2008 09:35:15 GMT  
 
 [ 7 post ] 

 Relevant Pages 

1. replacing imsl routines

2. Calling IMSL Subroutines in Digital Visual Fortran

3. VB or VF to call imsl fortran routines

4. How to call IMSL subroutine in Powerstation 4.0

5. calling imsl routines from c programs

6. Old IMSL calls

7. REPLACE not replacing under Win 95??

8. How to replace or create a file using the open/create/replace.vi

9. Issues with Non-COBOL characters and REPLACE / REPLACING

10. string.replace() can't replace newline characters???

11. IMSL & NAG on Sun Workstations

12. [Q: DLL - USE IMSL : EXTERNAL subroutine]

 

 
Powered by phpBB® Forum Software