Bilinear Interpolation with extrapolation
octaedro
use in a dynamic programming choice problem. It seems to work fine the
interpolation part but when extrapolating for some extreme values it
performs really bad throwing numebers that are not feasible. Can
anybody tell me what am I doing wrong or it there is some routine that
works and it is not propriety.
Thank you
*****************************************************************
*****************************************************************
SUBROUTINE sub_linintp2(n1,n2,y,xx1,xx2,yy,x1,x2)
! SUB_LININTP2.F90 performs 2-D linear interpolation over a given
abcisa
! ARGUMENTS :
! Y : Interpolated Value YY : Coordinates Data
! XX : Abcisas Data X : Evaluation Point
USE nrtype
IMPLICIT NONE
REAL(DP), INTENT(OUT) :: y
INTEGER, INTENT(IN) :: n1,n2
REAL(DP), INTENT(IN) :: x1,x2
REAL(DP), INTENT(IN) :: xx1(n1),xx2(n2),yy(n1,n2)
REAL(DP) :: t,u
INTEGER :: in1,in2
! THIS IS A PRELIMINARY EXTRAPOLATION STEP
! EXTRAPOLATION FROM BELOW
IF ((x1 .lt. xx1(1)) .and. (x2 .lt. xx2(1))) THEN
t=(x1-xx1(1))/(xx1(2)-xx1(1))
u=(x2-xx2(1))/(xx2(2)-xx2(1))
y=(1.0-t)*(1.0-u)*yy(1,1)+t*(1.0-u)*yy(2,1)+t*u*yy(2,2)+(1.0-t)*u*yy
(1,2)
RETURN
ELSE IF (x1 .lt. xx1(1)) THEN
t=(x1-xx1(1))/(xx1(2)-xx1(1))
DO in2=1,n2-1
IF ((x2 .ge. xx2(in2)) .and. (x2 .le. xx2(in2+1))) THEN
u=(x2-xx2(in2))/(xx2(in2+1)-xx2(in2))
y=(1.0-t)*(1.0-u)*yy(1,in2)+t*(1.0-u)*yy(2,in2)+t*u*yy
(2,in2+1)+(1.0-t)*u*yy(1,in2+1)
RETURN
END IF
END DO
ELSE IF (x2 .lt. xx2(1)) THEN
u=(x2-xx2(1))/(xx2(2)-xx2(1))
DO in1=1,n1-1
IF ((x1 .ge. xx1(in1)) .and. (x1 .le. xx1(in1+1))) THEN
t=(x1-xx1(in1))/(xx1(in1+1)-xx1(in1))
y=(1.0-t)*(1.0-u)*yy(in1,1)+t*(1.0-u)*yy(in1+1,1)+t*u*yy
(in1+1,2)+(1.0-t)*u*yy(in1,2)
RETURN
END IF
END DO
END IF
! EXTRAPOLATION FROM ABOVE
IF ((x1 .gt. xx1(n1)) .and. (x2 .gt. xx2(n2))) THEN
t=(x1-xx1(n1))/(xx1(n1)-xx1(n1-1))
u=(x2-xx2(n2))/(xx2(n2)-xx2(n2-1))
y=(1.0-t)*(1.0-u)*yy(n1-1,n2-1)+t*(1.0-u)*yy(n1,n2-1)+t*u*yy
(n1,n2)+(1.0-t)*u*yy(n1-1,n2)
RETURN
ELSE IF (x1 .gt. xx1(n1)) THEN
t=(x1-xx1(n1))/(xx1(n1)-xx1(n1-1))
DO in2=1,n2-1
IF ((x2 .ge. xx2(in2)) .and. (x2 .le. xx2(in2+1))) THEN
u=(x2-xx2(in2))/(xx2(in2+1)-xx2(in2))
y=(1.0-t)*(1.0-u)*yy(n1-1,in2)+t*(1.0-u)*yy(n1,in2)+t*u*yy
(n1,in2+1)+(1.0-t)*u*yy(n1-1,in2+1)
RETURN
END IF
END DO
ELSE IF (x2 .gt. xx2(n2)) THEN
u=(x2-xx2(n2))/(xx2(n2)-xx2(n2-1))
DO in1=1,n1-1
IF ((x1 .ge. xx1(in1)) .and. (x1 .le. xx1(in1+1))) THEN
t=(x1-xx1(in1))/(xx1(in1+1)-xx1(in1))
y=(1.0-t)*(1.0-u)*yy(in1,n2-1)+t*(1.0-u)*yy(in1+1,n2-1)+t*u*yy
(in1+1,n2)+(1.0-t)*u*yy(in1,n2)
RETURN
END IF
END DO
END IF
! INTERPOLATION INSIDE THE "UNIT" SQUARE
DO in1=1,n1-1
DO in2=1,n2-1
IF (((x1 .ge. xx1(in1)) .and. (x1 .le. xx1(in1+1))) .and.
((x2 .ge. xx2(in2)) .and. (x2 .le. xx2(in2+1)))) THEN
t=(x1-xx1(in1))/(xx1(in1+1)-xx1(in1))
u=(x2-xx2(in2))/(xx2(in2+1)-xx2(in2))
y=(1.0-t)*(1.0-u)*yy(in1,in2)+t*(1.0-u)*yy
(in1+1,in2)+t*u*yy(in1+1,in2+1)+(1.0-t)*u*yy(in1,in2+1)
RETURN
END IF
END DO
END DO
END SUBROUTINE sub_linintp2
Gordon Sande
> I wrote a simple code to interpolate two dimensional functions that I
> use in a dynamic programming choice problem. It seems to work fine the
> interpolation part but when extrapolating for some extreme values it
> performs really bad throwing numebers that are not feasible. Can
> anybody tell me what am I doing wrong or it there is some routine that
> works and it is not propriety.
>
> Thank you
Most polynomial interpolation schemes are intended to work well
within some range at the price of exploding wildly once outside
that range. The usual interpolation schmemes are just shorcut ways
of finding the polynomial.
Try evaluating a Tchebycheff polynomial outside of [-1,+1] and see
for yourself. There is a theorem of classical mathematics (i.e real
stuff from circa 1890!) that these are the extreme case.
Maybe it is working perfectly and it is your expectation that is wrong.
To do extrapolation you usually need to fix the asymptopic form and work
back from there.
And by the way you should have used sci.mat,num-analysis and waded through
all the kiddies peddling solution manuals. I have never figured out why
s.m.n-a is the chosen bulletin board for that stuff. :-(
glen herrmannsfeldt
< I wrote a simple code to interpolate two dimensional functions that I
< use in a dynamic programming choice problem. It seems to work fine the
< interpolation part but when extrapolating for some extreme values it
< performs really bad throwing numebers that are not feasible. Can
< anybody tell me what am I doing wrong or it there is some routine that
< works and it is not propriety.
I just looked at the first case, and it seems fine to me.
I would think that it wouldn't take so many if/then/else cases, but
once the appropriate interpolation/extrapolation points were selected,
the expansion would be the same in all cases.
Can you show the 'wrong' results that it generates?
-- glen
Arjan
Management summary: don't!
Beyond the boundaries, you may adopt the value at the boundary.
Or give a warning/error-code.
Arjan
glen herrmannsfeldt
< On 2009-07-18 13:58:52 -0300, octaedro <jorge.al...@gmail.com> said:
<> I wrote a simple code to interpolate two dimensional functions that I
<> use in a dynamic programming choice problem. It seems to work fine the
<> interpolation part but when extrapolating for some extreme values it
<> performs really bad throwing numebers that are not feasible.
< Most polynomial interpolation schemes are intended to work well
< within some range at the price of exploding wildly once outside
< that range. The usual interpolation schmemes are just shorcut ways
< of finding the polynomial.
It seems to be linear extrapolation which should not explode
wildly as higher order expansions do. It should diverge linearly!
Though the routine seems to be much more complicated than it
should be, which makes it more difficult to notice problems.
It should just be two loops to find the points used in the
interpolation/extrapolation, two points in each direction for
linear, and then one actual expansion to generate the result.
-- glen
fj
As Arjan and Gordon, I advise you to never extrapolate. As mentioned
by Arjan, you might adopt the value at the boundary in case of
extrapolation.
From the programming point of view, as Glen noticed, your solution is
neither elegant nor efficient because you try to interpolate (or
extrapolate) along the two coordinates simultaneously. The normal way
consists in finding the right interval of xx1 containing x1 and,
independently, the right interval of xx2 containing x2. Look at your
two enclosed loops : the average CPU cost is proportional to (n1*n2)/2
instead of n1/2 + n2/2 when dealing with the two coordinates
separately.
Notice also that you might take advantage of the fact that the arrays
xx1 and xx2 are probably ordered. If the size of the arrays xx1 and
xx2 is sufficiently high (> 20), then a dichotomic search is more
efficient than a loop : Log2(n) instead of n/2.
At last, if you need to interpolate very often and if the coordinates
x1 and x2 change slowly, then it may be still more efficient to keep
the indexes i1 and i2 of the previous call and to start the search
from these indexes ...
mecej4
"octaedro" <jorge.al...@gmail.com> wrote in message
news:6eca28cf-6af7-4fc6...@v37g2000prg.googlegroups.com...
>I wrote a simple code to interpolate two dimensional functions that I
> use in a dynamic programming choice problem. It seems to work fine the
> interpolation part but when extrapolating for some extreme values it
> performs really bad throwing numebers that are not feasible. Can
> anybody tell me what am I doing wrong or it there is some routine that
> works and it is not propriety.
>
> Thank you
>
Matlab comes with an example from which one learns the dangers of polynomial
extrapolation. Using the census data for the US, for a certain choice of
degree the extrapolated value of today's population is NEGATIVE! -- and that
is before the Chinese thought of annihilating us.
Read also about Wilkinson's "perfidious polynomial".
-- mecej4
octaedro
stuff without enough theory background. I am aware about the problems
of extrapolation but in my case it is important how you deal with
those misfits as my problem has to do with wealth accumulation and
labor productivity those that may end up out of the pre-defined grid
may matter.
I agree with the lack of elegance of my routine but it was working
(most of the times)
I try to write this one but I don´t know why it does not work
SUBROUTINE sub_linint2d(n1,n2,y,xx1,xx2,yy,x1,x2)
USE nrtype
IMPLICIT NONE
REAL(DP), INTENT(OUT) :: y
INTEGER, INTENT(IN) :: n1,n2
REAL(DP), INTENT(IN) :: x1,x2
REAL(DP), INTENT(IN) :: xx1(n1),xx2(n2),yy(n1,n2)
REAL(DP) :: t,u
INTEGER :: in1,in2
in1=MINLOC(xx1-x1,DIM=1)
in2=MINLOC(xx2-x2,DIM=1)
IF (x1 .le. xx1(in1) .and. in1 .gt. 1) in1=in1-1
IF (x2 .le. xx2(in2) .and. in2 .gt. 1) in2=in2-1
t=(x1-xx1(in1))/(xx1(in1+1)-xx1(in1))
u=(x2-xx2(in2))/(xx2(in2+1)-xx2(in2))
y=(1.0-t)*(1.0-u)*yy(in1,in2)+t*(1.0-u)*yy(in1+1,in2)+t*u*yy
(in1+1,in2+1)+(1.0-t)*u*yy(in1,in2+1)
END SUBROUTINE sub_linint2d
glen herrmannsfeldt
< Thank you very much for your suggestions. It is hard to do numerical
< stuff without enough theory background. I am aware about the problems
< of extrapolation but in my case it is important how you deal with
< those misfits as my problem has to do with wealth accumulation and
< labor productivity those that may end up out of the pre-defined grid
< may matter.
It seems that you have non-uniform spacing on xx1 and xx2.
If the first or last two points are very closely spaced then
even linear extrapolation can give very bad results.
(snip)
< in1=MINLOC(xx1-x1,DIM=1)
< in2=MINLOC(xx2-x2,DIM=1)
I don't believe this is what you want. MINLOC gives
the location of the lowest valued array element,
positive or negative. Subtracting x1 doesn't change that.
-- glen
mecej4
>>Thank you very much for your suggestions. It is hard to do numerical
>>stuff without enough theory background. I am aware about the problems
>>of extrapolation but in my case it is important how you deal with
>>those misfits as my problem has to do with wealth accumulation and
>>labor productivity those that may end up out of the pre-defined grid
>>may matter.
>>
>>I agree with the lack of elegance of my routine but it was working
>>(most of the times)
<--CUT Code-->
If you can somehow avoid extrapolation, do so. If not, it is well worth the
effort to compute confidence intervals for the extrapolated values.
Look up PCHIP, "shape-preserving splines" and "Akima Interpolation". For the
last, Matlab files for the 1-D problem are available at
http://www.mathworks.com/matlabcentral/fileexchange/1814
-- mecej4
Arjan
> of extrapolation but in my case it is important how you deal with
> those misfits as my problem has to do with wealth accumulation and
> labor productivity those that may end up out of the pre-defined grid
> may matter.
No matter the field of application, whether you are estimating future
sea level rise, population growth, stock-option pricing or wealth
accumulation,
a proper set of boundaries is essential. In all of these types of
scenario
analysis, unexpected results may arise. In all cases extrapolation
beyond
the assumed boundaries will give dramatic, non-realistic effects.
Once you find yourself extrapolating, re-run your exercise
with broader boundaries!
Arjan
fj
Replace
in1=MINLOC(xx1-x1,DIM=1)
in2=MINLOC(xx2-x2,DIM=1)
by
in1=IFINR(n1,xx1,x1)
in2=IFINR(n2,xx2,x2)
IFINR being a function I wrote a long time ago in strict FORTRAN 77
(before the arrival of F90). Sorry : the comments are in French. But
the algorithm is not difficult to understand.
FUNCTION IFINR (N ,T ,X )
C&F FONCTION : RECHERCHE PAR DICHOTOTOMIE DU NUMERO
C&F DE L'INTERVALLE CONTENANT X :
C&F T(IFINR) < X < T(IFINR+1)
C&F SI X EST EN DEHORS DE LA TABLE , LE RESULTAT SERA :
C&F IFINR=1 SI X < T(1)
C&F IFINR=N-1 SI X > T(N)
C&A ARGUMENTS :
C&A N : NOMBRE DE POINTS DE LA TABLE
C&A T : TABLE RANGEE DANS UN ORDRE CROISSANT
C&A X : ABSCISSE RECHERCHEE
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION T(*)
IF(N.LE.0) THEN
IFINR=1
RETURN
ENDIF
IF(X.LE.T(1)) THEN
IFINR=1
ELSE IF(X.GE.T(N)) THEN
IFINR=N-1
ELSE
NINF=1
NSUP=N
DO 50 K=1,N+1
IF(NSUP-NINF.EQ.1) THEN
IFINR=NINF
RETURN
ENDIF
IFINR=(NINF+NSUP)/2
IF(T(IFINR).GE.X) THEN
NSUP=IFINR
ELSE
NINF=IFINR
ENDIF
50 CONTINUE
ENDIF
END