JiscMail Logo
Email discussion lists for the UK Education and Research communities

Help for COMP-FORTRAN-90 Archives


COMP-FORTRAN-90 Archives

COMP-FORTRAN-90 Archives


COMP-FORTRAN-90@JISCMAIL.AC.UK


View:

Message:

[

First

|

Previous

|

Next

|

Last

]

By Topic:

[

First

|

Previous

|

Next

|

Last

]

By Author:

[

First

|

Previous

|

Next

|

Last

]

Font:

Proportional Font

LISTSERV Archives

LISTSERV Archives

COMP-FORTRAN-90 Home

COMP-FORTRAN-90 Home

COMP-FORTRAN-90  May 2018

COMP-FORTRAN-90 May 2018

Options

Subscribe or Unsubscribe

Subscribe or Unsubscribe

Log In

Log In

Get Password

Get Password

Subject:

Re: automatic sample size determination, while construct, gfortran compiler option (end of file)

From:

Mudelsee M <[log in to unmask]>

Reply-To:

Fortran 90 List <[log in to unmask]>

Date:

Tue, 15 May 2018 08:22:32 +0200

Content-Type:

multipart/mixed

Parts/Attachments:

Parts/Attachments

text/plain (183 lines) , kernel.f90 (1 lines)

Dear Bill,

thanks for the reply. I agree that malformed files (i.e., others than 
file "a" I sent with my inquiry message) should be avoided, but that is 
life if you give courses and participants bring own data.

However: I tried your suggestion by adding/replacing text with

use,intrinsic :: iso_fortran_env
if (iocheck == IOSTAT_END) exit

but the results for file "d" are still negative ("n = 820") on my 
machines here:

o gfortran that is included in gcc-4.6.0 (32 bit, Windows)
o an old (2007) IBM workstation under Windows XPP 64 bit
o an old (2010) PC under Windows 7 Prof 64 bit

(I doubt that is has something to do with 32/64 bit.)

For the sake of completeness, I attach one code here (there are several 
codes I use for teaching, which all have the same data size 
determination). There your suggestion has been built into subroutine 
chsett5.

Cheers

Manfred (puzzled)

Am 14.05.2018 um 21:19 schrieb Bill Long:
> Note that the wc system utility also thinks that d.txt has 820 lines:
> 
> % wc d.txt
>       820    1642   11402 d.txt
> 
> This is because the last “line” of the file is malformed.  By checking for iostat .ne. 0 you are tripping on that error, rather than the EOF indicator.  If the code is modified to
> 
> use,intrinsic :: iso_fortran_env
>       integer :: i       ! counter
>       integer :: iocheck ! check
>       integer :: n       ! sample size
>       real :: t          ! time
>       real :: x          ! variable
>       character(len=*),parameter :: file1 = "d.txt"
> 
> 
> !
> ! Determine data size automatically
> ! =================================
> !
>       open (unit=1, file=file1, status='old')
>       i=1
>       do
>          read (unit=1, fmt=*, iostat=iocheck) t, x
>          if (iocheck == IOSTAT_END) exit
>          i=i+1
>       end do
>       close (unit=1, status='keep')
>       n=i-1
> 
>       print *, n
> end
> 
> 
> then it prints 821:
> 
> % gf d.f90
> % ./a.out
>           821
> 
> [gf is my shorthand for gfortran]
> 
> Cheers,
> Bill
> 
> 
> 
>> On May 14, 2018, at 1:59 AM, Mudelsee M <[log in to unmask]> wrote:
>>
>> Dear colleagues,
>>
>> I am looking for help in the interpretation and correction of an issue (which seemingly has something to do with the end-of-file situation) that arises when I compile a code under gfortran.
>>
>> Below code extract shows the while construct used to determine the sample size of a time series. The idea is that after iocheck is .ne. 0 (end reached), the counting is stopped. Automatic sample size determination is convenient in applied time series analysis, so I wish to keep this.
>>
>> The issue now is that a gfortran-compiled (just the -O2 option) executable may (depending on how the input data file is, see the file description below) yield an n value that is less by one than the true value. From the 5 attached test files (a, b, c, d, e), this happens for file1 = d (which produces n = 820 instead of 821); the other files correctly yield n = 821.
>>
>> This issue does not occur when the code is compiled under ifort (I do not remember by heart the compiling options).
>>
>> The 5 files (attached) contain the same data but diverge with regard to the end portion: carriage-return (CR), line-feed (LF), space, decimal point. The ominous file "d" has:
>> o no decimal point (i.e., an integer is read for x)
>> o no space after the last value
>> o no extra line
>> The attached scan shows the end portions of the files within Notepad++ editor (the "[NOTE SPACE]" is inserted by me for clarification).
>>
>> I have the intuition that an added gfortran compiling option may help (i.e., makes also "d" yield correctly n = 821), but I am not an expert on gfortran, and I thought there may be guys out there who quickly spot how to correct this.
>>
>> I hope that the amount of information is sufficient and that it is possible to understand the issue without being given the full source code.
>>
>> Many thanks!
>>
>> Manfred
>>
>>
>> === code extract ===
>>
>>
>>       integer :: i       ! counter
>>       integer :: iocheck ! check
>>       integer :: n       ! sample size
>>       real :: t          ! time
>>       real :: x          ! variable
>>       character(len=*) :: file1
>>
>> !
>> ! Determine data size automatically
>> ! =================================
>> !
>>       open (unit=1, file=file1, status='old')
>>       i=1
>>       do while (.true.)
>>          read (unit=1, fmt=*, iostat=iocheck) t, x
>>          if (iocheck .ne. 0) exit
>>          i=i+1
>>       end do
>>       close (unit=1, status='keep')
>>       n=i-1
>>
>> ===================
>>
>>
>>
>> -- 
>> Dr. Manfred Mudelsee
>>
>> Chief Executive Officer
>> Climate Risk Analysis
>> Kreuzstrasse 27
>> Heckenbeck
>> 37581 Bad Gandersheim
>> Germany
>>
>> Telephone: +49 5563 9998140
>> Email: [log in to unmask]
>> URL: http://www.climate-risk-analysis.com
>> Skype: mudelsee1
>> LinkedIn: https://de.linkedin.com/in/mudelsee
>> Twitter: @MMudelsee
>>
>> Climate Time Series and Risk Analyses
>> Book: http://www.manfredmudelsee.com/book/
>> Courses: http://www.climate-risk-analysis.com/courses/
>> <a.txt><b.txt><c.txt><d.txt><e.txt><Scannen0001.pdf>
> 
> Bill Long                                                                       [log in to unmask]
> Principal Engineer, Fortran Technical Support &   voice:  651-605-9024
> Bioinformatics Software Development                      fax:  651-605-9143
> Cray Inc./ 2131 Lindau Lane/  Suite 1000/  Bloomington, MN  55425
> 
> 

-- 
Dr. Manfred Mudelsee

Chief Executive Officer
Climate Risk Analysis
Kreuzstrasse 27
Heckenbeck
37581 Bad Gandersheim
Germany

Telephone: +49 5563 9998140
Email: [log in to unmask]
URL: http://www.climate-risk-analysis.com
Skype: mudelsee1
LinkedIn: https://de.linkedin.com/in/mudelsee
Twitter: @MMudelsee

Climate Time Series and Risk Analyses
Book: http://www.manfredmudelsee.com/book/
Courses: http://www.climate-risk-analysis.com/courses/



! ! kernel 3.2 ! ====== ! ! o Gasser M\"uller kernel smoothing ! o polynomial kernel, modified at time interval boundaries ! o 0th, 1st or 2nd derivative estimation ! o bandwidth: global (predefined) or local (optimized) ! o persistence time: either predefine or write residuals to /tauest/ and perform externally estimation ! o MBB bootstrap se band ! o no timescale simulations ! ! o envisaged but not done yet: allow predefined local bandwidth hy ! !======================================================================================================================== ! ! Delevopment directories (current and previous) ! ======================= ! ! E:/kernel/ (current) ! ! F:/CRA/Projects/Cardiff-Bark/Training/kernel/ ! F:/CRA(Products/Course/software/kernel/kernel/ ! F:/CRA/Projects/Daphne-II/SI-Paper/kernel/ ! F:/IPCC/AR5/Ch5/kernel/ ! F:/Cenozoic/stack/kernel/ ! !======================================================================================================================== ! ! Language ! ======== ! ! FORTRAN 77, Fortran 90 ! !======================================================================================================================== ! ! Machine and system ! ================== ! ! (1) ! IBM Intelli Workstation Z Pro ! 4 x Intel Xeon 5160 (3 GHz) CPU ! 8 GB RAM ! Windows XP Professional x64 Edition Service Pack 2 ! ! (2) ! Samsung X360 ! 2 x Intel x86 (1.2 GHz) CPU ! 2 GB RAM ! Windows XP Professional Service Pack 3 ! !======================================================================================================================== ! ! Compiler ! ======== ! ! gfortran (no options, -O2) ! ! Note: also tested Absoft 6.2, Console Application, no options ! Stack: 0 x 20000000 reserve, 0 x 200000 commit ! Heap: 0 x 2500000 reserve, 0 x 2500 commit ! !======================================================================================================================== ! ! Author (Fortran 90 part) ! ====== ! ! Manfred Mudelsee ! Climate Risk Analysis ! Schneiderberg 26 ! 30167 Hannover ! Germany ! Email: [log in to unmask] ! URL: http://www.climate-risk-analysis.com ! !======================================================================================================================== ! ! Change log ! ========== ! ! Version Date Comments ! ! ! 3.2 January output time points selection augmented, ! residuals defined, ! (subroutines chsettpar1, output) ! ! 3.1 July 2014 output bandwidth also for global smoothing, ! output +- 2 se in addition to +- 1 se ! ! 3.0 October 2013 data file selectable ! bandwidth: global (predefined) or local (optimized) ! persistence time: either predefine or write residuals to ! /tauest/ and perform externally estimation ! 0th, 1st or 2nd derivative estimation ! variables til and tiu 'parameter' attribute removed ! ! 2.4 July 2013 module kernelparameters, ! variable sig 'parameter' attribute removed ! latest version (23 July 2013) of kernel.f90 is at ! F:/CRA/Products/course/soft/kernel/kernel/ ! ! course version (only data file BU4-s, ! see F:/CRA/Projects/Daphne-II/SI-Paper/kernel/) ! ! 2.3 April 2012 Experiment 2 added (combining StalAge and iscam) ! ! 2.2 April 2012 hyglobal set to 0.25 ka (from formerly 0.3 ka) ! because of former problems with boundary kernel ! for lower bound of time intervals for ! Bu11s and Bu4-s (subroutine GLKERN) ! ! 'present = AD 1950' to be set ! (new subroutine setpresent; ! module parameters: new tlow, tupp) ! ! 2.1 March 2012 Timescale simulations done (yes/no) ! on identical bootstrap versions of x ! ! MBB replaces ARB ! ! 2.0 January 2012 Adaptation to Daphne-II ! o notation & names ! o read tboot from data file ! o addition of GLKERN subroutine ! o tau may be set ! ! 1.0 December 2011 Original version from /Cenozoic/ ! !======================================================================================================================== ! ! Numerical Recipes modules ! ========================= ! include 'nrtype.f90' ! gfortran: uncomment this line include 'nrutil.f90' ! gfortran: uncomment this line include 'nr.f90' ! gfortran: uncomment this line ! !======================================================================================================================== ! module parameters use nrtype implicit none ! Constants, parameters and adjustable parameters (not for kernel estimation). save ! integer, parameter :: nb=400 ! number of bootstrap resamplings (default: nb = 400) real(sp) :: du=-999.0_sp ! spacing of time grid u real(sp) :: tlow=-999.0_sp ! lower bound of time interval real(sp) :: tupp=-999.0_sp ! upper bound of time interval integer :: nu=-999 ! dimension of time grid u real(sp), allocatable, dimension(:) :: u ! time grid end module parameters ! !======================================================================================================================== ! module kernelparameters use nrtype implicit none ! Constants, parameters and adjustable parameters for kernel estimation. save integer, parameter :: ihom=1 ! heteroscedasticity (ihom = 1) / homoscedasticity (ihom = 0) ! ! (default: ihom = 1) integer, parameter :: irnd=1 ! spacing type (1: non-random; 0: random) ! (default: irnd = 1) integer, parameter :: ismo=0 ! bandwidth selection (1: predefined; 0: not predefined, to be estimated) ! (default: ismo = 0) integer, parameter :: ismob=1 ! bandwidth selection for bootstrap ! (1: overtake from estimation; 0: re-estimate) ! (default: ismob = 1) integer :: ismog=-999 ! bandwidth selection for global smoothing (1: predefined; 0: not predefined) ! (default: ismog = 1) integer, parameter :: ismofix=1 ! bandwidth selection for global/local smoothing ! (1: predefined; 0: not predefined) ! (default: ismofix = 1, i.e., overtake previously determined bandwidth)) integer :: kord=-999 ! order of polynomial kernel ! ! adjustable ! ! (default: kord = nue + 2) integer :: kordx=2 ! order of polynomial kernel for residuals calculation ! (only choice: kordx=nuex + 2 = 2) integer, parameter :: m1=2000 ! integral approximation, higher accuracy <-> higher m1 ! ! (default: m1 = 2000, note: subroutine lokern gives default m1 = 400) integer :: nue=-999 ! order of derivative ! ! adjustable integer, parameter :: nuex=0 ! order of derivative for residuals calculation ! (only choice: nuex = 0) real(sp) :: sig=-1.0_sp ! residual variance ! ! (default: sig = -1.0) real(sp) :: til=1.0_sp ! integral approximation, lower bound ! ! (default: til = 1.0) real(sp) :: tiu=0.0_sp ! integral approximation, upper bound ! ! (default: tiu = 0.0) end module kernelparameters ! !======================================================================================================================== ! module data1 use nrtype implicit none ! Data. save real(sp), allocatable, dimension(:) :: e ! unweighted residual real(sp), allocatable, dimension(:) :: t ! time real(sp), allocatable, dimension(:) :: x ! delta 18-O real(sp), allocatable, dimension(:) :: xboot ! bootstrap version of x real(sp), allocatable, dimension(:) :: xtrend ! trend, original time points end module data1 ! !======================================================================================================================== ! module kernelbandwidth use nrtype use parameters, only: nu implicit none ! Bandwidth for kernel estimation. save real(sp), allocatable, dimension(:) :: hx ! on t(1:n), x(1:n) for residual determination real(sp), allocatable, dimension(:) :: hy ! on u(1:nu), y(1:nu) for smoothing local real(sp) :: hyglobal=-999.0_sp ! on u(1:nu), y(1:nu) for smoothing global end module kernelbandwidth ! !======================================================================================================================== ! module kernelwork use nrtype use parameters, only: nu implicit none ! Work arrays for kernel estimation. save real(sp), allocatable, dimension(:) :: s real(sp), allocatable, dimension(:,:) :: w1 real(sp), allocatable, dimension(:) :: wmx real(sp), allocatable, dimension(:) :: wmy real(sp), allocatable, dimension(:,:) :: wn end module kernelwork ! !======================================================================================================================== ! module persist use nrtype implicit none ! Constants for tau estimation. save character(len=13), parameter :: taufile='tauest/te.dat' character(len=1) :: tauflag='n' ! 'y': t(i) e(i) data written to taufile, no bootstrap ! ! 'n': bootstrap integer :: l ! MBB block length real(sp) :: tau=-999.0_sp ! persistence time end module persist ! !======================================================================================================================== ! module result1 use nrtype use parameters, only: nb,nu implicit none ! Result. real(sp), allocatable, dimension(:) :: y ! trend on time grid u real(sp), allocatable, dimension(:) :: seytbootno ! se(trend) on time grid u, t* no real(sp), allocatable, dimension(:) :: ybootno ! bootstrap version of y (subroutine smooth), t* no save end module result1 ! !======================================================================================================================== ! module setting use nrtype implicit none ! Setting. save character(len=1) :: smtype='g' ! 'g': smoothing: global ! ! 'l': local character (len=40) :: infile ! input data file integer :: n=-999 ! sample size (within [tlow; tupp]) character (len=80) :: outfile ! output data file end module setting ! !======================================================================================================================== ! program kernel1 use persist, only: tauflag use setting, only: smtype implicit none ! ! 1. Initial stuff ! ============= ! call ranseed ! seeds random number generator call welcome call chsett1 ! changes setting: infile call chsett2 ! changes setting: smtype call chsettpar0 ! changes setting: order of derivative if (smtype .eq. 'g') call chsett3 ! changes setting: hyglobal if (tauflag .eq. 'n') call chsett4 ! changes setting: tau ! ! 2. Estimation and bootstrap resampling or residuals writing ! ======================================================== ! call chsett5 ! changes setting: n call allocate1 ! allocates and initializes (module data1) call allocate2 ! allocates and initializes (module kernelwork) call allocate3 ! allocates and initializes hx call read1 ! reads t, x call chsettpar1 ! sets tlow, tupp, nu, allocates and sets u call allocate4 ! allocates hy, y, wmy, seytbootno, ybootno call chsett6 ! changes setting: l call smooth ! global smoothing (determines y, seytbootno) call chsett7 ! changes setting: outfile call output ! outputs u, y, hy, seytbootno, y + seytbootno, y - seytbootno call deallocate4 ! deallocates hy, y, wmy, seytbootno, ybootno, u call deallocate3 ! deallocates hx call deallocate2 ! deallocates (module kernelwork) call deallocate1 ! deallocates (module data1) end program kernel1 ! !======================================================================================================================== ! ! Fortran 90 subroutines ! !======================================================================================================================== ! subroutine allocate1 use nrtype use data1 use parameters, only: nb use setting, only: n implicit none ! Allocates and initializes (module data1). real(sp), parameter :: zero=0.0_sp allocate(e(n)) allocate(t(n)) allocate(x(n)) allocate(xboot(n)) allocate(xtrend(n)) e(1:n)=zero t(1:n)=zero x(1:n)=zero xboot(1:n)=zero xtrend(1:n)=zero end subroutine allocate1 ! !======================================================================================================================== ! subroutine allocate2 use nrtype use kernelparameters, only: m1 use kernelwork use setting, only: n implicit none ! Allocates and initializes (module kernelwork). real(sp), parameter :: zero=0.0_sp allocate(s(0:n)) allocate(w1(1:m1,1:3)) allocate(wmx(1:n)) allocate(wn(0:n,1:5)) call inits w1(1:m1,1:3)=zero wmx(1:n)=zero wn(0:n,1:5)=zero end subroutine allocate2 ! !======================================================================================================================== ! subroutine allocate3 use nrtype use kernelbandwidth, only: hx use setting, only: n implicit none ! Allocates and initializes hx. allocate(hx(n)) hx=-999.0_sp end subroutine allocate3 ! !======================================================================================================================== ! subroutine allocate4 use nrtype use kernelbandwidth, only: hy use kernelwork, only: wmy use parameters, only: nu use result1, only: y,seytbootno,ybootno implicit none ! Allocates hy, y, wmy, seytbootno, ybootno. allocate(hy(nu)) hy=-999.0_sp allocate(y(nu)) y=-999.0_sp allocate(wmy(nu)) wmy=0.0_sp allocate(seytbootno(nu)) seytbootno=-999.0_sp allocate(ybootno(nu)) ybootno=-999.0_sp end subroutine allocate4 ! !======================================================================================================================== ! subroutine avevar_n(n,x,ave,var) use nrtype implicit none integer, intent(in) :: n real(sp), dimension(n), intent(in) :: x real(sp), intent(out) :: ave,var ! Calculates VAR with n instead of n-1. real(sp), dimension(n) :: s ave=sum(x(:))/n s(:)=x(:)-ave var=dot_product(s,s) var=var/n end subroutine avevar_n ! !======================================================================================================================== ! subroutine bootstrap1(n,x,l,x_resample) use nrtype use nr, only: ran implicit none integer, intent(in) :: n real(sp), dimension(n), intent(in) :: x integer, intent(in) :: l real(sp), dimension(n), intent(out) :: x_resample ! ! MBB: Moving block bootstrap resampling (K\"unsch 1989). ! Block length = l. ! Blocks are randomly selected and concatenated until n data ! are resampled (resize last block by removing last observations in it). ! ! Note: MBB is faster (factor 1.5 to 2) than SB. ! integer, dimension(n) :: indxx integer :: i=0 integer :: j=0 integer :: k=0 integer :: idum=1 integer :: n_block_start n_block_start=n-l+1 if (l == 1) then do i=1,n indxx(i)=int(n*ran(idum))+1 end do else k=1 outer: do i=1,n indxx(k)=int(n_block_start*ran(idum))+1 k=k+1 if (k > n) exit outer do j=1,l-1 indxx(k)=indxx(k-1)+1 k=k+1 if (k > n) exit outer end do end do outer end if x_resample(:)=x(indxx(:)) end subroutine bootstrap1 ! !======================================================================================================================== ! subroutine chsett1 use setting, only: infile implicit none ! Changes setting: infile. print '(a)', ' Data file (directory ''/data/'') name with extension): [INPUT]' read (unit=5, fmt='(a)') infile end subroutine chsett1 ! !======================================================================================================================== ! subroutine chsett2 use setting, only: smtype implicit none ! Changes setting: smtype. print '(a)', ' ' print '(a)', ' ' print '(a)', ' ' print '(a)', ' Select smoothing type: global [g]' print '(a)', ' local [l]' read (unit=5, fmt='(a1)') smtype end subroutine chsett2 ! !======================================================================================================================== ! subroutine chsett3 use nrtype use kernelbandwidth, only: hyglobal use kernelparameters, only: ismog implicit none ! Changes setting: hyglobal. real(sp), parameter :: zero=0.0_sp real(sp), parameter :: one=1.0_sp real(sp) :: dummy1 print '(a)', ' ' print '(a)', ' ' print '(a)', ' ' print '(a)', ' Global bandwidth: predefine [INPUT positive value]' print '(a)', ' optimize [INPUT arbitrary negative value]' read (unit=5, fmt=*) dummy1 if (dummy1 .le. zero) then ismog=0 else ismog=1 hyglobal=dummy1*one end if end subroutine chsett3 ! !======================================================================================================================== ! subroutine chsett4 use nrtype use persist, only: tau implicit none ! Changes setting: tau. real(sp), parameter :: one=1.0_sp print '(a)', ' ' print '(a)', ' ' print '(a)', ' ' print '(a)', ' Select persistence time tau: [INPUT]' read (unit=5, fmt=*) tau tau=tau*one end subroutine chsett4 ! !======================================================================================================================== ! subroutine chsett5 use, intrinsic :: iso_fortran_env ! MAY 2018 !use nrtype use setting, only: infile,n implicit none ! Changes setting: n. integer :: i integer :: iocheck integer :: openerror ! real(sp) :: dummy1 ! dummy ! real(sp) :: dummy2 ! dummy ! real :: dummy1 ! dummy real :: dummy2 ! dummy ! ! ! 1. Open file ! ========= ! open (unit=1, file='data/'//infile, status='old', & form='formatted', action='read', iostat=openerror) if (openerror /= 0 ) then print *,'input file: ',infile,' Error during data file opening---please check' return end if ! ! ! 2. Count lines ! =========== ! i=1 do read (unit=1, fmt=*, iostat=iocheck) dummy1,dummy2 ! if (iocheck .ne. 0) exit if (iocheck == IOSTAT_END) exit i=i+1 end do close (unit=1, status='keep') n=i-1 ! ! test MAY 2018 ! print '(a,a,a,i4)', & 'infile: ',infile,' n = ',n read (*,'()') end subroutine chsett5 ! !======================================================================================================================== ! subroutine chsett6 use nrtype use data1, only: t use persist, only: l,tau use setting, only: n implicit none ! ! Changes setting: l. ! References: Carlstein et al. 1986, Sherman et al. 1998, Mudelsee 2010. ! real(sp) :: a ! equivalent autocorrelation coefficient real(sp) :: dbar ! average spacing real(sp) :: dummy1 real(sp) :: dummy2 real(sp) :: dummy3 real(sp), parameter :: zero=0.0_sp real(sp), parameter :: onethird=0.3333333333333333333333333333333333333333333333333333333333333333333333333333333_sp real(sp), parameter :: half=0.5_sp real(sp), parameter :: twothird=0.6666666666666666666666666666666666666666666666666666666666666666666666666666666_sp real(sp), parameter :: one=1.0_sp real(sp), parameter :: two=2.0_sp real(sp), parameter :: six=6.0_sp if (tau .le. zero) then l=1 else dbar=(t(n)-t(1))/(n-one) a=exp(-dbar/tau) dummy1=six**half*a dummy2=one-a**two dummy3=(dummy1/dummy2)**twothird*n**onethird l=anint(dummy3) l=max(1,l) l=min(n-1,l) end if end subroutine chsett6 ! !======================================================================================================================== ! subroutine chsett7 use setting, only: infile,outfile implicit none ! Changes setting: outfile. character(len=40) :: dummyc1 dummyc1=adjustr(infile) outfile=trim(adjustl(dummyc1(1:36)))//'_result.txt' ! ! test ! print *, 'Outfile = ',outfile ! read (*,'()') end subroutine chsett7 ! !======================================================================================================================== ! subroutine chsettpar0 use kernelparameters, only: kord,nue use persist, only: tauflag implicit none ! Changes setting: order of derivative. nue=0 if (tauflag .eq. 'n') then print '(a)', ' ' print '(a)', ' ' print '(a)', ' ' print '(a)', ' Select order of derivative: 0th (trend) [0]' print '(a)', ' 1st (slope) [l]' print '(a)', ' 2nd (curvature) [2]' read (unit=5, fmt=*) nue end if kord=nue+2 end subroutine chsettpar0 ! !======================================================================================================================== ! subroutine chsettpar1 use nrtype use parameters, only: du,nu,tlow,tupp,u use persist, only: tauflag use data1, only: t use setting, only: n implicit none ! Sets tlow, tupp, nu, u. character(len=1) :: answer1 character(len=1000) :: filet integer :: i integer :: iocheck integer :: iu real(sp) :: dbar ! average spacing real(sp) :: dummy1 ! dummy real(sp), parameter :: one=1.0_sp tlow=t(1) tupp=t(n) dbar=(tupp-tlow)/(n-one) print '(a)', ' ' print '(a)', ' ' print '(a)', ' ' print '(a,1x,f18.6)', & ' Selected time series has start tlow =',tlow print '(a,1x,f18.6)', & ' end tupp =',tupp print '(a,1x,f18.6)', & ' average spacing dbar =',dbar print '(a)', ' ' if (tauflag .eq. 'y') then print '(a,a)', ' For persistence time estimation',',' print '(a)', ' the same time spacing is used as for kernel estimation' print '(a)', ' ' answer1='o' else print '(a)', & ' Select output time points: [tlow; tupp] and even spacing [e]' print '(a)', & ' [tlow; tupp] and original time [o]' print '(a)', & ' read from file [f]' read (unit=5, fmt='(a1)') answer1 end if if (answer1 .eq. 'e') then print '(a)', ' Select time spacing: [INPUT]' read (unit=5, fmt=*) du du=du*one nu=int((tupp-tlow)/du)+1 allocate(u(nu)) u=(/ (tlow+du*(iu-1), iu=1,nu) /) else if (answer1 .eq. 'o') then nu=n allocate(u(nu)) u=t else if (answer1 .eq. 'f') then print '(a)', ' Output time file: [INPUT]' read (unit=5, fmt='(a)') filet filet=trim(adjustl(filet)) open (unit=1, file=filet, status='old', form='formatted', action='read') i=1 do while (.true.) read (unit=1, fmt=*, iostat=iocheck) dummy1 if (iocheck .ne. 0) exit i=i+1 end do close (unit=1, status='keep') nu=i-1 allocate(u(nu)) open (unit=1, file=filet, status='old', form='formatted', action='read') do i=1,nu read (unit=1, fmt=*) u(i) u(i)=u(i)*one if (i .gt. 1 .and. u(i) .le. u(i-1)) & print '(a,i4,a,i4,a,a)', & 'Subroutine chsettpar1 --- nonincreasing t(',i,') of ',filet end do close (unit=1, status='keep') end if end subroutine chsettpar1 ! !======================================================================================================================== ! subroutine deallocate1 use data1 implicit none ! Deallocates (module data1). deallocate(e) deallocate(t) deallocate(x) deallocate(xboot) deallocate(xtrend) end subroutine deallocate1 ! !======================================================================================================================== ! subroutine deallocate2 use kernelwork implicit none ! Deallocates (module kernelwork). deallocate(s) deallocate(w1) deallocate(wmx) deallocate(wn) end subroutine deallocate2 ! !======================================================================================================================== ! subroutine deallocate3 use kernelbandwidth, only:hx implicit none ! Deallocates hx. deallocate(hx) end subroutine deallocate3 ! !======================================================================================================================== ! subroutine deallocate4 use kernelbandwidth, only:hy use kernelwork, only: wmy use parameters, only: u use result1, only: y,seytbootno,ybootno implicit none ! Deallocates hy, y, wmy, seytbootno, ybootno, u. deallocate(hy) deallocate(y) deallocate(wmy) deallocate(seytbootno) deallocate(ybootno) deallocate(u) end subroutine deallocate4 ! !======================================================================================================================== ! subroutine gasdev(harvest) use nrtype use nr, only : ran implicit none real(sp), intent(out) :: harvest ! Gaussian distribution (Numerical Recipes, modified: uses ran). real(sp) :: rsq,v1,v2 real(sp), save :: g logical, save :: gaus_stored=.false. integer :: idum=1 if (gaus_stored) then harvest=g gaus_stored=.false. else do ! call ran(v1) ! call ran(v2) v1=ran(idum) v2=ran(idum) v1=2.0_sp*v1-1.0_sp v2=2.0_sp*v2-1.0_sp rsq=v1**2+v2**2 if (rsq > 0.0 .and. rsq < 1.0) exit end do rsq=sqrt(-2.0_sp*log(rsq)/rsq) harvest=v1*rsq g=v2*rsq gaus_stored=.true. end if end subroutine gasdev ! !======================================================================================================================== ! subroutine inits use nrtype use kernelwork, only: s use setting, only: n implicit none ! Initializes s (module kernelwork). ! Note: subroutines glkern and lokern may change s; therefore, ! (1) s must not have 'parameter' attribute and ! (2) s has to be re-initialized using this subroutine. real(sp), parameter :: zero=0.0_sp real(sp), parameter :: one=1.0_sp s(0)=one s(1:n)=zero end subroutine inits ! !======================================================================================================================== ! subroutine initti use nrtype use kernelparameters, only: til,tiu implicit none ! Initializes til, tiu (module parameters). ! Note: subroutines glkern and lokern may change til or tiu; therefore, ! (1) til and tiu must not have 'parameter' attribute and ! (2) til and tiu have to be re-initialized using this subroutine. real(sp), parameter :: zero=0.0_sp real(sp), parameter :: one=1.0_sp til=one tiu=zero end subroutine initti ! !======================================================================================================================== ! subroutine output use nrtype use data1, only: t,x use kernelbandwidth, only: hy,hyglobal use kernelparameters, only: ismog,nue use parameters, only: nu,u use persist, only: l,tau use result1, only: y,seytbootno use setting, only: infile,n,outfile,smtype implicit none ! Outputs u, y, hy, seytbootno, y + seytbootno, y - seytbootno. integer :: j integer :: nout ! for output integer :: openerror real(sp), allocatable, dimension(:) :: uout ! for output real(sp), allocatable, dimension(:) :: yout ! for output real(sp), allocatable, dimension(:) :: hyout ! for output real(sp), allocatable, dimension(:) :: seytbootnoout ! for output real(sp), allocatable, dimension(:) :: ylow1out ! for output real(sp), allocatable, dimension(:) :: yupp1out ! for output real(sp), allocatable, dimension(:) :: ylow2out ! for output real(sp), allocatable, dimension(:) :: yupp2out ! for output real(sp), allocatable, dimension(:) :: tout ! for output real(sp), allocatable, dimension(:) :: xout ! for output real(sp), allocatable, dimension(:) :: xres ! residual, for output real(sp), parameter :: zero=0.0_sp real(sp), parameter :: two=2.0_sp ! ! 1. Open file ! ========= ! open (unit=2, file='result/'//outfile, status='unknown', & form='formatted', action='write', iostat=openerror) if (openerror /= 0 ) then print *,'output file: ',outfile,' Error during output file opening---please check' return end if ! ! 2. Write header ! ============ ! write (unit=2, fmt='(1x,5(a))') & 'Output file (directory ''/result/'')',',',' kernel 3.2 (January 2017)',',', & ' www.climate-risk-analysis.com' write (unit=2, fmt='(1x,a)') & ' ' write (unit=2, fmt='(1x,a,a)') & 'Input data (directory ''/data/'') = ',infile if (smtype .eq. 'g') then write (unit=2, fmt='(1x,a,a)') & 'Smoothing type = ','global' if (ismog .eq. 1) then write (unit=2, fmt='(1x,a,f0.6)') & 'Predefined global bandwidth = ',hyglobal else if (ismog .eq. 0) then write (unit=2, fmt='(1x,a,f0.6)') & 'Optimized global bandwidth = ',hyglobal end if else if (smtype .eq. 'l') then write (unit=2, fmt='(1x,a,a)') & 'Smoothing type = ','local' end if if (nue .eq. 0) then write (unit=2, fmt='(1x,a,a)') & 'Order of derivative: = ','0th (trend)' else if (nue .eq. 1) then write (unit=2, fmt='(1x,a,a)') & 'Order of derivative: = ','1st (slope)' else if (nue .eq. 2) then write (unit=2, fmt='(1x,a,a)') & 'Order of derivative: = ','2nd (curvature)' end if write (unit=2, fmt='(1x,a,f0.6)') & 'Predefined persistence time = ',tau write (unit=2, fmt='(1x,a,i0.0)') & 'MBB block length = ',l write (unit=2, fmt='(1x,a)') & ' ' write (unit=2, fmt='(1x,a)') & 'Time = time (kernel estimate)' write (unit=2, fmt='(1x,a,a,a)') & 'Estimate = kernel estimate (trend',',',' slope or curvature)' write (unit=2, fmt='(1x,a)') & 'Bandwidth = optimized or predefined local or global bandwidth' write (unit=2, fmt='(1x,a)') & 'se = standard error (estimate)' write (unit=2, fmt='(1x,a)') & 'Lower-1se = estimate - 1 * se' write (unit=2, fmt='(1x,a)') & 'Upper-1se = estimate + 1 * se' write (unit=2, fmt='(1x,a)') & 'Lower-2se = estimate - 2 * se' write (unit=2, fmt='(1x,a)') & 'Upper-2se = estimate + 2 * se' write (unit=2, fmt='(1x,a)') & 'Timedata = time input data' write (unit=2, fmt='(1x,a)') & 'Data = input data' write (unit=2, fmt='(1x,a)') & 'Residual = data - estimate (defined for time = timedata)' write (unit=2, fmt='(1x,a)') & ' ' write (unit=2, fmt='(1x,a)') & 'Value -999.0 means N/A ' write (unit=2, fmt='(1x,a)') & ' ' write (unit=2, fmt='(1x,11(x,a))') & 'Time ', & 'Estimate ', & 'Bandwidth ', & 'se ', & 'Lower-1se ', & 'Upper-1se ', & 'Lower-2se ', & 'Upper-2se ', & 'Timedata ', & 'Data ', & 'Residual' ! ! 3. Define output data ! ================== ! if (nu .gt. n) then nout=nu allocate(uout(nout)) allocate(yout(nout)) allocate(hyout(nout)) allocate(seytbootnoout(nout)) allocate(ylow1out(nout)) allocate(yupp1out(nout)) allocate(ylow2out(nout)) allocate(yupp2out(nout)) allocate(tout(nout)) allocate(xout(nout)) allocate(xres(nout)) uout=u yout=y if (smtype .eq. 'g') then hyout=hyglobal else if (smtype .eq. 'l') then hyout=hy end if seytbootnoout=seytbootno ylow1out=yout-seytbootno yupp1out=yout+seytbootno ylow2out=yout-two*seytbootno yupp2out=yout+two*seytbootno tout(1:n)=t(1:n) xout(1:n)=x(1:n) tout(n+1:nu)=-999.0_sp xout(n+1:nu)=-999.0_sp xres=-999.0_sp else if (nu .le. n) then nout=n allocate(uout(nout)) allocate(yout(nout)) allocate(hyout(nout)) allocate(seytbootnoout(nout)) allocate(ylow1out(nout)) allocate(yupp1out(nout)) allocate(ylow2out(nout)) allocate(yupp2out(nout)) allocate(tout(nout)) allocate(xout(nout)) allocate(xres(nout)) uout(1:nu)=u(1:nu) yout(1:nu)=y(1:nu) uout(nu+1:n)=-999.0_sp yout(nu+1:n)=-999.0_sp if (smtype .eq. 'g') then hyout(1:nu)=hyglobal hyout(nu+1:n)=-999.0_sp else if (smtype .eq. 'l') then hyout(1:nu)=hy(1:nu) hyout(nu+1:n)=-999.0_sp end if seytbootnoout(1:nu)=seytbootno(1:nu) seytbootnoout(nu+1:n)=-999.0_sp ylow1out(1:nu)=yout(1:nu)-seytbootno(1:nu) yupp1out(1:nu)=yout(1:nu)+seytbootno(1:nu) ylow2out(1:nu)=yout(1:nu)-two*seytbootno(1:nu) yupp2out(1:nu)=yout(1:nu)+two*seytbootno(1:nu) ylow1out(nu+1:n)=-999.0_sp yupp1out(nu+1:n)=-999.0_sp ylow2out(nu+1:n)=-999.0_sp yupp2out(nu+1:n)=-999.0_sp tout=t xout=x xres=-999.0_sp if (nu .eq. n) then if (maxval(abs(u(1:nu)-t(1:n))) .eq. zero) xres=x-y ! then define residuals end if end if ! ! 4. Write data ! ========== ! do j=1,nout write (unit=2, fmt='(11(1x,e14.6e3))') & uout(j),yout(j),hyout(j), & seytbootnoout(j),ylow1out(j),yupp1out(j),ylow2out(j),yupp2out(j), & tout(j),xout(j),xres(j) end do close (unit=2, status='keep') deallocate(uout) deallocate(yout) deallocate(hyout) deallocate(seytbootnoout) deallocate(ylow1out) deallocate(yupp1out) deallocate(ylow2out) deallocate(yupp2out) deallocate(tout) deallocate(xout) deallocate(xres) end subroutine output ! !======================================================================================================================== ! subroutine ranseed use nrtype use nr, only: ran implicit none ! Seeds random number generator. integer :: seed=0 real(sp) :: dummy=-999.0_sp seed=-206761862 dummy=ran(seed) end subroutine ranseed ! !======================================================================================================================== ! subroutine read1 use nrtype use data1 use setting, only: infile,n implicit none ! Reads t, x. integer :: i real(sp), parameter :: one=1.0_sp open (unit=1, file='data/'//infile, status='old', & form='formatted', action='read') do i=1,n read (unit=1, fmt=*) t(i),x(i) t(i)=t(i)*one x(i)=x(i)*one if (i .gt. 1 .and. t(i) .le. t(i-1)) & print '(a,i4,a,i4,a,a)', & 'Subroutine read1 --- nonincreasing t(',i,') of ',infile end do close (unit=1, status='keep') end subroutine read1 ! !======================================================================================================================== ! subroutine smooth use nrtype use data1 use kernelbandwidth use kernelparameters use kernelwork use parameters, only: nb,nu,u use persist use result1 use setting, only: n,smtype,infile implicit none ! ! Local/global smoothing with nonparametric regression using Gasser-M\"{u}ller kernel. ! ! IMPORTANT NOTE: before calling smoothing subroutines (glkern, lokern), ! array s has to be re-initialized by calling subroutine inits; and ! variables til and tiu have to be re-initialized by calling subroutine initti. ! ! Determines y, seytbootno. ! ! References: Fishman GS (1996) Monte Carlo: Concepts, Algorithms, and Applications. Springer, New York, 698 pp. ! Gasser T, M\"{u}ller H-G (1979) In: Smoothing Techniques for Curve Estimation, Eds. ! Gasser T, Rosenblatt M, Springer, Berlin, 23-68. ! Gasser T, M\"{u}ller H-G (1984) Scandinavian Journal of Statistics 11:171-185. ! Gasser T, Kneip A, K\"{o}hler W (1991) Journal of the American Statistical Association 86:643-652. ! Herrmann E (1997) Journal of Computational and Graphical Statistics 6:35-54. ! Mudelsee M (2010) Climate Time Series Analysis: Classical Statistical and Bootstrap Methods. ! Springer, Dordrecht, 474 pp. ! Simonoff J (1996) Smoothing Methods in Statistics. Springer, New York, 338 pp. ! ! URL (FORTRAN 77 subroutines): http://www.biostat.unizh.ch/research/software/kernel.html (3 May 2007) ! integer :: b=0 ! counter (bootstrap resampling) integer :: iresult=-1 ! result (exclusion test; '0': passed) integer :: itest=0 ! counter (exclusion test) integer :: j=0 ! counter real(sp), dimension(nu) :: sybootno ! updating ave (Fishman 1996: page 67f therein), t* no real(sp), dimension(nu) :: tybootno ! updating var (Fishman 1996: page 67f therein), t* no real(sp), parameter :: zero=0.0_sp real(sp), parameter :: one=1.0_sp real(sp), parameter :: two=2.0_sp ! ! 1. Bandwidth optimization done on 0th (trend), time grid u ! ======================================================= ! ! Global smoothing: determine hyglobal (only in case of ismog = 0) ! Local smoothing: determine hy (use ismo = 0) ! ! For both: use nuex = 0 (0th or trend) and kordx = nuex + 2 = 2 ! if (smtype .eq. 'g') then if (ismog .eq. 0) then call inits call initti call glkern(t(1:n),x(1:n),n,u(1:nu),nu,ihom,nuex,kordx,irnd,ismog,m1,til,tiu, & s(0:n),sig,wn(0:n,1:5),w1(1:m1,1:3), & hyglobal,y(1:nu)) end if else if (smtype .eq. 'l') then call inits call initti call lokern(t(1:n),x(1:n),n,u(1:nu),nu,ihom,nuex,kordx,irnd,ismo,m1,til,tiu, & s(0:n),sig,wn(0:n,1:5),w1(1:m1,1:3), & wmy(1:nu),hy(1:nu),y(1:nu)) end if ! ! 2. Kernel estimation of 0th/1st/2nd (y), time grid u ! ================================================= ! ! Global smoothing: overtake hyglobal (use ismofix = 1) ! Local smoothing: overtake hy (use ismofix = 1) ! if (smtype .eq. 'g') then call inits call initti call glkern(t(1:n),x(1:n),n,u(1:nu),nu,ihom,nue,kord,irnd,ismofix,m1,til,tiu, & s(0:n),sig,wn(0:n,1:5),w1(1:m1,1:3), & hyglobal,y(1:nu)) else if (smtype .eq. 'l') then call inits call initti call lokern(t(1:n),x(1:n),n,u(1:nu),nu,ihom,nue,kord,irnd,ismofix,m1,til,tiu, & s(0:n),sig,wn(0:n,1:5),w1(1:m1,1:3), & wmy(1:nu),hy(1:nu),y(1:nu)) end if ! ! ! 3. Persistence time (option 1)/standard error (option 2) determination (seytbootno), time grid u ! ================================================================================================== ! ! ! 3.1 Kernel estimation of 0th (xtrend), timescale t ! ============================================== ! ! Global smoothing: use ismofix = 1 and overtake hyglobal ! Local smoothing: use ismo = 0 and it does not matter that hx is overwritten, ! the assumption is that the smoothing solution ! is the same for same data, t(i), x(i), ! regardless where output, t(i) or u(i), is done; ! validity of this assumption has been verified ! using test below (result: xtrend agrees with y) ! ! For both: use nuex = 0 (0th or trend) and kordx = nuex + 2 = 2 ! if (smtype .eq. 'g') then call inits call initti call glkern(t(1:n),x(1:n),n,t(1:n),n,ihom,nuex,kordx,irnd,ismofix,m1,til,tiu, & s(0:n),sig,wn(0:n,1:5),w1(1:m1,1:3), & hyglobal,xtrend(1:n)) else if (smtype .eq. 'l') then call inits call initti call lokern(t(1:n),x(1:n),n,t(1:n),n,ihom,nuex,kordx,irnd,ismo,m1,til,tiu, & s(0:n),sig,wn(0:n,1:5),w1(1:m1,1:3), & wmx(1:n),hx(1:n),xtrend(1:n)) end if ! ! ! 3.2 Unweighted residuals ! ==================== ! e(1:n)=x(1:n)-xtrend(1:n) ! ! ! 3.3 Option 1: Persistence time estimation ! ===================================== ! if (tauflag .eq. 'y') then open (unit=3, file=taufile, status='unknown', & form='formatted', action='write') do j=1,n write (unit=3, fmt='(2(2x,f20.6))') t(j),e(j) end do close (unit=3, status='keep') print '(a,a,a,i4)', 'TAUEST data file written into directory /tauest/ for ',infile,' with n = ',n print '(a)', ' ' print '(a)', 'program kernel terminates now' stop end if ! ! ! 3.4 Option 2: Bootstrap loop (start) ! ================================ ! if (tauflag .eq. 'n') then ! ! ! 3.4.1 Initialize updating variables ! ============================= ! sybootno(1:nu)=zero tybootno(1:nu)=zero ! ! ! 3.4.2 Bootstrap simulation # 1 ! ======================== ! ! Use ismob = 1, that is, ! do not change hyglobal (global smoothing), ! do not change hy (local smoothing) ! b=1 call bootstrap1(n,e(1:n),l,xboot(1:n)) xboot(1:n)=xtrend(1:n)+xboot(1:n) if (smtype .eq. 'g') then call inits call initti call glkern(t(1:n),xboot(1:n),n,u(1:nu),nu,ihom,nue,kord,irnd,ismob,m1,til,tiu, & s(0:n),sig,wn(0:n,1:5),w1(1:m1,1:3), & hyglobal,ybootno(1:nu)) else if (smtype .eq. 'l') then call inits call initti call lokern(t(1:n),xboot(1:n),n,u(1:nu),nu,ihom,nue,kord,irnd,ismob,m1,til,tiu, & s(0:n),sig,wn(0:n,1:5),w1(1:m1,1:3), & wmy(1:nu),hy(1:nu),ybootno(1:nu)) end if sybootno(1:nu)=ybootno(1:nu) ! ! ! 3.4.3 Bootstrap simulations > 1 ! ========================= ! ! Use ismob = 1, that is, ! do not change hyglobal (global smoothing), ! do not change hy (local smoothing) ! do b=2,nb call bootstrap1(n,e(1:n),l,xboot(1:n)) xboot(1:n)=xtrend(1:n)+xboot(1:n) if (smtype .eq. 'g') then call inits call initti call glkern(t(1:n),xboot(1:n),n,u(1:nu),nu,ihom,nue,kord,irnd,ismob,m1,til,tiu, & s(0:n),sig,wn(0:n,1:5),w1(1:m1,1:3), & hyglobal,ybootno(1:nu)) else if (smtype .eq. 'l') then call inits call initti call lokern(t(1:n),xboot(1:n),n,u(1:nu),nu,ihom,nue,kord,irnd,ismob,m1,til,tiu, & s(0:n),sig,wn(0:n,1:5),w1(1:m1,1:3), & wmy(1:nu),hy(1:nu),ybootno(1:nu)) end if tybootno(1:nu)=tybootno(1:nu)+(one-one/b)*(ybootno(1:nu)-sybootno(1:nu)/(b-one))**two sybootno(1:nu)=sybootno(1:nu)+ybootno(1:nu) ! ! ! 3.4.4 Option 2: Bootstrap loop (end) ! ============================== ! end do ! ! ! 3.4.5 Calculation of seytbootno ! ========================= ! seytbootno(1:nu)=sqrt(tybootno(1:nu)/(nb-one)) end if ! !-test-start------------------------------------------------------------------------------------------------------------- ! open (unit=2, file='test.dat', status='unknown', & ! form='formatted', action='write') ! write (unit=2, fmt='(a)') 'u y hy seytbootno t x xtrend e xboot hx' ! do j=1,nu ! if (j .le. n) then ! write (unit=2, fmt='(10(2x,f20.6))') & ! u(j),y(j),hy(j), & ! seytbootno(j), & ! t(j),x(j),xtrend(j),e(j),xboot(j),hx(j) ! else ! write (unit=2, fmt='(4(2x,f20.6))') & ! u(j),y(j),hy(j), & ! seytbootno(j) ! end if ! end do ! close (unit=2, status='keep') ! print '(a,a)', 'test.dat file written for: ',infile ! read (*,'()') !-test-end--------------------------------------------------------------------------------------------------------------- ! end subroutine smooth ! !======================================================================================================================== ! subroutine welcome use persist implicit none ! Welcomes. print '(a)', ' ' print '(a)', ' ' print '(a)', ' ' print '(a)', ' kernel' print '(a)', ' ======' print '(a)', ' ' print '(a)', ' Version 3.2 (January 2017)' print '(a)', ' ' print '(a)', ' ' print '(a)', ' ' print '(a)', ' Author (Fortran 90 part): Manfred Mudelsee' print '(a)', ' Climate Risk Analysis' print '(a)', ' www.climate-risk-analysis.com' print '(a)', ' ' print '(a)', ' ' print '(a)', ' ' if (tauflag .eq. 'y') then print '(a)', ' Current mode: persistence time estimation' else if (tauflag .eq. 'n') then print '(a)', ' Current mode: kernel estimation' end if print '(a)', ' ' print '(a)', ' To change mode (persistence time vs kernel estimation):' print '(a)', ' Modify module ''persist'' and recompile' print '(a)', ' ' print '(a)', ' ' print '(a)', ' ' end subroutine welcome ! !======================================================================================================================== ! ! FORTRAN 77 subroutines ! !======================================================================================================================== ! SUBROUTINE GLKERN(T,X,N,TT,M,IHOM,NUE,KORD,IRND, & ISMO,M1,TL,TU,S,SIG,WN,W1,B,Y) !----------------------------------------------------------------------* ! SHORT-VERSION: OCT 1996 ! ! PURPOSE: ! ! GENERAL SUBROUTINE FOR KERNEL SMOOTHING: ! COMPUTATION OF ITERATIVE PLUG-IN ALGORITHM FOR GLOBAL BANDWIDTH ! SELECTION FOR KERNELS WITH (NUE,KORD) = (0,2),(0,4),(1,3) OR ! (2,4). ! !----------------------------------------------------------------------* ! THE RAW DATA SHOULD BE GIVEN BY THE POINTS ! (T(1),X(1)),...,(T(N),X(N)) ! ! THE RESULTING ESTIMATOR OF THE NUE-TH DERIVATIVE OF THE ! REGRESSION CURVE IS GIVEN THROUGH THE POINTS ! (TT(1),Y(1)),...,(TT(M),Y(M)) ! ! THE PLUG-IN BANDWIDTH IS GIVEN BY B !----------------------------------------------------------------------* ! ! PARAMETERS : ! ! INPUT T(N) INPUT GRID (T(1)<T(2)<...<T(N)) ! INPUT X(N) DATA ! INPUT N LENGTH OF X ! ! INPUT TT(M) OUTPUT GRID, SHOULD BE ORDERED ! INPUT M LENGTH OF TT ! ! INPUT IHOM HOMOSCEDASTICITY OF VARIANCE ! 0: HOMOSZEDASTIC ERROR VARIABLES, ! <> 0: IF THE VARIANCE SHOULD ESTIMATED AS ! SMOOTH FUNCTION. ! ****** DEFAULT VALUE: IHOM=0 ! ! INPUT NUE ORDER OF DERIVATIVE (0-4) OF THE REGRESSION ! FUNCTION WHICH SHALL BE ESTIMATED ! ****** DEFAULT VALUE: NUE=0 ! ! INPUT KORD ORDER OF KERNEL (<=6), FOR ISMO=0 ONLY ! NUE=0, KORD=2 OR KORD=4 ! OR NUE=1, KORD=3 OR NUE=2, KORD=4 ARE ALLOWED ! ****** DEFAULT VALUE: KORD=NUE+2 ! ! INPUT IRND 0: IF RANDOM GRID POINTS T MAY OCCUR ! <>0 ELSE (ONLY NECESSARY IF S SHOULD BE ! COMPUTED) ! ****** DEFAULT VALUE IRND=0 ! ! INPUT ISMO 0:ESTIMATING THE OPTIMAL GLOBAL BANDWIDTH ! <>0 USING GLOBAL INPUT BANDWIDTH B ! ****** DEFAULT VALUE ISMO=0 ! ! INPUT M1 >=10, LENGTH OF W1, LARGE VALUES WILL INCREASE ! THE ACCURACY OF THE INTEGRAL APPROXIMATION ! ****** DEFAULT VALUE: M1=400 ! !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! IN/OUTPUT TL/TU LOWER/UPPER BOUND FOR INTEGRAL APPROXIMATION ! AND VARIANCE ESTIMATION (IF SIG=0 AND IHOM=0), ! IF TU<=TL, [TL,TU] ARE COMPUTED AS ABOUT ! THE 87% MIDDLE PART OF [T(1),T(N)] ! ****** DEFAULT VALUES: TL=1.0, TU=0.0 ! ! IN/OUTPUT S(0:N) IF S(N)<=S(0) THIS ARRAY IS COMPUTED AS ! MIDPOINTS OF T, FOR NON-RANDOM DESIGN AND AS ! SMOOTHED QUANTILES FOR RANDOM DESIGN ! ****** DEFAULT VALUES: S(0)=1.0, S(N)=0.0 ! AND THE OTHER S(I) CAN BE UNDEFINED ! ! IN/OUTPUT SIG RESIDUAL VARIANCE, ESTIMATED FOR SIG=0 OR ! IHOM<>0, ELSE GIVEN BY INPUT ! ****** DEFAULT VALUE: SIG=-1.0 ! ! IN/OUTPUT B GLOBAL PLUG-IN BANDWIDTH ! ****** B CAN BE UNDIFINED IF ISMO=0 ! ! ! WORK WN(0:N,5) WORK ARRAY FOR KERNEL SMOOTHING ROUTINE ! OR NUE=1, KORD=3 OR NUE=2, KORD=4 ARE ALLOWED ! ****** WILL BE SET IN SUBROUTINE ! WORK W1(M1,3) WORK ARRAY FOR INTEGRAL APPROXIMATION ! ****** WILL BE SET IN SUBROUTINE ! !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! OUTPUT Y(M) KERNEL ESTIMATE WITH BOP(=B0 FOR ISMO<>0) ! ****** WILL BE SET IN SUBROUTINE !----------------------------------------------------------------------- ! USED SUBROUTINES: COFF, RESEST, KERNEL WITH FURTHER SUBROUTINES ! WHICH ARE CONTAINED IN THE FILE subs.f !----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N),T(N),TT(M),Y(M),WN(0:N,5),S(0:N) DIMENSION W1(M1,3) DIMENSION BIAS(2,0:2),VARK(2,0:2),FAK2(2:4) !- !-------- 1. INITIALISATIONS AND SOME ERROR-CHECKS DATA BIAS/.2,.04762,.4286,.1515,1.33,.6293/ DATA VARK/.6,1.250,2.143,11.93,35.0,381.6/ DATA FAK2/4.,36.,576./ NYG=0 INPUTS=0 !-------- IF NO ERRORS SHOULD BE WRITTEN ON STANDARD OUTPUT, SET IPRINT=1 !-------- IF ERRORS AND VERY DETAILED WARNINGS SHOULD BE WRITTEN ON !-------- STANDARD OUTPUT, SET IPRINT < 0 IPRINT=-1 ! IF(NUE.GT.4.OR.NUE.LT.0) THEN IF(IPRINT.EQ.0) & PRINT *,'glkern: Order of derivative not allowed' STOP END IF IF(NUE.GT.2.AND.ISMO.EQ.0) THEN IF(IPRINT.EQ.0) & PRINT *,'glkern: Order of derivative not allowed' STOP END IF IF(N.LE.2) THEN IF(IPRINT.EQ.0) & PRINT *,'glkern: Number of data too small' STOP END IF IF(M.LT.1) THEN IF(IPRINT.EQ.0) & PRINT *,'glkern: No output points' STOP END IF IF(M1.LT.10) THEN IF(IPRINT.EQ.0) & PRINT *,'glkern: Variable M1 is choosen too small' STOP END IF KK=(KORD-NUE)/2 IF(2*KK+NUE.NE.KORD) THEN IF(IPRINT.EQ.0) & PRINT *,'glkern: Kernel order not allowed, set to ',NUE+2 KORD=NUE+2 END IF IF(KORD.GT.4.AND.ISMO.EQ.0) THEN IF(IPRINT.EQ.0) & PRINT *,'glkern: Kernel order not allowed, set to ',NUE+2 KORD=NUE+2 END IF IF(KORD.GT.6.OR.KORD.LE.NUE) THEN IF(IPRINT.EQ.0) & PRINT *,'glkern: Kernel order not allowed, set to ',NUE+2 KORD=NUE+2 END IF IF(ISMO.NE.0.AND.B.LE.0) THEN IF(IPRINT.EQ.0) & PRINT *,'glkern: Plug-in bandwidth is used' ISMO=0 END IF RVAR=SIG !- !-------- 2. COMPUTATION OF S-SEQUENCE S0=1.5*T(1)-0.5*T(2) SN=1.5*T(N)-0.5*T(N-1) IF(S(N).LE.S(0)) THEN INPUTS=1 DO 20 I=1,N-1 20 S(I)=.5*(T(I)+T(I+1)) S(0)=S0 S(N)=SN IF(ISMO.NE.0.AND.IRND.NE.0) GOTO 160 ELSE IF(ISMO.NE.0) GOTO 160 END IF !- !-------- 3. COMPUTATION OF MINIMAL, MAXIMAL ALLOWED BANDWIDTH BMAX=(SN-S0)*.5 BMIN=(SN-S0)/DBLE(N)*DBLE(KORD-1)*.6 !- !-------- 4. WARNINGS IF TT-GRID LARGER THAN T-GRID IF(TT(1).LT.S0.AND.TT(M).GT.SN.AND.IPRINT.LT.0) PRINT *, & 'glkern: Extrapolation at both boundaries not optimized' IF(TT(1).LT.S0.AND.TT(M).LE.SN.AND.IPRINT.LT.0) PRINT *, & 'glkern: Extrapolation at left boundary not optimized' IF(TT(1).GE.S0.AND.TT(M).GT.SN.AND.IPRINT.LT.0) PRINT *, & 'glkern: Extrapolation at right boundary not optimized' !- !-------- 5. COMPUTE TL,TU AND THEIR T-GRID AS INNER PART FOR ! INTEGRAL APPROXIMATION IN THE ITERATIONS ITT=0 51 IF (TU.LE.TL) THEN TL=.933*S0+.067*SN TU=.067*S0+.933*SN ITT=ITT+1 END IF TL=MAX(TL,S0) TU=MIN(TU,SN) IL=1 IU=N WN(1,1)=0.0 WN(N,1)=0.0 DO 50 I=1,N IF(T(I).LE.TL.OR.T(I).GE.TU) WN(I,1)=0.0 IF(T(I).GT.TL.AND.T(I).LT.TU) WN(I,1)=1.0 IF(T(I).LT.TL) IL=I+1 50 IF(T(I).LE.TU) IU=I NN=IU-IL+1 IF(NN.EQ.0.AND.ITT.EQ.0) THEN TU=TL-1.0 GOTO 51 END IF IF(NN.EQ.0.AND.ITT.EQ.1) THEN TU=SN TL=S0 GOTO 51 END IF !- !-------- 6. COMPUTE T-GRID FOR INTEGRAL APPROXIMATION DO 60 I=1,M1 W1(I,2)=1.0 60 W1(I,1)=TL+(TU-TL)*DBLE(I-1)/DBLE(M1-1) !- !-------- 7. CALCULATION OF WEIGHT FUNCTION ALPHA=1.D0/DBLE(13) DO 70 I=IL,IU XI=(T(I) - TL)/ALPHA/(TU-TL) IF(XI.GT.1) GOTO 71 70 WN(I,1)=(10.0-15*XI+6*XI*XI)*XI*XI*XI 71 DO 72 I=IU,IL,-1 XI=(TU-T(I))/ALPHA/(TU-TL) IF(XI.GT.1) GOTO 73 72 WN(I,1)=(10.0-15*XI+6*XI*XI)*XI*XI*XI 73 DO 74 I=1,M1 XI=(W1(I,1)-TL)/ALPHA/(TU-TL) IF(XI.GT.1) GOTO 75 74 W1(I,2)=(10.0-15*XI+6*XI*XI)*XI*XI*XI 75 DO 76 I=M1,1,-1 XI=(TU-W1(I,1))/ALPHA/(TU-TL) IF(XI.GT.1) GOTO 77 76 W1(I,2)=(10.0-15*XI+6*XI*XI)*XI*XI*XI 77 CONTINUE !- !-------- 8. COMPUTE CONSTANTS FOR ITERATION EX=1./DBLE(KORD+KORD+1) KK2=(KORD-NUE) KK=KK2/2 !- !-------- 9. ESTIMATING VARIANCE AND SMOOTHED PSEUDORESIDUALS IF((SIG.LE..0).AND.(IHOM.EQ.0)) & CALL RESEST(T(IL),X(IL),NN,WN(IL,2),R2,SIG) IF(IHOM.NE.0) THEN CALL RESEST(T,X,N,WN(1,2),SNR,SIG) BRES=MAX(BMIN,.2*NN**(-.2)*(S(IU)-S(IL-1))) DO 91 I=1,N WN(I,3)=T(I) 91 WN(I,2)=WN(I,2)*WN(I,2) CALL KERNEL(T,WN(1,2),N,BRES,0,KK2,NYG,S, & WN(IL,3),NN,WN(IL,4)) ELSE CALL COFF(WN(1,4),N,SIG) END IF !- !-------- 10. ESTIMATE/COMPUTE INTEGRAL CONSTANT 100 VI=0. DO 101 I=IL,IU 101 VI=VI+WN(I,1)*N*(S(I)-S(I-1))**2*WN(I,4) !- !-------- 11. REFINEMENT OF S-SEQUENCE FOR RANDOM DESIGN IF(INPUTS.EQ.1.AND.IRND.EQ.0) THEN DO 110 I=0,N WN(I,5)=DBLE(I)/DBLE(N+1) WN(I,2)=(DBLE(I)+.5)/DBLE(N+1) 110 WN(I,3)=WN(I,2) EXS=-DBLE(3*KORD+1)/DBLE(6*KORD+3) EXSVI=DBLE(KORD)/DBLE(6*KORD+3) BS=0.1*(VI/(SN-S0)**2)**EXSVI*N**EXS CALL KERNEL(WN(1,5),T,N,BS,0,2,NYG,WN(0,3),WN(0,2),N+1,S(0)) 111 ISORT=0 VI=0.0 DO 112 I=1,N VI=VI+WN(I,1)*N*(S(I)-S(I-1))**2*WN(I,4) IF(S(I).LT.S(I-1)) THEN SSI=S(I-1) S(I-1)=S(I) S(I)=SSI ISORT=1 END IF 112 CONTINUE IF(ISORT.EQ.1) GOTO 111 IF(ISMO.NE.0) GOTO 160 END IF B=BMIN*2. !- !-------- 12. COMPUTE INFLATION CONSTANT AND EXPONENT AND LOOP OF ITERATIONS CONST=DBLE(2*NUE+1)*FAK2(KORD)*VARK(KK,NUE)*VI & /(DBLE(2*KORD-2*NUE)*BIAS(KK,NUE)**2*DBLE(N)) FAC=1.1*(1.+(NUE/10.)+0.05*(KORD-NUE-2.)) & *N**(2./DBLE((2*KORD+1)*(2*KORD+3))) ITENDE=1+2*KORD+KORD*(2*KORD+1) DO 120 IT=1,ITENDE !- !-------- 13. ESTIMATE DERIVATIVE OF ORDER KORD IN ITERATIONS B2=B*FAC B2=MAX(B2,BMIN/DBLE(KORD-1)*DBLE(KORD+1)) B2=MIN(B2,BMAX) CALL KERNEL(T,X,N,B2,KORD,KORD+2,NYG,S,W1(1,1),M1,W1(1,3)) !- !-------- 14. ESTIMATE INTEGRALFUNCTIONAL IN ITERATIONS XMY2=.75*(W1(1,2)*W1(1,3)*W1(1,3)+W1(M1,2)*W1(M1,3)*W1(M1,3)) DO 140 I=2,M1-1 140 XMY2=XMY2+W1(I,2)*W1(I,3)*W1(I,3) XMY2=XMY2*(TU-TL)/DBLE(M1) !- !-------- 15. FINISH OF ITERATIONS B=(CONST/XMY2)**EX B=MAX(BMIN,B) B=MIN(BMAX,B) 120 CONTINUE !- !-------- 16 COMPUTE SMOOTHED FUNCTION WITH PLUG-IN BANDWIDTH 160 CALL KERNEL(T,X,N,B,NUE,KORD,NYG,S,TT,M,Y) !- !-------- 17. VARIANCE CHECK IF(IHOM.NE.0) SIG=RVAR IF(RVAR.EQ.SIG.OR.R2.LT..88.OR.IHOM.NE.0.OR.NUE.GT.0) RETURN II=0 IIL=0 J=2 TLL=MAX(TL,TT(1)) TUU=MIN(TU,TT(M)) DO 170 I=IL,IU IF(T(I).LT.TLL.OR.T(I).GT.TUU) GOTO 170 II=II+1 IF(IIL.EQ.0) IIL=I 171 IF(TT(J).LT.T(I)) THEN J=J+1 IF(J.LE.M) GOTO 171 END IF WN(II,3)=X(I)-Y(J)+(Y(J)-Y(J-1))*(TT(J)-T(I))/(TT(J)-TT(J-1)) 170 CONTINUE IF(IIL.EQ.0.OR.II-IIL.LT.10) THEN CALL RESEST(T(IL),WN(1,3),NN,WN(1,4),SNR,RVAR) ELSE CALL RESEST(T(IIL),WN(1,3),II,WN(1,4),SNR,RVAR) END IF Q=SIG/RVAR IF(Q.LE.2.) RETURN IF(Q.GT.5..AND.R2.GT..95) RVAR=RVAR*.5 SIG=RVAR CALL COFF(WN(1,4),N,SIG) GOTO 100 END SUBROUTINE LOKERN(T,X,N,TT,M,IHOM,NUE,KORD,IRND, & ISMO,M1,TL,TU,S,SIG,WN,W1,WM,BAN,Y) !----------------------------------------------------------------------* !----------------------------------------------------------------------- ! SHORT-VERSION: JANUARY 1997 ! ! PURPOSE: ! ! GENERAL SUBROUTINE FOR KERNEL SMOOTHING: ! COMPUTATION OF ITERATIVE PLUG-IN ALGORITHM FOR LOCAL BANDWIDTH ! SELECTION FOR KERNELS WITH (NUE,KORD) = (0,2),(0,4),(1,3) OR ! (2,4). ! !----------------------------------------------------------------------* ! THE RAW DATA SHOULD BE GIVEN BY THE POINTS ! (T(1),X(1)),...,(T(N),X(N)) ! ! THE RESULTING ESTIMATOR OF THE NUE-TH DERIVATIVE OF THE ! REGRESSION CURVE IS GIVEN THROUGH THE POINTS ! (TT(1),Y(1)),...,(TT(M),Y(M)) ! ! THE LOCAL PLUG-IN BANDWIDTH ARRAY IS GIVEN BY BAN(1),...,BAN(M) !----------------------------------------------------------------------* ! ! ! PARAMETERS : ! ! INPUT T(N) INPUT GRID (T(1)<T(2)<...<T(N)) ! INPUT X(N) DATA ! INPUT N LENGTH OF X ! ! INPUT TT(M) OUTPUT GRID ! INPUT M LENGTH OF TT ! ! INPUT IHOM HOMOSCEDASTICITY OF VARIANCE ! 0: HOMOSCEDASTIC ERROR VARIABLES ! <> 0: IF THE VARIANCE SHOULD ESTIMATED AS ! SMOOTH FUNCTION. ! ****** DEFAULT VALUE: IHOM=0 ! ! INPUT NUE ORDER OF DERIVATIVE (0-4) ! ****** DEFAULT VALUE: NUE=0 ! ! INPUT KORD ORDER OF KERNEL (<=6), FOR ISMO=0 ONLY ! NUE=0, KORD=2 OR KORD=4 ! OR NUE=1, KORD=3 OR NUE=2, KORD=4 ARE ALLOWED ! ****** DEFAULT VALUE: KORD=NUE+2 ! ! INPUT IRND 0: IF RANDOM GRID POINTS T MAY OCCUR ! <>0 ELSE (ONLY NECESSARY IF S SHOULD BE ! COMPUTED) ! ****** DEFAULT VALUE IRND=0 ! ! INPUT ISMO 0:ESTIMATING THE OPTIMAL LOCAL BANDWIDTH ! <>0 USING LOCAL INPUT BANDWIDTH-ARRAY IN BAN ! ****** DEFAULT VALUE ISMO=0 ! ! INPUT M1 >=10, LENGTH OF W1, LARGE VALUES WILL INCREASE ! THE ACCURACY OF THE INTEGRAL APPROXIMATION ! ****** DEFAULT VALUE: M1=400 ! !----------------------------------------------------------------------- ! IN/OUTPUT TL/TU LOWER/UPPER BOUND FOR INTEGRAL APPROXIMATION ! AND VARIANCE ESTIMATION (IF SIG=0 AND IHOM=0), ! IF TU<=TL, [TL,TU] ARE COMPUTED AS ABOUT ! THE 87% MIDDLE PART OF [T(1),T(N)] ! ****** DEFAULT VALUES: TL=1.0, TU=0.0 ! ! IN/OUTPUT S(0:N) IF S(N)<=S(0) THIS ARRAY IS COMPUTED AS ! MIDPOINTS OF T, FOR NON-RANDOM DESIGN AND AS ! SMOOTHED QUANTILES FOR RANDOM DESIGN ! ****** DEFAULT VALUES: S(0)=1.0, S(N)=0.0 ! AND THE OTHER S(I) CAN BE UNDEFINED ! ! IN/OUTPUT SIG RESIDUAL VARIANCE, ESTIMATED FOR SIG=0 OR ! IHOM<>0, ELSE GIVEN BY INPUT ! ****** DEFAULT VALUE: SIG=-1.0 ! ! IN/OUTPUT BAN(M) LOCAL PLUG-IN BANDWIDTH ARRAY ! ****** WILL BE SET IN SUBROUTINE FOR ISMO=0 ! ! WORK WN(0:N,5) WORK ARRAY FOR KERNEL SMOOTHING ROUTINE ! ****** WILL BE SET IN SUBROUTINE ! WORK W1(M1,3) WORK ARRAY FOR INTEGRAL APPROXIMATION ! ****** WILL BE SET IN SUBROUTINE ! WORK WM(M) WORK ARRAY FOR INTEGRAL APPROXIMATION ! ****** WILL BE SET IN SUBROUTINE !----------------------------------------------------------------------- ! OUTPUT Y(M) KERNEL ESTIMATE WITH BOP(=B0 FOR ISMO<>0) ! ****** WILL BE SET IN SUBROUTINE !----------------------------------------------------------------------- ! USED SUBROUTINES: COFF, RESEST, KERNEL WITH FURTHER SUBROUTINES ! WHICH ARE CONTAINED IN THE FILE subs.f !----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N),T(N),TT(M),Y(M),WN(0:N,5),S(0:N) DIMENSION W1(M1,3),BAN(M),WM(M) DIMENSION BIAS(2,0:2),VARK(2,0:2),FAK2(2:4) !- !-------- 1. INITIALISATIONS AND SOME ERROR-CHECKS DATA BIAS/.2,.04762,.4286,.1515,1.33,.6293/ DATA VARK/.6,1.250,2.143,11.93,35.0,381.6/ DATA FAK2/4.,36.,576./ NYG=0 INPUTS=0 !-------- IF NO ERRORS SHOULD BE WRITTEN ON STANDARD OUTPUT, SET IPRINT > 0 !-------- IF ERRORS AND VERY DETAILED WARNINGS SHOULD BE WRITTEN ON !-------- STANDARD OUTPUT, SET IPRINT < 0 IPRINT=+1 ! IF(NUE.GT.4.OR.NUE.LT.0) THEN IF(IPRINT.EQ.0) & PRINT *,'lokern: Order of derivative not allowed' STOP END IF IF(NUE.GT.2.AND.ISMO.EQ.0) THEN IF(IPRINT.EQ.0) & PRINT *,'lokern: Order of derivative not allowed' STOP END IF IF(N.LE.2) THEN IF(IPRINT.EQ.0) & PRINT *,'lokern: Number of data too small' STOP END IF IF(M.LT.1) THEN IF(IPRINT.EQ.0) & PRINT *,'lokern: No output points' STOP END IF IF(M1.LT.10) THEN IF(IPRINT.EQ.0) & PRINT *,'lokern: Variable M1 is choosen too small' STOP END IF KK=(KORD-NUE)/2 IF(2*KK+NUE.NE.KORD) THEN IF(IPRINT.EQ.0) & PRINT *,'lokern: Kernel order not allowed, set to ',NUE+2 KORD=NUE+2 END IF IF(KORD.GT.4.AND.ISMO.EQ.0) THEN IF(IPRINT.EQ.0) & PRINT *,'lokern: Kernel order not allowed, set to ',NUE+2 KORD=NUE+2 END IF IF(KORD.GT.6.OR.KORD.LE.NUE) THEN IF(IPRINT.EQ.0) & PRINT *,'lokern: Kernel order not allowed, set to ',NUE+2 KORD=NUE+2 END IF RVAR=SIG !- !-------- 2. COMPUTATION OF S-SEQUENCE S0=1.5*T(1)-0.5*T(2) SN=1.5*T(N)-0.5*T(N-1) IF(S(N).LE.S(0)) THEN INPUTS=1 DO 20 I=1,N-1 20 S(I)=.5*(T(I)+T(I+1)) S(0)=S0 S(N)=SN IF(ISMO.NE.0.AND.IRND.NE.0) GOTO 230 ELSE IF(ISMO.NE.0) GOTO 230 END IF !- !-------- 3. COMPUTATION OF MINIMAL, MAXIMAL ALLOWED GLOBAL BANDWIDTH BMAX=(SN-S0)*.5 BMIN=(SN-S0)/DBLE(N)*DBLE(KORD-1)*.6 !- !-------- 4. WARNINGS IF TT-GRID LARGER THAN T-GRID IF(TT(1).LT.S0.AND.TT(M).GT.SN.AND.IPRINT.LT.0) PRINT *, & 'lokern: Extrapolation at both boundaries not optimized' IF(TT(1).LT.S0.AND.TT(M).LE.SN.AND.IPRINT.LT.0) PRINT *, & 'lokern: Extrapolation at left boundary not optimized' IF(TT(1).GE.S0.AND.TT(M).GT.SN.AND.IPRINT.LT.0) PRINT *, & 'lokern: Extrapolation at right boundary not optimized' !- !-------- 5. COMPUTE TL,TU AND THEIR T-GRID AS INNER PART FOR ! INTEGRAL APPROXIMATION IN THE ITERATIONS ITT=0 51 IF (TU.LE.TL) THEN TL=.933*S0+.067*SN TU=.067*S0+.933*SN ITT=ITT+1 END IF TL=MAX(S0,TL) TU=MIN(SN,TU) IL=1 IU=N WN(1,1)=0.0 WN(N,1)=0.0 DO 50 I=1,N IF(T(I).LE.TL.OR.T(I).GE.TU) WN(I,1)=0.0 IF(T(I).GT.TL.AND.T(I).LT.TU) WN(I,1)=1.0 IF(T(I).LT.TL) IL=I+1 50 IF(T(I).LE.TU) IU=I NN=IU-IL+1 IF(NN.EQ.0.AND.ITT.EQ.0) THEN TU=TL-1.0 GOTO 51 END IF IF(NN.EQ.0.AND.ITT.EQ.1) THEN TU=SN TL=S0 GOTO 51 END IF !- !-------- 6. COMPUTE T-GRID FOR INTEGRAL APPROXIMATION DO 60 I=1,M1 W1(I,2)=1.0 60 W1(I,1)=TL+(TU-TL)*DBLE(I-1)/DBLE(M1-1) !- !-------- 7. CALCULATION OF WEIGHT FUNCTION ALPHA=1.D0/DBLE(13) DO 70 I=IL,IU XI=(T(I) - TL)/ALPHA/(TU-TL) IF(XI.GT.1) GOTO 71 70 WN(I,1)=(10.0-15*XI+6*XI*XI)*XI*XI*XI 71 DO 72 I=IU,IL,-1 XI=(TU-T(I))/ALPHA/(TU-TL) IF(XI.GT.1) GOTO 73 72 WN(I,1)=(10.0-15*XI+6*XI*XI)*XI*XI*XI 73 DO 74 I=1,M1 XI=(W1(I,1)-TL)/ALPHA/(TU-TL) IF(XI.GT.1) GOTO 75 74 W1(I,2)=(10.0-15*XI+6*XI*XI)*XI*XI*XI 75 DO 76 I=M1,1,-1 XI=(TU-W1(I,1))/ALPHA/(TU-TL) IF(XI.GT.1) GOTO 77 76 W1(I,2)=(10.0-15*XI+6*XI*XI)*XI*XI*XI 77 CONTINUE !- !-------- 8. COMPUTE CONSTANTS FOR ITERATION EX=1./DBLE(KORD+KORD+1) KK2=(KORD-NUE) KK=KK2/2 !- !-------- 9. ESTIMATING VARIANCE AND SMOOTHED PSEUDORESIDUALS IF((SIG.LE..0).AND.(IHOM.EQ.0)) & CALL RESEST(T(IL),X(IL),NN,WN(IL,2),R2,SIG) IF(IHOM.NE.0) THEN CALL RESEST(T,X,N,WN(1,2),SNR,SIG) BRES=MAX(BMIN,.2*NN**(-.2)*(S(IU)-S(IL-1))) DO 91 I=1,N WN(I,3)=T(I) 91 WN(I,2)=WN(I,2)*WN(I,2) CALL KERNEL(T,WN(1,2),N,BRES,0,KK2,NYG,S, & WN(1,3),N,WN(1,4)) ELSE CALL COFF(WN(1,4),N,SIG) END IF !- !-------- 10. ESTIMATE/COMPUTE INTEGRAL CONSTANT 100 VI=0. DO 101 I=IL,IU 101 VI=VI+WN(I,1)*N*(S(I)-S(I-1))**2*WN(I,4) !- !-------- 11. REFINEMENT OF S-SEQUENCE FOR RANDOM DESIGN IF(INPUTS.EQ.1.AND.IRND.EQ.0) THEN DO 110 I=0,N WN(I,5)=DBLE(I)/DBLE(N+1) WN(I,2)=(DBLE(I)+.5)/DBLE(N+1) 110 WN(I,3)=WN(I,2) EXS=-DBLE(3*KORD+1)/DBLE(6*KORD+3) EXSVI=DBLE(KORD)/DBLE(6*KORD+3) BS=0.1*(VI/(SN-S0)**2)**EXSVI*N**EXS CALL KERNEL(WN(1,5),T,N,BS,0,2,NYG,WN(0,3),WN(0,2),N+1,S(0)) 111 ISORT=0 VI=0.0 DO 112 I=1,N VI=VI+WN(I,1)*N*(S(I)-S(I-1))**2*WN(I,4) IF(S(I).LT.S(I-1)) THEN SSI=S(I-1) S(I-1)=S(I) S(I)=SSI ISORT=1 END IF 112 CONTINUE IF(ISORT.EQ.1) GOTO 111 IF(ISMO.NE.0) GOTO 230 END IF B=BMIN*2. !- !-------- 12. COMPUTE INFLATION CONSTANT AND EXPONENT AND LOOP OF ITERATIONS CONST=DBLE(2*NUE+1)*FAK2(KORD)*VARK(KK,NUE)*VI & /(DBLE(2*KORD-2*NUE)*BIAS(KK,NUE)**2*DBLE(N)) FAC=1.1*(1.+(NUE/10.)+0.05*(KORD-NUE-2.)) & *N**(2./DBLE((2*KORD+1)*(2*KORD+3))) ITENDE=1+2*KORD+KORD*(2*KORD+1) DO 120 IT=1,ITENDE !- !-------- 13. ESTIMATE DERIVATIVE OF ORDER KORD IN ITERATIONS B2=B*FAC B2=MAX(B2,BMIN/DBLE(KORD-1)*DBLE(KORD+1)) B2=MIN(B2,BMAX) CALL KERNEL(T,X,N,B2,KORD,KORD+2,NYG,S,W1(1,1),M1,W1(1,3)) !- !-------- 14. ESTIMATE INTEGRALFUNCTIONAL IN ITERATIONS XMY2=.75*(W1(1,2)*W1(1,3)*W1(1,3)+W1(M1,2)*W1(M1,3)*W1(M1,3)) DO 140 I=2,M1-1 140 XMY2=XMY2+W1(I,2)*W1(I,3)*W1(I,3) XMY2=XMY2*(TU-TL)/DBLE(M1) !- !-------- 15. FINISH OF ITERATIONS B=(CONST/XMY2)**EX B=MAX(BMIN,B) B=MIN(BMAX,B) 120 CONTINUE !-------- 16 COMPUTE SMOOTHED FUNCTION WITH GLOBAL PLUG-IN BANDWIDTH 160 CALL KERNEL(T,X,N,B,NUE,KORD,NYG,S,TT,M,Y) !- !-------- 17. VARIANCE CHECK IF(IHOM.NE.0) SIG=RVAR IF(RVAR.EQ.SIG.OR.R2.LT..88.OR.IHOM.NE.0.OR.NUE.GT.0) GOTO 180 II=0 IIL=0 J=2 TLL=MAX(TL,TT(1)) TUU=MIN(TU,TT(M)) DO 170 I=IL,IU IF(T(I).LT.TLL.OR.T(I).GT.TUU) GOTO 170 II=II+1 IF(IIL.EQ.0) IIL=I 171 IF(TT(J).LT.T(I)) THEN J=J+1 IF(J.LE.M) GOTO 171 END IF WN(II,3)=X(I)-Y(J)+(Y(J)-Y(J-1))*(TT(J)-T(I))/(TT(J)-TT(J-1)) 170 CONTINUE CALL RESEST(T(IIL),WN(1,3),II,WN(1,4),SNR,RVAR) Q=SIG/RVAR CALL COFF(WN(1,4),N,SIG) IF(Q.LE.2.) GOTO 180 IF(Q.GT.5..AND.R2.GT..95) RVAR=RVAR*.5 SIG=RVAR CALL COFF(WN(1,4),N,SIG) GOTO 100 !- !-------- 18. LOCAL INITIALIZATIONS 180 BVAR=B NUEV=0 KORDV=2 NYL=1 !- !-------- 19. COMPUTE INNER BANDWIDTHS G1=0.86*(1.+DBLE(KORD-NUE-2)*.05)*B G1=G1*DBLE(N)**(4./DBLE(2*KORD+1)/(2*KORD+5)) G1=MAX(G1,BMIN/DBLE(KORD-1)*DBLE(KORD+1)) G1=MIN(G1,BMAX) G2=1.4*(1.+DBLE(KORD-NUE-2)*0.05)*B G2=G2*DBLE(N)**(2./DBLE(2*KORD+1)/DBLE(2*KORD+3)) G2=MAX(G2,BMIN) G2=MIN(G2,BMAX) !- !-------- 20. ESTIMATE/COMPUTE INTEGRAL CONSTANT VI LOCALLY DO 200 I=1,N 200 WN(I,4)=DBLE(N)*WN(I,4)*(S(I)-S(I-1)) DO 201 J=1,M BAN(J)=BVAR WM(j)=tt(j) IF(TT(J).LT.S(0)+G1) THEN DIST=((TT(J)-G1-S(0))/G1)**2 BAN(J)=BVAR*(1.0+1.0*DIST) BAN(J)=MIN(BAN(J),BMAX) WM(J)=TT(J)+.5*DIST*G1 ELSE IF(TT(J).GT.S(N)-G1) THEN DIST=((TT(J)-S(N)+G1)/G1)**2 BAN(J)=BVAR*(1.0+1.0*DIST) BAN(J)=MIN(BAN(J),BMAX) WM(J)=TT(J)-.5*DIST*G1 END IF 201 CONTINUE CALL KERNEL(T,WN(1,4),N,BVAR,NUEV,KORDV,NYL,S,WM,M,BAN) !- !-------- 21. ESTIMATION OF KORDTH DERIVATIVE LOCALLY WSTEP=(TT(M)-TT(1))/DBLE(M1-2) DO 210 J=2,M1 W1(J,2)=TT(1)+DBLE(J-2)*WSTEP 210 W1(J,1)=TT(1)+DBLE(J-1.5)*WSTEP W1(1,1)=TT(1)+.5*WSTEP CALL KERNEL(T,X,N,G1,KORD,KORD+2,NYG,S,W1(2,2),M1-1,W1(2,3)) DO 211 J=2,M1 211 W1(J,3)=W1(J,3)*W1(J,3) DO 212 J=1,M Y(J)=G2 IF(TT(J).LT.S(0)+G1) THEN Y(J)=G2*(1.0+1.0*((TT(J)-G1-S(0))/G1)**2) Y(J)=MIN(Y(J),BMAX) ELSE IF(TT(J).GT.S(N)-G1) THEN Y(J)=G2*(1.0+1.0*((TT(J)-S(N)+G1)/G1)**2) Y(J)=MIN(Y(J),BMAX) END IF 212 CONTINUE CALL KERNP(W1(2,2),W1(2,3),M1-1,G2,NUEV,KORDV,NYL, & W1(1,1),WM,M,Y) !- !-------- 22. FINISH DO 220 J=1,M XH=BMIN**(2*KORD+1)*ABS(Y(J))*VI/CONST XXH=CONST*ABS(BAN(J))/VI/BMAX**(2*KORD+1) IF(BAN(J).LT.XH) THEN BAN(J)=BMIN ELSE IF(Y(J).LT.XXH) THEN BAN(J)=BMAX ELSE BAN(J)=(CONST*BAN(J)/Y(J)/VI)**EX END IF END IF 220 CONTINUE !- !-------- 23. COMPUTE SMOOTHED FUNCTION WITH LOCAL PLUG-IN BANDWIDTH 230 DO 231 J=1,M 231 Y(J)=BAN(J) CALL KERNEL(T,X,N,B,NUE,KORD,NYL,S,TT,M,Y) RETURN END !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! SEVERAL KERNEL SMOOTHING SUBROUTINES WHICH ARE USED BY GLKERN.F ! AND LOKERN.F, VERSION JANUARY 1997 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! THIS FILE CONTAINS: ! ! SUBROUTINE RESEST(T,X,N,RES,SNR,SIGMA2) ! FOR VARIANCE ESTIMATION ! ! SUBROUTINE KERNEL(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) ! DRIVER SUBROUTINE FOR KERNEL REGRESSION ESTIMATION ! CALLS FAST OR CONVENTIAL KERNEL ROUTINE ! ! SUBROUTINE KERNP(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) ! DRIVER SUBROUTINE FOR KERNEL REGRESSION ESTIMATION ! WITHOUT USE OF BOUNDARY KERNELS ! --------------------------------------------------------------------- ! SUBROUTINE KERNFA(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) ! FAST ALGORITHM FOR KERNEL ESTIMATION ! SUBROUTINE DREG(SW,A1,A2,IORD,X,SL,SR,T,B,IFLOP) ! USED BY SUBROUTINE KERNFA,KERNFP ! SUBROUTINE LREG(SW,A3,IORD,D,DOLD,Q,C) ! USED BY SUBROUTINE KERNFA,KERNFP ! SUBROUTINE FREG(SW,NUE,KORD,IBOUN,Y,C,ICALL,A) ! USED BY SUBROUTINE KERNFA,KERNFP ! SUBROUTINE KERNFP(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) ! FAST ALGORITHM FOR KERNP ESTIMATION WITHOUT BOUNDARY ! --------------------------------------------------------------------- ! SUBROUTINE KERNCL(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) ! CONVENTIONAL ALGORITHM FOR KERNEL ESTIMATION ! SUBROUTINE SMO(S,X,N,TAU,WID,NUE,IORD,IBOUN,IST,S1,C,Y) ! SINGLE ESTIMATION STEP, USED BY KERNCL ! SUBROUTINE KERNCP(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) ! CONVENTIONAL ALGORITHM WITHOUT BOUNDARY KERNELS ! SUBROUTINE SMOP(S,X,N,TAU,WID,NUE,IORD,IBOUN,IST,S1,C,Y) ! SINGLE ESTIMATION STEP, USED BY KERNCP ! --------------------------------------------------------------------- ! SUBROUTINE COFFI(NUE,KORD,C) ! KERNEL COEFFICIENT OF POLYNOMIAL KERNELS USED BY ! KERNCL,KERNCP AND KERNFP ! SUBROUTINE COFFB(NUE,KORD,Q,IBOUN,C) ! KERNEL COEFFICIENT OF POLYNOMIAL BOUNDARY KERNELS ! USED BY KERNFA AND KERNCL ! --------------------------------------------------------------------- ! SUBROUTINE COFF(X,N,FA) ! SIMPLE SUBROUTINE FOR ARRAY INITIALIZATION !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE RESEST(T,X,N,RES,SNR,SIGMA2) !----------------------------------------------------------------------- ! VERSION: JUNE, 1996 ! ! PURPOSE: ! ! COMPUTES ONE-LEAVE-OUT RESIDUALS FOR NONPARAMETRIC ESTIMATION ! OF RESIDUAL VARIANCE (LOCAL LINEAR APPROXIMATION FOLLOWED BY ! REWEIGHTING) ! ! PARAMETERS: ! ! INPUT T(N) ABSCISSAE (ORDERED: T(I)<=T(I+1)) ! INPUT X(N) DATA ! INPUT N LENGTH OF DATA ( >2 ) ! OUTPUT RES(N) RESIDUALS AT T(1),...,T(N) ! OUTPUT SNR EXPLAINED VARIANCE OF THE TRUE CURVE ! OUTPUT SIGMA2 ESTIMATION OF SIGMA**2 (RESIDUAL VARIANCE) ! !----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N),T(N),RES(N) !- SIGMA2=0. EX=X(1)*(T(2)-T(1)) EX2=X(1)*EX DO 1 I=2,N-1 TT=T(I+1)-T(I-1) IF(TT.NE.0.) G1=(T(I+1)-T(I))/TT IF(TT.EQ.0.) G1=.5 G2=1.-G1 RES(I)=(X(I)-G1*X(I-1)-G2*X(I+1))/SQRT(1.+G1*G1+G2*G2) SIGMA2=SIGMA2+RES(I)*RES(I) SX=X(I)*TT EX=EX+SX EX2=EX2+X(I)*SX 1 CONTINUE TT=T(3)-T(2) IF(TT.NE.0.) G1=(T(1)-T(2))/TT IF(TT.EQ.0.) G1=.5 G2=1.-G1 RES(1)=(X(1)-G1*X(3)-G2*X(2))/SQRT(1.+G1*G1+G2*G2) TT=T(N-1)-T(N-2) IF(TT.NE.0.) G1=(T(N-1)-T(N))/TT IF(TT.EQ.0.) G1=.5 G2=1.-G1 RES(N)=(X(N)-G1*X(N-2)-G2*X(N-1))/SQRT(1.+G1*G1+G2*G2) SIGMA2=(SIGMA2+RES(1)*RES(1)+RES(N)*RES(N))/N !- SX=X(N)*(T(N)-T(N-1)) DN=2.*(T(N)-T(1)) EX=(EX+SX)/DN EX2=(EX2+X(N)*SX)/DN IF(EX2.EQ.0) SNR=0. IF(EX2.GT.0) SNR=1-SIGMA2/(EX2-EX*EX) RETURN END SUBROUTINE KERNEL(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) !----------------------------------------------------------------------- ! SHORT-VERSION MAY, 1995 ! ! DRIVER SUBROUTINE FOR KERNEL SMOOTHING, CHOOSES BETWEEN ! STANDARD AND O(N) ALGORITHM ! ! PARAMETERS : ! ! INPUT T(N) INPUT GRID (REGRESSION DESIGN) ! INPUT X(N) DATA, GIVEN ON T(N) ! INPUT N LENGTH OF X ! INPUT B ONE SIDED BANDWIDTH (FOR NY=1 MEAN BANDWIDTH) ! INPUT NUE ORDER OF DERIVATIVE (0-4) ! INPUT KORD ORDER OF KERNEL (<=6); DEFAULT IS KORD=NUE+2 ! INPUT NY 0: GLOBAL BANDWIDTH (DEFAULT) ! 1: VARIABLE BANDWIDTHS, GIVEN IN Y AS INPUT ! INPUT S(0:N) INTERPOLATION SEQUENCE ! INPUT TT(M) OUTPUT GRID. MUST BE PART OF INPUT GRID FOR ! IEQ=0 ! INPUT M NUMBER OF POINTS WHERE FUNCTION IS ESTIMATED, ! OR LENGTH OF TT. DEFAULT IS M=400 ! INPUT Y(M) BANDWITH SEQUENCE FOR NY=1, DUMMY FOR NY=0 ! OUTPUT Y(M) ESTIMATED REGRESSION FUNCTION ! !----------------------------------------------------------------------- DOUBLE PRECISION T(N),X(N),B,S(0:N),TT(M),Y(M),CHAN INTEGER N,NUE,KORD,NY,M !- !------ COMPUTING CHANGE POINT CHAN=(5.+KORD)*MAX(1.,SQRT(FLOAT(N)/FLOAT(M))) !------ IF(B*(N-1)/(T(N)-T(1)).LT.CHAN) THEN CALL KERNCL(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) ELSE CALL KERNFA(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) END IF !- RETURN END SUBROUTINE KERNP(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) !----------------------------------------------------------------------- ! SHORT-VERSION JANUARY, 1997 ! ! DRIVER SUBROUTINE FOR KERNEL SMOOTHING, CHOOSES BETWEEN ! STANDARD AND O(N) ALGORITHM WITHOUT USING BOUNDARY KERNELS ! ! PARAMETERS : ! ! INPUT T(N) INPUT GRID (REGRESSION DESIGN) ! INPUT X(N) DATA, GIVEN ON T(N) ! INPUT N LENGTH OF X ! INPUT B ONE SIDED BANDWIDTH (FOR NY=1 MEAN BANDWIDTH) ! INPUT NUE ORDER OF DERIVATIVE (0-4) ! INPUT KORD ORDER OF KERNEL (<=6); DEFAULT IS KORD=NUE+2 ! INPUT NY 0: GLOBAL BANDWIDTH (DEFAULT) ! 1: VARIABLE BANDWIDTHS, GIVEN IN Y AS INPUT ! INPUT S(0:N) INTERPOLATION SEQUENCE ! INPUT TT(M) OUTPUT GRID. MUST BE PART OF INPUT GRID FOR ! IEQ=0 ! INPUT M NUMBER OF POINTS WHERE FUNCTION IS ESTIMATED, ! OR LENGTH OF TT. DEFAULT IS M=400 ! INPUT Y(M) BANDWITH SEQUENCE FOR NY=1, DUMMY FOR NY=0 ! OUTPUT Y(M) ESTIMATED REGRESSION FUNCTION ! !----------------------------------------------------------------------- DOUBLE PRECISION T(N),X(N),B,S(0:N),TT(M),Y(M),CHAN INTEGER N,NUE,KORD,NY,M !- !------ COMPUTING CHANGE POINT CHAN=(5.+KORD)*MAX(1.,SQRT(FLOAT(N)/FLOAT(M))) !------ IF(B*(N-1)/(T(N)-T(1)).LT.CHAN) THEN CALL KERNCP(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) ELSE CALL KERNFP(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) END IF !- RETURN END SUBROUTINE KERNFA(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) !----------------------------------------------------------------------- ! SHORT-VERSION: MAY, 1995 ! ! PURPOSE: ! ! COMPUTATION OF KERNEL ESTIMATE USING O(N) ALGORITHM BASED ON ! LEGENDRE POLYNOMIALS, GENERAL SPACED DESIGN AND LOCAL ! BANDWIDTH ALLOWED. (NEW INITIALISATIONS OF THE LEGENDRE SUMS ! FOR NUMERICAL REASONS) ! ! PARAMETERS : ! ! INPUT T(N) INPUT GRID ! INPUT X(N) DATA, GIVEN ON T(N) ! INPUT N LENGTH OF X ! INPUT B ONE SIDED BANDWIDTH ! INPUT NUE ORDER OF DERIVATIVE (0-4) ! INPUT KORD ORDER OF KERNEL (<=6) ! INPUT NY 0, GLOBAL BANDWIDTH; 1, LOCAL BANDWIDTH IN Y ! INPUT S(0:N) HALF POINT INTERPOLATION SEQUENCE ! INPUT TT(M) OUTPUT GRID ! INPUT M NUMBER OF POINTS TO ESTIMATE ! INPUT Y(M) BANDWITH SEQUENCE FOR NY=1, DUMMY FOR NY=0 ! OUTPUT Y(M) ESTIMATED FUNCTION ! ! !----------------------------------------------------------------------- INTEGER N,NUE,KORD,NY,M INTEGER J,K,IORD,INIT,ICALL,I,IBOUN INTEGER JL,JR,JNR,JNL DOUBLE PRECISION X(N),T(N),S(0:N),TT(M),Y(M),B DOUBLE PRECISION C(7),SW(7),XF(7),DOLD DOUBLE PRECISION A(7,7),A1(7),A2(7),A3(7,7),CM(7,6) DOUBLE PRECISION BMIN,BMAX,BB,WWL,WWR,WID,WR,WIDO !- !------ COMPUTE CONSTANTS FOR LATER USE S0=1.5*T(1)-0.5*T(2) SN=1.5*T(N)-0.5*T(N-1) BMIN=(SN-S0)*.6D0/DBLE(N)*DBLE(KORD-1) BMAX=(S(N)-S(0))*.5 IF(KORD.EQ.2) BMIN=BMIN*0.1D0 IORD=KORD+1 DO 2 K=3,IORD A1(K)=DBLE(2*K-1)/DBLE(K) 2 A2(K)=DBLE(1-K)/DBLE(K) !- INIT=0 ICALL=0 DOLD=0.D0 !- !------ SMOOTHING LOOP DO 100 I=1,M BB=B IF (NY .EQ. 1) BB=Y(I) IF(BB.LT.BMIN) BB=BMIN IF(BB.GT.BMAX) BB=BMAX IBOUN=0 !- !------ COMPUTE LEFT BOUNDARY KERNEL IF(TT(I).LT.S(0)+BB) THEN WWL=S(0) WWR=S(0)+BB+BB WID=WWR-TT(I) IBOUN=1 CALL COFFB(NUE,KORD,(TT(I)-S(0))/WID,IBOUN,C) END IF !- !------ COMPUTE RIGHT BOUNDARY KERNEL IF(TT(I)+BB.GT.S(N)) THEN WWL=S(N)-(BB+BB) WWR=S(N) WID=TT(I)-WWL IBOUN=-1 CALL COFFB(NUE,KORD,(S(N)-TT(I))/WID,IBOUN,C) END IF !- !------ NO BOUNDARY IF(IBOUN.EQ.0) THEN WID=BB WWL=TT(I)-BB WWR=TT(I)+BB END IF !- !------ INITIALISATION FOR INIT=0 IF(INIT.EQ.0) THEN DO 44 K=1,IORD 44 SW(K)=0. JL=1 DO 48 J=1,N IF(S(J-1).LT.WWL) THEN JL=J+1 ELSE IF(S(J).GT.WWR) GOTO 488 CALL DREG(SW,A1,A2,IORD,X(J),S(J-1),S(J),TT(I),WID,1) END IF 48 CONTINUE 488 JR=J-1 WR=WWR INIT=1 GOTO 6666 ELSE INIT=INIT+1 END IF !- !------ COMPARE OLD SUM WITH NEW SMOOTHING INTERVALL TT(I)-B,TT(I)+B IF(S(JR-1).GE.WWL) THEN JNR=JR JNL=JL IF(S(JR).GT.WWR) THEN DO 201 J=JR,JL,-1 CALL DREG(SW,A1,A2,IORD,X(J),S(J-1),S(J),TT(I-1),WIDO,-1) JNR=J-1 IF(S(JNR).LE.WWR) GOTO 2011 201 CONTINUE 2011 END IF IF(S(JL-1).LT.WWL) THEN DO 301 J=JL,JR CALL DREG(SW,A1,A2,IORD,X(J),S(J-1),S(J),TT(I-1),WIDO,-1) JNL=J+1 IF(S(J).GE.WWL) GOTO 3011 301 CONTINUE 3011 END IF !- !------ UPDATING OF SW CALL LREG(SW,A3,IORD,(TT(I)-TT(I-1))/WID,DOLD,WIDO/WID,CM) IF(JNR.EQ.JR) THEN DO 401 J=JR+1,N IF(S(J).GT.WWR) GOTO 4011 CALL DREG(SW,A1,A2,IORD,X(J),S(J-1),S(J),TT(I),WID,1) JNR=J 401 CONTINUE 4011 END IF JR=JNR IF(JL.EQ.JNL) THEN DO 402 J=JL-1,1,-1 IF(S(J-1).LT.WWL) GOTO 4022 CALL DREG(SW,A1,A2,IORD,X(J),S(J-1),S(J),TT(I),WID,1) JNL=J 402 CONTINUE 4022 END IF JL=JNL ELSE !- !------ NEW INITIALISATION OF SW DO 22 K=1,IORD 22 SW(K)=0. DO 202 J=JR,N IF(S(J-1).LT.WWL) THEN JL=J+1 ELSE IF(S(J).GT.WWR) GOTO 2022 CALL DREG(SW,A1,A2,IORD,X(J),S(J-1),S(J),TT(I),WID,1) END IF 202 CONTINUE 2022 JR=J-1 WR=WWR END IF 6666 CONTINUE !- !------ IF BANDWIDTH IS TOO SMALL NO SMOOTHING IF(WWL.GE.S(JR-1).AND.WWR.LE.S(JR)) THEN Y(I)=X(JR) IF(NUE.GT.0) Y(I)=0.D0 ELSE !- !------ ADD FIRST AND LAST POINT OF THE SMOOTHING INTERVAL DO 501 K=1,IORD 501 XF(K)=SW(K) IF(JL.NE.1) & CALL DREG(XF,A1,A2,IORD,X(JL-1),WWL,S(JL-1),TT(I),WID,1) IF(JR.NE.N) & CALL DREG(XF,A1,A2,IORD,X(JR+1),S(JR),WWR,TT(I),WID,1) !- !------ NOW THE SUMS ARE BUILT THAT ARE NEEDED TO COMPUTE THE ESTIMATE CALL FREG(XF,NUE,KORD,IBOUN,Y(I),C,ICALL,A) IF(NUE.GT.0) Y(I)=Y(I)/(WID**NUE) END IF !- !------ NEW INITIALISATION ? IF(JL.GT.JR.OR.WWL.GT.WR.OR.INIT.GT.100) INIT=0 WIDO=WID !- 100 CONTINUE !- RETURN END SUBROUTINE KERNFP(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) !----------------------------------------------------------------------- ! SHORT-VERSION: JANUARY, 1997 ! ! PURPOSE: ! ! COMPUTATION OF KERNEL ESTIMATE USING O(N) ALGORITHM BASED ON ! LEGENDRE POLYNOMIALS, GENERAL SPACED DESIGN AND LOCAL ! BANDWIDTH ALLOWED. (NEW INITIALISATIONS OF THE LEGENDRE SUMS ! FOR NUMERICAL REASONS) WITHOUT BOUNDARY KERNELS, JUST NORMALIZING ! ! PARAMETERS : ! ! INPUT T(N) INPUT GRID ! INPUT X(N) DATA, GIVEN ON T(N) ! INPUT N LENGTH OF X ! INPUT B ONE SIDED BANDWIDTH ! INPUT NUE ORDER OF DERIVATIVE (0-4) ! INPUT KORD ORDER OF KERNEL (<=6) ! INPUT NY 0, GLOBAL BANDWIDTH; 1, LOCAL BANDWIDTH IN Y ! INPUT S(0:N) HALF POINT INTERPOLATION SEQUENCE ! INPUT TT(M) OUTPUT GRID ! INPUT M NUMBER OF POINTS TO ESTIMATE ! INPUT Y(M) BANDWITH SEQUENCE FOR NY=1, DUMMY FOR NY=0 ! OUTPUT Y(M) ESTIMATED FUNCTION ! ! !----------------------------------------------------------------------- INTEGER N,NUE,KORD,NY,M INTEGER J,K,IORD,INIT,ICALL,I,IBOUN INTEGER JL,JR,JNR,JNL DOUBLE PRECISION X(N),T(N),S(0:N),TT(M),Y(M),B DOUBLE PRECISION C(7),SW(7),XF(7),DOLD,QQ,Q,XNOR DOUBLE PRECISION A(7,7),A1(7),A2(7),A3(7,7),CM(7,6) DOUBLE PRECISION BMIN,BMAX,BB,WWL,WWR,WID,WR,WIDO !- !------ COMPUTE CONSTANTS FOR LATER USE S0=1.5*T(1)-0.5*T(2) SN=1.5*T(N)-0.5*T(N-1) BMIN=(SN-S0)*.6D0/DBLE(N)*DBLE(KORD-1) BMAX=(S(N)-S(0))*.5 IF(KORD.EQ.2) BMIN=BMIN*0.1D0 IORD=KORD+1 CALL COFFI(NUE,KORD,C) DO 2 K=3,IORD A1(K)=DBLE(2*K-1)/DBLE(K) 2 A2(K)=DBLE(1-K)/DBLE(K) !- INIT=0 ICALL=0 DOLD=0.D0 !- !------ SMOOTHING LOOP DO 100 I=1,M BB=B IF (NY .EQ. 1) BB=Y(I) IF(BB.LT.BMIN) BB=BMIN IF(BB.GT.BMAX) BB=BMAX IBOUN=0 !- !------ COMPUTE LEFT BOUNDARY IF(TT(I).LT.S(0)+BB) THEN WWL=S(0) WWR=S(0)+BB+BB WID=WWR-TT(I) IBOUN=1 END IF !- !------ COMPUTE RIGHT BOUNDARY IF(TT(I)+BB.GT.S(N)) THEN WWL=S(N)-(BB+BB) WWR=S(N) WID=TT(I)-WWL IBOUN=-1 END IF !- !------ NO BOUNDARY IF(IBOUN.EQ.0) THEN WID=BB WWL=TT(I)-BB WWR=TT(I)+BB XNOR=1.D0 END IF !- !------ COMPUTE NORMALIZING CONSTANT IF(IBOUN.NE.0) THEN IF(IBOUN.EQ.1) Q=(TT(I)-S(0))/WID IF(IBOUN.EQ.-1) Q=(S(N)-TT(I))/WID QQ=Q*Q XNOR=C(1)*(1.D0+Q) DO 3 K=3,IORD,2 Q=Q*QQ 3 XNOR=XNOR+C(K)*(1.D0+Q) IBOUN=0 END IF !- !------ INITIALISATION FOR INIT=0 IF(INIT.EQ.0) THEN DO 44 K=1,IORD 44 SW(K)=0. JL=1 DO 48 J=1,N IF(S(J-1).LT.WWL) THEN JL=J+1 ELSE IF(S(J).GT.WWR) GOTO 488 CALL DREG(SW,A1,A2,IORD,X(J),S(J-1),S(J),TT(I),WID,1) END IF 48 CONTINUE 488 JR=J-1 WR=WWR INIT=1 GOTO 6666 ELSE INIT=INIT+1 END IF !- !------ COMPARE OLD SUM WITH NEW SMOOTHING INTERVALL TT(I)-B,TT(I)+B IF(S(JR-1).GE.WWL) THEN JNR=JR JNL=JL IF(S(JR).GT.WWR) THEN DO 201 J=JR,JL,-1 CALL DREG(SW,A1,A2,IORD,X(J),S(J-1),S(J),TT(I-1),WIDO,-1) JNR=J-1 IF(S(JNR).LE.WWR) GOTO 2011 201 CONTINUE 2011 END IF IF(S(JL-1).LT.WWL) THEN DO 301 J=JL,JR CALL DREG(SW,A1,A2,IORD,X(J),S(J-1),S(J),TT(I-1),WIDO,-1) JNL=J+1 IF(S(J).GE.WWL) GOTO 3011 301 CONTINUE 3011 END IF !- !------ UPDATING OF SW CALL LREG(SW,A3,IORD,(TT(I)-TT(I-1))/WID,DOLD,WIDO/WID,CM) IF(JNR.EQ.JR) THEN DO 401 J=JR+1,N IF(S(J).GT.WWR) GOTO 4011 CALL DREG(SW,A1,A2,IORD,X(J),S(J-1),S(J),TT(I),WID,1) JNR=J 401 CONTINUE 4011 END IF JR=JNR IF(JL.EQ.JNL) THEN DO 402 J=JL-1,1,-1 IF(S(J-1).LT.WWL) GOTO 4022 CALL DREG(SW,A1,A2,IORD,X(J),S(J-1),S(J),TT(I),WID,1) JNL=J 402 CONTINUE 4022 END IF JL=JNL ELSE !- !------ NEW INITIALISATION OF SW DO 22 K=1,IORD 22 SW(K)=0. DO 202 J=JR,N IF(S(J-1).LT.WWL) THEN JL=J+1 ELSE IF(S(J).GT.WWR) GOTO 2022 CALL DREG(SW,A1,A2,IORD,X(J),S(J-1),S(J),TT(I),WID,1) END IF 202 CONTINUE 2022 JR=J-1 WR=WWR END IF 6666 CONTINUE !- !------ IF BANDWIDTH IS TOO SMALL NO SMOOTHING IF(WWL.GE.S(JR-1).AND.WWR.LE.S(JR)) THEN Y(I)=X(JR) IF(NUE.GT.0) Y(I)=0.D0 ELSE !- !------ ADD FIRST AND LAST POINT OF THE SMOOTHING INTERVAL DO 501 K=1,IORD 501 XF(K)=SW(K) IF(JL.NE.1) & CALL DREG(XF,A1,A2,IORD,X(JL-1),WWL,S(JL-1),TT(I),WID,1) IF(JR.NE.N) & CALL DREG(XF,A1,A2,IORD,X(JR+1),S(JR),WWR,TT(I),WID,1) !- !------ NOW THE SUMS ARE BUILT THAT ARE NEEDED TO COMPUTE THE ESTIMATE CALL FREG(XF,NUE,KORD,IBOUN,Y(I),C,ICALL,A) IF(NUE.GT.0) Y(I)=Y(I)/(WID**NUE) IF(XNOR.NE.1.D0) Y(I)=Y(I)/XNOR END IF !- !------ NEW INITIALISATION ? IF(JL.GT.JR.OR.WWL.GT.WR.OR.INIT.GT.100) INIT=0 WIDO=WID !- 100 CONTINUE !- RETURN END SUBROUTINE DREG(SW,A1,A2,IORD,X,SL,SR,T,B,IFLOP) !----------------------------------------------------------------------- ! VERSION: MAY, 1995 ! ! PURPOSE: ! ! COMPUTES NEW LEGENDRE SUMS (FOR REGRESSION) ! ! PARAMETERS: ! ************** INPUT ******************* ! ! SW(IORD) : OLD SUM OF DATA WEIGHTS FOR LEGENDRE POLYNOM. ! A1(7) : CONSTANTS OF RECURSIVE FORMULA FOR LEGENDRE POL. ! A2(7) : " ! IORD : ORDER OF KERNEL POLYNOMIAL ! X : DATA POINT ! SL : LEFT S-VALUE ! SR : RIGHT S-VALUE ! T : POINT WHERE THE SMOOTHED VALUE IS TO BE ESTIMATED ! B : BANDWIDTH ! IFLOP : 1: ADDITION, ELSE SUBTRACTION ! ! ! ************** OUTPUT ******************* ! ! SW(IORD) : NEW SUM OF DATA WEIGHTS FOR LEGENDRE POLYNOM. ! !----------------------------------------------------------------------- INTEGER IORD,IFLOP,K DOUBLE PRECISION SW(7),A1(7),A2(7),X,SL,SR,T,B DOUBLE PRECISION P(7,2) !- !------ COMPUTE LEGENDRE POLYNOMIALS P(1,1)=(T-SL)/B P(1,2)=(T-SR)/B P(2,1)=1.5D0*P(1,1)*P(1,1)-.5D0 P(2,2)=1.5D0*P(1,2)*P(1,2)-.5D0 DO 1 K=3,IORD P(K,1)=A1(K)*P(K-1,1)*P(1,1)+A2(K)*P(K-2,1) 1 P(K,2)=A1(K)*P(K-1,2)*P(1,2)+A2(K)*P(K-2,2) !- !------ COMPUTE NEW LEGENDRE SUMS IF(IFLOP.EQ.1) THEN DO 2 K=1,IORD 2 SW(K)=SW(K)+(P(K,1)-P(K,2))*X ELSE DO 3 K=1,IORD 3 SW(K)=SW(K)+(P(K,2)-P(K,1))*X END IF RETURN END SUBROUTINE LREG(SW,A3,IORD,D,DOLD,Q,C) !------------------------------------------------------------------ ! VERSION: MAY, 1995 ! ! PURPOSE: ! ! UPDATE OF SW-SEQUENCE ACCORDING TO NEW BANDWIDTH AND NEW DATA ! (VERSION FOR REGRESSION) ! ! PARAMETERS : ! ************** INPUT ******************* ! ! SW(IORD) : SUM OF DATA WEIGHTS FOR LEGENDRE POLYNOM. ! IORD : ORDER OF KERNEL POLYNOMIAL ! D : DIST. TO THE NEXT POINT DIVIDED BY BANDW. ! DOLD : D PREVIOUS STEP ! Q : NEW BANDWIDTH DIVIDED BY OLD BANDWIDTH ! ************** WORK ******************* ! ! A3(7,7) : MATRIX (P*Q*P)**(-1) ! C(7,6) : MATRIX OF COEFFICIENTS ! ************** OUTPUT ****************** ! ! SW(IORD) : UPDATED VERSION OF SW ! !--------------------------------------------------------------------- INTEGER IORD,K,I,L DOUBLE PRECISION D,DOLD,Q,DD,WW,QQ,XX DOUBLE PRECISION A3(7,7),C(7,6),SW(7) !- !- BUILD UP MATRIX IF(DOLD.NE.D.OR.DOLD.EQ.0) THEN DOLD=D DD=D*D !- IF(IORD.EQ.7) THEN C(7,6)=13.D0*D C(7,5)=71.5D0*DD C(7,4)=(214.5D0*DD+9.D0)*D C(7,3)=(375.375D0*DD+77.D0)*DD C(7,2)=((375.375D0*DD+247.5D0)*DD+5.D0)*D C(7,1)=((187.6875D0*DD+346.5D0)*DD+40.5D0)*DD END IF !- IF(IORD.GE.6) THEN C(6,5)=11.D0*D C(6,4)=49.5D0*DD C(6,3)=(115.5D0*DD+7.D0)*D C(6,2)=(144.375D0*DD+45.D0)*DD C(6,1)=((86.625D0*DD+94.5D0)*DD+3.D0)*D END IF !- IF(IORD.GE.5) THEN C(5,4)=9.D0*D C(5,3)=31.5D0*DD C(5,2)=(52.5D0*DD+5.D0)*D C(5,1)=(39.375D0*DD+21.D0)*DD END IF !- IF(IORD.GE.4) THEN C(4,3)=7.D0*D C(4,2)=17.5D0*DD C(4,1)=(17.5D0*DD+3.D0)*D END IF !- IF(IORD.GE.3) THEN C(3,2)=5.D0*D C(3,1)=7.5D0*DD END IF !- C(2,1)=3.D0*D END IF IF(Q.LT..9999.OR.Q.GT.1.0001) THEN !- !------- BUILT UP MATRIX A3=P*Q*P**-1 A3(1,1)=Q DO 1 K=2,IORD 1 A3(K,K)=A3(K-1,K-1)*Q WW=Q*Q-1.D0 DO 2 K=1,IORD-2 WW=WW*Q 2 A3(K+2,K)=(K+.5D0)*WW !- IF(IORD.GE.5) THEN QQ=A3(2,2) A3(5,1)=Q*(1.875D0+QQ*(-5.25D0+QQ*3.375D0)) END IF IF(IORD.GE.6) A3(6,2)=QQ*(4.375D0+QQ*(-11.25D0+QQ*6.875D0)) IF(IORD.EQ.7) THEN A3(7,1)=Q*(-2.1875D0+QQ*(11.8125D0+QQ*(-18.5625D0+QQ* & 8.9375D0))) A3(7,3)=Q*QQ*(7.875D0+QQ*(-19.25D0+QQ*11.375D0)) END IF !- !------- COMPUTE A*C AND NEW LEGENDRE SUMS DO 10 I=IORD,2,-1 XX=0. DO 20 K=1,I WW=0. DO 30 L=K,I-1,2 30 WW=WW+A3(L,K)*C(I,L) IF(MOD(I-K,2).EQ.0) WW=WW+A3(I,K) XX=XX+WW*SW(K) 20 CONTINUE SW(I)=XX 10 CONTINUE SW(1)=A3(1,1)*SW(1) ELSE DO 111 I=IORD,2,-1 DO 112 K=1,I-1 112 SW(I)=SW(I)+C(I,K)*SW(K) 111 CONTINUE END IF RETURN END SUBROUTINE FREG(SW,NUE,KORD,IBOUN,Y,C,ICALL,A) !------------------------------------------------------------------ ! SHORT-VERSION: MAY, 1995 ! ! PURPOSE: ! ! FINAL COMPUTATION OF A SMOOTHED VALUE VIA LEGENDRE POLYNOMIALS ! ! PARAMETERS : ! ************** INPUT ******************* ! ! SW(KORD+1): SUM OF DATA WEIGHTS FOR LEGENDRE POLYNOM. ! NUE : ORDER OF DERIVATIVE ! KORD : ORDER OF KERNEL ! IBOUN : 0: INTERIOR KERNEL, ELSE BOUNDARY KERNEL ! C(KORD+1) : SEQUENCE OF POLYN. COEFF. FOR BOUND. KERNEL ! ICALL : PARAMETER USED TO INITIALISE COMPUTATION ! : OF A MATRIX ! ************** WORK ****************** ! ! A(7,7) : MATRIX OF COEFFICIENTS ! ************** OUTPUT ****************** ! ! Y : COMPUTED ESTIMATE ! !-------------------------------------------------------------------- INTEGER NUE,KORD,IBOUN,ICALL,I,J DOUBLE PRECISION SW(7),C(7),A(7,7),Y,WW !- !------- DEFINITION OF LEGENDRE COEFFICIENTS FOR BOUNDARY IF(ICALL.EQ.0.AND.IBOUN.NE.0) THEN A(2,2)=2./3. A(1,3)=.6 A(3,3)=.4 A(2,4)=4./7. A(4,4)=8./35. A(1,5)=27./63. A(3,5)=28./63. A(5,5)=8./63. A(2,6)=110./231. A(4,6)=72./231. A(6,6)=16./231. A(1,7)=143./429. A(3,7)=182./429. A(5,7)=88./429. A(7,7)=16./429. ICALL=1 END IF IF(IBOUN.NE.0) THEN !- !------- COMPUTATION OF THE SMOOTHED VALUE AT BOUNDARY Y=C(1)*SW(1)+C(2)*A(2,2)*SW(2) DO 1 J=3,KORD+1 WW=A(J,J)*SW(J) DO 2 I=J-2,1,-2 2 WW=WW+A(I,J)*SW(I) Y=Y+C(J)*WW 1 CONTINUE ELSE !- !------- COMPUTATION OF THE SMOOTHED VALUE AT INTERIOR IF(NUE.EQ.0) THEN IF(KORD.EQ.2) Y=-.1*SW(3)+.6*SW(1) IF(KORD.EQ.4) Y=(SW(5)-4.*SW(3)+9.*SW(1))/12. IF(KORD.EQ.6) Y=-7.2115379E-02*SW(7)+.25961537*SW(5) & -.4375*SW(3)+.75*SW(1) END IF IF(NUE.EQ.1) THEN IF(KORD.EQ.3) Y=(3.*SW(4)-10.*SW(2))/14. IF(KORD.EQ.5) Y=(-15.*SW(6)+48.*SW(4)-55.*SW(2))/44. END IF IF(NUE.EQ.2) THEN IF(KORD.EQ.4) Y=(-5.*SW(5)+14.*SW(3)-9.*SW(1))/6. IF(KORD.EQ.6) Y=2.01923*SW(7)-5.76923*SW(5)+5.25*SW(3) & -1.5*SW(1) END IF IF(NUE.EQ.3) Y=4.772727*SW(6)-12.272727*SW(4)+7.5*SW(2) IF(NUE.EQ.4) Y=-36.34615*SW(7)+88.84615*SW(5)-52.5*SW(3) END IF !- RETURN END SUBROUTINE KERNCL(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) !----------------------------------------------------------------------- ! SHORT-VERSION JANUARY 1995 ! ! KERNEL SMOOTHING, CONVENTIONAL ALGORITHM,GENERAL ! DESIGN, LOCAL BANDWIDTH ALLOWED ! ! PARAMETERS : ! ! INPUT T(N) INPUT GRID ! INPUT X(N) DATA, GIVEN ON T(N) ! INPUT N LENGTH OF X ! INPUT B ONE SIDED BANDWIDTH (FOR NY=1 MEAN BANDWIDTH) ! INPUT NUE ORDER OF DERIVATIVE (0-4) ! INPUT KORD ORDER OF KERNEL (<=6) ! INPUT S(0:N) HALF POINT INTERPOLATION SEQUENCE ! INPUT NY 0: GLOBAL, 1: LOCAL BANDWIDTH INPUTED IN Y ! INPUT TT(M) OUTPUT GRID ! INPUT M NUMBER OF POINTS TO ESTIMATE ! INPUT Y(M) BANDWITH SEQUENCE FOR NY=1, DUMMY FOR NY=0 ! OUTPUT Y(M) ESTIMATED REGRESSION FUNCTION ! !----------------------------------------------------------------------- DOUBLE PRECISION X(N),T(N),S(0:N),TT(M),Y(M) DOUBLE PRECISION C(7),C1(7) INTEGER N,NUE,KORD,NY,M,IST,I,IBOUN,IORD DOUBLE PRECISION B, BB, BMAX, WID, S1, BMIN !- !------ COMPUTE KERNEL COEFFICIENTS FOR INTERIOR AND SOME CONSTANTS CALL COFFI(NUE,KORD,C) IORD=KORD+1 BB=B S0=1.5*T(1)-0.5*T(2) SN=1.5*T(N)-0.5*T(N-1) BMIN=(SN-S0)*.6D0/DBLE(N)*DBLE(KORD-1) BMAX=(S(N)-S(0))*.5 IF(KORD.EQ.2) BMIN=0.1D0*BMIN IST=1 !- !------- LOOP OVER OUTPUT GRID DO 100 I=1,M IF(NY.NE.0) BB=Y(I) IF(BB.GT.BMAX) BB=BMAX IF(BB.LT.BMIN) BB=BMIN WID=BB S1=TT(I)-BB IBOUN=0 !- !------- COMPUTE LEFT BOUNDARY KERNEL IF(S1.LT.S(0)) THEN S1=S(0) WID=BB+BB+S(0)-TT(I) CALL COFFB(NUE,KORD,(TT(I)-S(0))/WID,1,C1) IBOUN=1 END IF !- !------- COMPUTE RIGHT BOUNDARY KERNEL IF(TT(I)+BB.GT.S(N)) THEN S1=S(N)-(BB+BB) WID=TT(I)-S1 CALL COFFB(NUE,KORD,(S(N)-TT(I))/WID,-1,C1) IBOUN=-1 END IF !- !------ SEARCH FIRST S-POINT OF SMOOTHING INTERVAL 2 IF(S(IST).LE.S1) THEN IST=IST+1 GOTO 2 END IF 3 IF(S(IST-1).GT.S1) THEN IST=IST-1 GOTO 3 END IF !- !------- IF BANDWIDTH IS TOO SMALL NO SMOOTHING IF(S(IST).GE.TT(I)+WID.OR.IST.EQ.N) THEN Y(I)=X(IST) IF(NUE.GT.0) Y(I)=0. ELSE !- !----- COMPUTE SMOOTHED DATA AT TT(I) IF (IBOUN.NE.0) THEN CALL SMO(S,X,N,TT(I),WID,NUE,IORD,IBOUN,IST,S1,C1,Y(I)) ELSE CALL SMO(S,X,N,TT(I),WID,NUE,IORD,IBOUN,IST,S1,C,Y(I)) END IF END IF 100 CONTINUE !- RETURN END SUBROUTINE SMO(S,X,N,TAU,WID,NUE,IORD,IBOUN,IST,S1,C,Y) !----------------------------------------------------------------------- ! SHORT-VERSION JANUARY 1995 ! ! PERFORMS ONE SMOOTHING STEP ! ! PARAMETERS : ! ! INPUT S(0:N) HALF POINT INTERPOLATION SEQUENCE ! INPUT X(N) DATA ! INPUT N LENGTH OF X ! INPUT TAU POINT WHERE FUNCTION IS ESTIMATED ! INPUT WID ONE SIDED BANDWIDTH ! INPUT NUE ORDER OF DERIVATIVE (0-4) ! INPUT IORD ORDER OF KERNEL POLYNOMIAL ! INPUT IBOUN TYPE OF BOUNDARY ! INPUT IST INDEX OF FIRST POINT OF SMOOTHING INTERVAL ! INPUT S1 LEFT BOUNDARY OF SMOOTHING INTERVAL ! INPUT C(7) KERNEL COEFFICIENTS ! OUTPUT Y SMOOTHED VALUE AT TAU ! WORK WO(7) WORK ARRAY ! !----------------------------------------------------------------------- DOUBLE PRECISION X(N),S(0:N),WO(7) DOUBLE PRECISION C(7) INTEGER N,NUE,IORD,IBOUN,IST,JEND,IBEG,INCR,I,J DOUBLE PRECISION TAU,WID,S1,Y,YY,YYY,W,WIDNUE !- Y=0. JEND=0 IBEG=2 IF(IBOUN.NE.0.OR.(NUE.NE.1.AND.NUE.NE.3)) IBEG=1 INCR=2 IF(IBOUN.NE.0) INCR=1 !- !------ COMPUTE INITIAL KERNEL VALUES IF(IBOUN.GT.0) THEN YY=(TAU-S1)/WID WO(IBEG)=YY DO 1 I=IBEG,IORD-INCR,INCR 1 WO(I+INCR)=WO(I)*YY ELSE DO 2 I=IBEG,IORD,INCR 2 WO(I)=1. END IF !- !------ LOOP OVER SMOOTHING INTERVAL DO 100 J=IST,N YY=(TAU-S(J))/WID IF(YY.LT.-1.) THEN YY=-1. JEND=1 END IF YYY=YY IF(IBOUN.EQ.0) THEN YY=YY*YY IF(NUE.EQ.1.OR.NUE.EQ.3) YYY=YY END IF !- !------ LOOP FOR COMPUTING WEIGHTS W=0. DO 3 I=IBEG,IORD,INCR W=W+C(I)*(WO(I)-YYY) WO(I)=YYY YYY=YYY*YY 3 CONTINUE Y=Y+W*X(J) IF(JEND.EQ.1) GOTO 110 100 CONTINUE !- !------- NORMALIZING FOR NUE>0 110 IF(NUE.GT.0) THEN WIDNUE=WID**NUE Y=Y/WIDNUE END IF !- RETURN END SUBROUTINE KERNCP(T,X,N,B,NUE,KORD,NY,S,TT,M,Y) !----------------------------------------------------------------------- ! SHORT-VERSION JANUARY 1997 ! ! KERNEL SMOOTHING, CONVENTIONAL ALGORITHM,GENERAL ! DESIGN, LOCAL BANDWIDTH ALLOWED, WITHOUT BOUNDARY KERNELS, ! JUST NORMALIZING ! ! PARAMETERS : ! ! INPUT T(N) INPUT GRID ! INPUT X(N) DATA, GIVEN ON T(N) ! INPUT N LENGTH OF X ! INPUT B ONE SIDED BANDWIDTH (FOR NY=1 MEAN BANDWIDTH) ! INPUT NUE ORDER OF DERIVATIVE (0-4) ! INPUT KORD ORDER OF KERNEL (<=6) ! INPUT S(0:N) HALF POINT INTERPOLATION SEQUENCE ! INPUT NY 0: GLOBAL, 1: LOCAL BANDWIDTH INPUTED IN Y ! INPUT TT(M) OUTPUT GRID ! INPUT M NUMBER OF POINTS TO ESTIMATE ! INPUT Y(M) BANDWITH SEQUENCE FOR NY=1, DUMMY FOR NY=0 ! OUTPUT Y(M) ESTIMATED REGRESSION FUNCTION ! !----------------------------------------------------------------------- DOUBLE PRECISION X(N),T(N),S(0:N),TT(M),Y(M) DOUBLE PRECISION C(7),C1(7) INTEGER N,NUE,KORD,NY,M,IST,I,IBOUN,IORD DOUBLE PRECISION B, BB, BMAX, WID, S1, BMIN !- !------ COMPUTE KERNEL COEFFICIENTS FOR INTERIOR AND SOME CONSTANTS CALL COFFI(NUE,KORD,C) IORD=KORD+1 BB=B S0=1.5*T(1)-0.5*T(2) SN=1.5*T(N)-0.5*T(N-1) BMIN=(SN-S0)*.6D0/DBLE(N)*DBLE(KORD-1) BMAX=(S(N)-S(0))*.5 IF(KORD.EQ.2) BMIN=0.1D0*BMIN IST=1 !- !------- LOOP OVER OUTPUT GRID DO 100 I=1,M IF(NY.NE.0) BB=Y(I) IF(BB.GT.BMAX) BB=BMAX IF(BB.LT.BMIN) BB=BMIN WID=BB S1=TT(I)-BB IBOUN=0 !- !------- COMPUTE LEFT BOUNDARY KERNEL IF(S1.LT.S(0)) THEN S1=S(0) WID=BB+BB+S(0)-TT(I) CALL COFFB(NUE,KORD,(TT(I)-S(0))/WID,1,C1) IBOUN=1 END IF !- !------- COMPUTE RIGHT BOUNDARY KERNEL IF(TT(I)+BB.GT.S(N)) THEN S1=S(N)-(BB+BB) WID=TT(I)-S1 IBOUN=-1 END IF !- !------ SEARCH FIRST S-POINT OF SMOOTHING INTERVAL 2 IF(S(IST).LE.S1) THEN IST=IST+1 GOTO 2 END IF 3 IF(S(IST-1).GT.S1) THEN IST=IST-1 GOTO 3 END IF !- !------- IF BANDWIDTH IS TOO SMALL NO SMOOTHING IF(S(IST).GE.TT(I)+WID.OR.IST.EQ.N) THEN Y(I)=X(IST) IF(NUE.GT.0) Y(I)=0. ELSE !- !----- COMPUTE SMOOTHED DATA AT TT(I) CALL SMOP(S,X,N,TT(I),WID,NUE,IORD,IBOUN,IST,S1,C,Y(I)) END IF 100 CONTINUE !- RETURN END SUBROUTINE SMOP(S,X,N,TAU,WID,NUE,IORD,IBOUN,IST,S1,C,Y) !----------------------------------------------------------------------- ! SHORT-VERSION JANUARY 1997 ! ! PERFORMS ONE SMOOTHING STEP ! ! PARAMETERS : ! ! INPUT S(0:N) HALF POINT INTERPOLATION SEQUENCE ! INPUT X(N) DATA ! INPUT N LENGTH OF X ! INPUT TAU POINT WHERE FUNCTION IS ESTIMATED ! INPUT WID ONE SIDED BANDWIDTH ! INPUT NUE ORDER OF DERIVATIVE (0-4) ! INPUT IORD ORDER OF KERNEL POLYNOMIAL ! INPUT IBOUN TYPE OF BOUNDARY ! INPUT IST INDEX OF FIRST POINT OF SMOOTHING INTERVAL ! INPUT S1 LEFT BOUNDARY OF SMOOTHING INTERVAL ! INPUT C(7) KERNEL COEFFICIENTS ! OUTPUT Y SMOOTHED VALUE AT TAU ! WORK WO(7) WORK ARRAY ! !----------------------------------------------------------------------- DOUBLE PRECISION X(N),S(0:N),WO(7) DOUBLE PRECISION C(7) INTEGER N,NUE,IORD,IBOUN,IST,JEND,IBEG,INCR,I,J DOUBLE PRECISION TAU,WID,S1,Y,YY,YYY,W,WIDNUE,WW !- Y=0. WW=0. JEND=0 IBEG=2 IF(NUE.NE.1.AND.NUE.NE.3) IBEG=1 INCR=2 !- !------ COMPUTE INITIAL KERNEL VALUES IF(IBOUN.GT.0) THEN YY=(TAU-S1)/WID WO(IBEG)=YY YY=YY*YY IF(NUE.EQ.1.OR.NUE.EQ.3) WO(IBEG)=YY DO 1 I=IBEG,IORD-INCR,INCR 1 WO(I+INCR)=WO(I)*YY ELSE DO 2 I=IBEG,IORD,INCR 2 WO(I)=1. END IF !- !------ LOOP OVER SMOOTHING INTERVAL DO 100 J=IST,N YY=(TAU-S(J))/WID IF(YY.LT.-1.) THEN YY=-1. JEND=1 END IF YYY=YY YY=YY*YY IF(NUE.EQ.1.OR.NUE.EQ.3) YYY=YY !- !------ LOOP FOR COMPUTING WEIGHTS W=0. DO 3 I=IBEG,IORD,INCR W=W+C(I)*(WO(I)-YYY) WO(I)=YYY YYY=YYY*YY 3 CONTINUE Y=Y+W*X(J) WW=WW+W IF(JEND.EQ.1) GOTO 110 100 CONTINUE !- !------- NORMALIZING FOR NUE>0 110 IF(WW.NE.0) Y=Y/WW IF(NUE.GT.0) THEN WIDNUE=WID**NUE Y=Y/WIDNUE END IF !- RETURN END SUBROUTINE COFFI(NUE,KORD,C) !----------------------------------------------------------------------- ! SHORT-VERSION: JANUARY 1995 ! ! PURPOSE: ! ! DEFINES POLYNOMIAL KERNEL COEFFICIENTS FOR INTERIOR. ! ! PARAMETERS: ! ! INPUT NUE ORDER OF DERIVATIVE (0-4) ! INPUT KORD ORDER OF KERNEL (NUE+I, I=2,4,6; KORD<=6) ! OUTPUT C(7) POLYNOMIAL KERNEL COEFFICIENTS ! !----------------------------------------------------------------------- DOUBLE PRECISION C(7) !- DO 10 I=1,7 10 C(I)=0. IF(NUE.EQ.0.AND.KORD.EQ.2) THEN C(1)=0.75D0 C(3)=-0.25D0 END IF ! IF(NUE.EQ.0.AND.KORD.EQ.4) THEN C(1)=1.40625D0 C(3)=-1.5625D0 C(5)=0.65625D0 END IF ! IF(NUE.EQ.0.AND.KORD.EQ.6) THEN C(1)=2.05078125D0 C(3)=-4.78515625D0 C(5)=5.16796875D0 C(7)=-1.93359375D0 END IF ! IF(NUE.EQ.1.AND.KORD.EQ.3) THEN C(2)=-1.875D0 C(4)=0.9375D0 END IF ! IF(NUE.EQ.1.AND.KORD.EQ.5) THEN C(2)=-8.203125D0 C(4)=11.484375D0 C(6)=-4.921875D0 END IF ! IF(NUE.EQ.2.AND.KORD.EQ.4) THEN C(1)=-6.5625D0 C(3)=13.125D0 C(5)=-6.5625D0 END IF ! IF(NUE.EQ.2.AND.KORD.EQ.6) THEN C(1)=-24.609375D0 C(3)=103.359375D0 C(5)=-132.890625D0 C(7)=54.140625D0 END IF ! IF(NUE.EQ.3.AND.KORD.EQ.5) THEN C(2)=88.59375D0 C(4)=-147.65625D0 C(6)=68.90625D0 END IF ! IF(NUE.EQ.4.AND.KORD.EQ.6) THEN C(1)=324.84375D0 C(3)=-1624.21875D0 C(5)=2273.90625D0 C(7)=-974.53125D0 END IF ! RETURN END SUBROUTINE COFFB(NUE,KORD,Q,IBOUN,C) !----------------------------------------------------------------------- ! SHORT-VERSION: JANUARY 1995 ! ! PURPOSE: ! ! COMPUTES COEFFICIENTS OF POLYNOMIAL BOUNDARY KERNELS, ! FOLLOWING GASSER + MUELLER PREPRINT 38 SFB 123 HEIDELBERG ! AND UNPUBLISHED RESULTS ! ! PARAMETERS: ! ! INPUT NUE ORDER OF DERIVATIVE (0-4) ! INPUT KORD ORDER OF KERNEL (NUE+I, I=2,4,6; KORD<=6) ! INPUT Q PERCENTAGE OF WID AT BOUNDARY ! INPUT IBOUN < 0 RIGHT BOUNDARY OF DATA ! > 0 LEFT BOUNDARY OF DATA ! OUTPUT C(7) POLYNOMIAL KERNEL COEFFICIENTS ! !----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION C(7) ! DO 10 I=1,7 10 C(I)=0. P=-Q P1=1.+Q P3=P1*P1*P1 ! IF(NUE.EQ.0.AND.KORD.EQ.2) THEN D=1./(P3*P1) C(1)=(6.+P*(12.+P*18.))*D C(2)=9.*(1.+P)*(1.+P)*D C(3)=(4.+P*8.)*D END IF ! IF(NUE.EQ.0.AND.KORD.EQ.4) THEN D=P1/(P3*P3*P3) P12=(1.+P)*(1.+P)*D C(1)=20.*(1.+P*(12.+P*(78.+P*(164.+P*(165.+P*(60.+P*10.))))))*D C(2)=100.*(1.+P*(5.+P))**2*P12 C(3)=200.*(1.+P*(12.+P*(33.+P*(36.+P*(14.+P+P)))))*D C(4)=175.*(1.+P*(10.+P*3.))*P12 C(5)=56.*(1.+P*(12.+P*(18.+P*4.)))*D END IF ! IF(NUE.EQ.0.AND.KORD.EQ.6) THEN P6=P3*P3 D=1./(P6*P6) P12=(1.+P)*(1.+P)*D C(1)=42.*(1.+P*(30.+P*(465.+P*(3000.+P*(10050.+P*(17772.+P & *(17430.+P*(9240.+P*(2625.+P*(350.+P*21.))))))))))*D C(2)=441.*(1.+P*(14.+P*(36.+P*(14.+P))))**2*P12 C(3)=1960.*(1.+P*(30.+P*(255.+P*(984.+P*(1902.+P*(1956.+P & *(1065.+P*(300.+P*(39.+P+P)))))))))*D C(4)=4410.*(1.+P*(28.+P*(156.+P*(308.+P*(188.+P*(42.+P*3.)))))) & *P12 C(5)=5292.*(1.+P*(30.+P*(185.+P*(440.+P*(485.+P*(250.+P*(57. & +P*4.)))))))*D C(6)=3234.*(1.+P*(28.+P*(108.+P*(56.+P*5.))))*P12 C(7)=792.*(1.+P*(30.+P*(150.+P*(200.+P*(75.+P*6.)))))*D END IF ! IF(NUE.EQ.1.AND.KORD.EQ.3) THEN D=-1./(P3*P3) P12=(1.+P)*(1.+P)*D C(1)=(60.+P*240.)*P12 C(2)=120.*(2.+P*(6.+P*(6.+P)))*D C(3)=300.*P12 C(4)=(120.+P*180.)*D END IF ! IF(NUE.EQ.1.AND.KORD.EQ.5) THEN D=-1./(P3*P3*P3*P1) P12=(1.+P)*(1.+P)*D C(1)=420.*(1.+P*(18.+P*(98.+P*(176.+P*(75.+P*10.)))))*P12 C(2)=2100.*(2.+P*(25.+P*(120.+P*(245.+P*(238.+P*(105.+P*(20. & +P)))))))*D C(3)=14700.*(1.+P*(4.+P))**2*P12 C(4)=5880.*(4.+P*(35.+P*(90.+P*(95.+P*(40.+P*6.)))))*D C(5)=17640.*(1.+P*(6.+P+P))*P12 C(6)=2520.*(2.+P*(15.+P*(20.+P*5.)))*D END IF ! IF(NUE.EQ.2.AND.KORD.EQ.4) THEN D=P1/(P3*P3*P3) P12=(1.+P)*(1.+P)*D C(1)=840.*(1.+P*(12.+P*(28.+P*(24.+P*5.))))*D C(2)=2100.*(3.+P*(10.+P))*P12 C(3)=1680.*(9.+P*(28.+P*(27.+P*6.)))*D C(4)=14700.*P12 C(5)=(5040.+P*6720.)*D END IF ! IF(NUE.EQ.2.AND.KORD.EQ.6) THEN P6=P3*P3 D=1./(P6*P6) P12=(1.+P)*(1.+P)*D C(1)=5040.*(2.+P*(60.+P*(489.+P*(1786.+P*(3195.+P*(2952.+P & *(1365.+P*(294.+P*21.))))))))*D C(2)=52920.*(3.+P*(42.+P*(188.+P*(308.+P*(156.+P*(28. & +P))))))*P12 C(3)=141120.*(6.+P*(68.+P*(291.+P*(570.+P*(555.+P*(264.+P & *(57.+P*4.)))))))*D C(4)=529200.*(2.+P*(7.+P+P))**2*P12 C(5)=90720.*(30.+P*(228.+P*(559.+P*(582.+P*(255. & +P*40.)))))*D C(6)=582120.*(3.+P*(14.+P*5.))*P12 C(7)=221760.*(2.+P*(12.+P*(15.+P*4.)))*D END IF ! IF(NUE.EQ.3.AND.KORD.EQ.5) THEN D=-1./(P3*P3*P3*P1) P12=(1.+P)*(1.+P)*D C(1)=15120.*(1.+P*(18.+P*(38.+P*6.)))*P12 C(2)=45360.*(4.+P*(35.+P*(80.+P*(70.+P*(20.+P)))))*D C(3)=352800.*(2.+P*(6.+P))*P12 C(4)=151200.*(8.+P*(25.+P*(24.+P*6.)))*D C(5)=952560.*P12 C(6)=70560.*(4.+P*5.)*D END IF ! IF(NUE.EQ.4.AND.KORD.EQ.6) THEN P6=P3*P3 D=1./(P6*P6) P12=(1.+P)*(1.+P)*D C(1)=332640.*(1.+P*(30.+P*(171.+P*(340.+P*(285.+P*(90.+P*7. & ))))))*D C(2)=1164240.*(5.+P*(56.+P*(108.+P*(28.+P))))*P12 C(3)=6652800.*(5.+P*(38.+P*(85.+P*(76.+P*(25.+P+P)))))*D C(4)=17463600.*(5.+P*(14.+P*3.))*P12 C(5)=4656960.*(25.+P*(78.+P*(75.+P*20.)))*D C(6)=76839840.*P12 C(7)=3991680.*(5.+P*6.)*D END IF ! IF(IBOUN.GT.0) RETURN J=2 IF(NUE.EQ.1.OR.NUE.EQ.3) J=1 DO 2 I=J,KORD,2 2 C(I)=-C(I) RETURN END SUBROUTINE COFF(X,N,FA) DOUBLE PRECISION X(N),FA INTEGER I,N DO 1 I=1,N 1 X(I)=FA RETURN END ! !======================================================================================================================== ! ! Numerical Recipes file ! !======================================================================================================================== ! include 'ran.f90' ! include 'd:\nr7790\ran.f90'

Top of Message | Previous Page | Permalink

JiscMail Tools


RSS Feeds and Sharing


Advanced Options


Archives

November 2018
October 2018
September 2018
August 2018
July 2018
May 2018
April 2018
March 2018
February 2018
January 2018
December 2017
November 2017
October 2017
September 2017
August 2017
July 2017
June 2017
May 2017
April 2017
March 2017
January 2017
December 2016
November 2016
October 2016
September 2016
August 2016
July 2016
June 2016
May 2016
April 2016
March 2016
February 2016
December 2015
November 2015
October 2015
September 2015
August 2015
June 2015
April 2015
March 2015
January 2015
December 2014
November 2014
October 2014
August 2014
July 2014
May 2014
April 2014
March 2014
February 2014
January 2014
December 2013
November 2013
July 2013
June 2013
May 2013
April 2013
February 2013
January 2013
December 2012
November 2012
October 2012
September 2012
August 2012
July 2012
June 2012
April 2012
March 2012
February 2012
January 2012
December 2011
November 2011
October 2011
September 2011
August 2011
July 2011
June 2011
May 2011
April 2011
March 2011
February 2011
January 2011
December 2010
November 2010
October 2010
August 2010
July 2010
June 2010
March 2010
February 2010
January 2010
December 2009
October 2009
August 2009
July 2009
June 2009
March 2009
February 2009
January 2009
December 2008
November 2008
October 2008
September 2008
August 2008
July 2008
June 2008
May 2008
April 2008
March 2008
February 2008
December 2007
November 2007
October 2007
September 2007
August 2007
July 2007
June 2007
May 2007
April 2007
March 2007
January 2007
2006
2005
2004
2003
2002
2001
2000
1999
1998
1997


JiscMail is a Jisc service.

View our service policies at https://www.jiscmail.ac.uk/policyandsecurity/ and Jisc's privacy policy at https://www.jisc.ac.uk/website/privacy-notice

Secured by F-Secure Anti-Virus CataList Email List Search Powered by the LISTSERV Email List Manager