Nothing
#:::::::::::
# dtrev
#:::::::::::
subroutine dtrev (vmu, t, ldt, n, z, score, varht, info, work)
# Acronym: Double-precision TRidiagonal EValuation.
# Purpose: To compute the GCV/GML function and the related variance
# estimate from the tridiagonal matrix `t' and data vector `z'.
# References: 1. Gu, Bates, Chen, and Wahba(1988), TR#823, Stat, UW-M.
# 2. Dongarra et al. (1979) LINPACK User's Guide. (Chap. 4)
character*1 vmu
integer n, info
double precision t(ldt,*), z(*), score, varht, work(*)
# On entry:
# vmu 'v': GCV.
# 'm': GML.
# 'u': unbiased risk estimate.
# t the positive definite tridiagonal matrix T,
# stored in packed form:
# t(1,2:n): off-diagonal
# t(2,1:n): diagonal.
# ldt leading dimension of t.
# n the dimension of the matrix.
# z the appropriately transformed data vector.
# varht known variance if vmu=='u'.
# On exit:
# score the GCV/GML/URE score.
# varht \hat\sigma^{2}.
# info -3: vmu is none of 'v', 'm', or 'u'.
# > -3: as from LINPACK's `dpbfa'.
# Work array:
# work of size at least (n).
# Routines called directly:
# Fortran -- dexp, dble, dlog
# Blas -- dasum, dcopy, ddot, dscal
# Linpack -- dpbfa, dpbsl
# Written: Chong Gu, Statistics, UW-Madison, latest version 12/29/91.
double precision nume, deno, tmp, alph, la, dasum, ddot
integer j
info = 0
# check vmu
if ( vmu != 'v' & vmu != 'm' & vmu != 'u' ) {
info = -3
return
}
la = t(1,1)
# standardize the matrix for numerical stability
alph = dble (n) / dasum (n, t(2,1), ldt)
call dscal (n, alph, t(2,1), ldt)
call dscal (n-1, alph, t(1,2), ldt)
# decomposition
call dpbfa (t, ldt, n, 1, info)
if ( info != 0 ) return
call dcopy (n, z, 1, work, 1)
call dpbsl (t, ldt, n, 1, work)
# GCV computation
if ( vmu == 'v' ) {
tmp = 1.d0 / t(2,n) / t(2,n)
deno = tmp
for (j=n-1;j>0;j=j-1) {
tmp = ( 1.d0 + t(1,j+1) * t(1,j+1) * tmp ) / t(2,j) / t(2,j)
deno = deno + tmp
}
nume = ddot (n, work, 1, work, 1) / dble (n)
deno = deno / dble (n)
varht = alph * la * nume / deno
score = nume / deno / deno
}
# GML computation
if ( vmu == 'm' ) {
deno = dlog (t(2,n))
for (j=n-1;j>0;j=j-1) deno = deno + dlog (t(2,j))
nume = ddot (n, z, 1, work, 1) / dble (n)
varht = alph * la * nume
score = nume * dexp (2.d0 * deno / dble (n))
}
# unbiased risk computation
if ( vmu == 'u' ) {
nume = ddot (n, work, 1, work, 1) / dble (n)
tmp = 1.d0 / t(2,n) / t(2,n)
deno = tmp
for (j=n-1;j>0;j=j-1) {
tmp = ( 1.d0 + t(1,j+1) * t(1,j+1) * tmp ) / t(2,j) / t(2,j)
deno = deno + tmp
}
deno = deno / dble (n)
score = alph * alph * la * la * nume - 2.d0 * varht * alph * la * deno
}
return
end
#...............................................................................
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.