-
Notifications
You must be signed in to change notification settings - Fork 9
/
Copy pathexponpen.R
53 lines (43 loc) · 1.36 KB
/
exponpen.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
exponpen <- function(basisobj, Lfdobj=int2Lfd(2))
{
# Computes the Exponential penalty matrix.
# Argument:
# BASISOBJ ... a basis.fd object of type "expon"
# LFDOBJ ... either the order of derivative or a
# linear differential operator to be penalized.
# Returns the penalty matrix.
# Last modified 9 February 2007
# Check BASISOBJ
if (!(inherits(basisobj, "basisfd"))) stop(
"First argument is not a basis object.")
type <- basisobj$type
if (type != "expon") stop ("Wrong basis type")
# Check LFDOBJ
Lfdobj <- int2Lfd(Lfdobj)
# Compute penalty matrix
if (is.integerLfd(Lfdobj)) {
nderiv <- Lfdobj$nderiv
ratevec <- basisobj$params
nrate <- length(ratevec)
penaltymatrix <- matrix(0,nrate,nrate)
tl <- basisobj$rangeval[1]
tu <- basisobj$rangeval[2]
for (irate in 1:nrate) {
ratei <- ratevec[irate]
for (jrate in 1:irate) {
ratej <- ratevec[jrate]
ratesum <- ratei + ratej
if (ratesum != 0) {
penaltymatrix[irate,jrate] <- (ratei*ratej)^nderiv *
(exp(ratesum*tu) - exp(ratesum*tl)) / ratesum
} else {
if (nderiv == 0) penaltymatrix[irate,jrate] <- tu - tl
}
penaltymatrix[jrate,irate] <- penaltymatrix[irate,jrate]
}
}
} else {
penaltymatrix <- inprod(basisobj, basisobj, Lfdobj, Lfdobj)
}
penaltymatrix
}