Logo Search packages:      
Sourcecode: octave-octgpr version File versions  Download package

infpgp.f

c Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
c 
c Author: Jaroslav Hajek <highegg@gmail.com>
c 
c This file is part of OctGPR.
c 
c OctGPR is free software; you can redistribute it and/or modify
c it under the terms of the GNU General Public License as published by
c the Free Software Foundation; either version 3 of the License, or
c (at your option) any later version.
c 
c This program is distributed in the hope that it will be useful,
c but WITHOUT ANY WARRANTY; without even the implied warranty of
c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c GNU General Public License for more details.
c 
c You should have received a copy of the GNU General Public License
c along with this software; see the file COPYING.  If not, see
c <http://www.gnu.org/licenses/>.
c 
      subroutine infpgp(ndim,nf,F,theta,nu,var,nlin,mu,QP,corr,
     +x0,y0,sig0,nder,yd0,work)
c purpose:      compute the prediction value, spatial derivatives 
c               and prediction variance of the PGP regressor in
c               a single spatial point. Should be used after
c               nllpgp, but *NOT* after nldpgp (on the same data).
c arguments:
c ndim (in)     number of dimensions of input space
c nf (in)       number of inducing vectors
c F (in)        array of inducing input vectors
c theta (in)    length scales
c nu (in)       relative white noise. nu = sqrt(var_white/var)
c var (out)     MLE estimated global variance
c nlin (in)     number of leading input variables to include in linear
c               mean trend. Set nlin = 0 to use a constant mean.
c mu (out)      at least (nlin+1)*(nlin+2). The first nlin+1 components
c               contain the components of the linear trend (constant
c               first). 
c QP (out)      at least nf*(nf+3). The array contains packed factorization 
c               details as returned from pakpgp
c nll (out)     the negative log-likelihood
c corr          subroutine to calculate correlation value and its
c               derivative. Must be declared as follows:
c               subroutine corr(t,f,d)
c               double precision t,f,d
c               f = correlation
c               d = derivative
c               end subroutine
c               The correlation should satisfy f(0) = 1 and f(+inf) = 0.
c               t >= 0 is the scaled squared norm of input vectors
c               difference, i.e. sum(theta*(X(:,i)-X(:,j))**2)
c x0 (in)       the spatial point to predict in
c y0 (out)      the prediction values.
c sig0 (out)    the prediction sigmas (noise included).
c nder(in)      number of derivatives requested. 0 to omit derivatives.
c yd0 (out)     the prediction derivatives. if nder <= 0, yd0 is not
c               referenced.
c work (out)    workspace; size at least nf*(2+min(nder,1))
c
      integer ndim,nf,nlin,info
      real*8 F(ndim,nf),theta(ndim),nu,var,x0(ndim)
      real*8 mu(0:nlin),QP(nf,*),y0,sig0,yd0(*),work(nf,*)
      external corr
      external dtrsv,dgemv,dwdis2,dsumsq,dsdacc
      integer i,j,iz
      real*8 dwdis2,dsumsq,sums,sum,l2pi,eps,tmp

      iz = nf+2
c calculate correlation vector r
      do i = 1,nf
        call corr(dwdis2(ndim,theta,F(1,i),x0),work(i,1),tmp)
        work(i,2) = work(i,1)
c only use last part of workspace if derivatives requested
        if (nder > 0) work(i,3) = tmp
      end do
c form LA \ r
      call dtrsv('U','T','N',nf,QP(1,2),nf,work(1,1),1)
c accumulate sum((L\r).^2) and (L\y)'*(L\r)
      call dsdacc(nf,work(1,1),QP(1,iz),sig0,tmp)
c add linear trend
      tmp = tmp + mu(0)
      do k = 1,nlin
        tmp = tmp + x0(k)*mu(k)
      end do
      y0 = tmp
c form LQ \ r
      call dtrsv('L','N','N',nf,QP(1,1),nf,work(1,2),1)
      sig0 = 1 + nu**2*(1+sig0) - dsumsq(nf,work(1,2),1)
c get deviation
      sig0 = sqrt(sig0 * var)
c calc derivatives only if necessary
      if (nder == 0) return
      do k = 1,nder
        yd0(k) = 0
      end do
      do i = 1,nf
c calculate the primary part
        tmp = QP(i,iz+1)*work(i,3)
c apply chain rule
        do k = 1,nder
          yd0(k) = yd0(k) + tmp*(x0(k) - F(k,i))
        end do
      end do
c correction
      do k = 1,nder
        yd0(k) = yd0(k) * (2 * theta(k)**2)
        if (k <= nlin) yd0(k) = yd0(k) + mu(k)
      end do
      end subroutine

Generated by  Doxygen 1.6.0   Back to index