Skip to content

Commit

Permalink
first commit
Browse files Browse the repository at this point in the history
  • Loading branch information
B. W. Lewis committed Apr 29, 2010
0 parents commit e9f0195
Show file tree
Hide file tree
Showing 7 changed files with 303 additions and 0 deletions.
10 changes: 10 additions & 0 deletions DESCRIPTION
@@ -0,0 +1,10 @@
Package: fls
Type: Package
Title: Flexible least squares and related weighted Kalman filter
Version: 0.1
Date: 2009-11-02
Author: Bryan W. Lewis
Maintainer: Bryan W. Lewis <bwaynelewis@gmail.com>
Description: Implementation of the Kalaba-Tesfatsion flexible least squares algorithm and a related weighted Kalman filter of Montana, et. al. See the documentation for additional references.
License: GPL
LazyLoad: yes
42 changes: 42 additions & 0 deletions R/fls.R
@@ -0,0 +1,42 @@
`fls` <-
function (A, b, mu=1, ncap=length(b), smoothed=TRUE)
{
m <- nrow (A)
n <- ncol (A)
M <- array (0,c(n,n,ncap))
E <- array (0,c(n,ncap))
X <- array (0,c(n,ncap))
R <- matrix(0,n,n)
diag(R) <- diag(R) + mu
for (j in 1:ncap) {
Z <- solve(qr(R + tcrossprod(A[j,]),LAPACK=TRUE),diag(1.0,n));
M[,,j] <- mu*Z # (5.7b)
v <- b[j]*A[j,]
if(j==1) p <- rep(0,n)
else p <- mu*E[,j-1]
w <- p + v
E[,j] <- Z %*% w # (5.7c)
R <- -mu*mu*Z
diag(R) <- diag(R) + 2*mu
}
# Calculate eqn (5.15) FLS estimate at ncap
Q <- -mu*M[,,ncap-1]
diag(Q) <- diag(Q) + mu
Ancap <- A[ncap,,drop=FALSE]
C <- Q + t(Ancap) %*% Ancap
d <- mu*E[,ncap-1,drop=FALSE] + b[ncap]*t(Ancap)
X[,ncap] <- C %*% d
X[,ncap] <- solve(qr(C,LAPACK=TRUE),d)
if (smoothed) {
# Use eqn (5.16) to obtain smoothed FLS estimates for
# X[,1], X[,2], ..., X[,ncap-1]
for (j in 1:(ncap-1)) {
l <- ncap - j
X[,l] <- E[,l] + M[,,l] %*% X[,l+1]
}
}
else {
X <- X[,ncap]
}
X
}
26 changes: 26 additions & 0 deletions R/flskf.R
@@ -0,0 +1,26 @@
`flskf` <-
function (A, b, mu=1, ncap=length(b))
{
m <- nrow (A)
n <- ncol (A)
X <- array (0,c(n,ncap))
R <- array (0,c(n,n))
diag (R) <- 1/mu
a <- t(A[1,,drop=FALSE])
phi <- c((1/mu) * crossprod(a) + 1)
w <- (1/phi) * R %*% a
rho <- b[1]
X[,1] <- rho * w
for (j in 2:ncap) {
R = R - phi * tcrossprod(w)
diag (R) <- diag (R) + 1/mu
a <- t(A[j,,drop=FALSE])
phi <- c(t(a) %*% R %*% a)
phi <- phi + 1
rho <- c(b[j] - crossprod (a, X[,j-1]))
w <- (1/phi) * R %*% a
X[,j] <- X[,j-1] + rho*w
}
X
}

5 changes: 5 additions & 0 deletions README
@@ -0,0 +1,5 @@
The Kalaba-Tesfatsion Flexible Least Squares method.
Numbered equations refer to the paper:
R. Kalaba and L. Tesfatsion, Time-Varying Linear Regression Via
Flexible Least Squares, Computers and Mathematics with Applications,
Vol. 17 (1989), pp. 1215-1245.
31 changes: 31 additions & 0 deletions man/fls-package.Rd
@@ -0,0 +1,31 @@
\name{fls-package}
\alias{fls-package}
\alias{fls}
\docType{package}
\title{
What the package does (short line)
}
\description{
More about what it does (maybe more than one line)
}
\details{
}
\author{
Who wrote it

Maintainer: Who to complain to <yourfault@somewhere.net>
~~ The author and/or maintainer of the package ~~
}
\references{
~~ Literature or other references for background information ~~
}
~~ Optionally other standard keywords, one per line, from file KEYWORDS in ~~
~~ the R documentation directory ~~
\keyword{ package }
\seealso{
~~ Optional links to other man pages, e.g. ~~
~~ \code{\link[<pkg>:<pkg>-package]{<pkg>}} ~~
}
\examples{
~~ simple examples of the most important functions ~~
}
101 changes: 101 additions & 0 deletions man/fls.Rd
@@ -0,0 +1,101 @@
\name{fls}
\Rdversion{1.1}
\alias{fls}
\title{Flexible Least Squares}
\description{
The Kalaba-Tesfatsion Flexible Least Squares method.
Numbered equations refer to the paper:
R. Kalaba and L. Tesfatsion, Time-Varying Linear Regression Via
Flexible Least Squares, Computers and Mathematics with Applications,
Vol. 17 (1989), pp. 1215-1245.
}
\usage{
fls(A, b, mu = 1, ncap = length(b))
}
%- maybe also 'usage' for other objects documented here.
\arguments{
\item{A}{
%% ~~Describe \code{A} here~~
}
\item{b}{
%% ~~Describe \code{b} here~~
}
\item{mu}{
%% ~~Describe \code{mu} here~~
}
\item{ncap}{
%% ~~Describe \code{ncap} here~~
}
}
\details{
%% ~~ If necessary, more details than the description above ~~
}
\value{
%% ~Describe the value returned
%% If it is a LIST, use
%% \item{comp1 }{Description of 'comp1'}
%% \item{comp2 }{Description of 'comp2'}
%% ...
}
\references{
%% ~put references to the literature/web site here ~
}
\author{
%% ~~who you are~~
}
\note{
%% ~~further notes~~
}

%% ~Make other sections like Warning with \section{Warning }{....} ~

\seealso{
%% ~~objects to See Also as \code{\link{help}}, ~~~
}
\examples{
##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.

## The function is currently defined as
function (A, b, mu=1, ncap=length(b))
{
m <- nrow (A)
n <- ncol (A)
M <- array (0,c(n,n,ncap))
E <- array (0,c(n,ncap))
X <- array (0,c(n,ncap))
R <- matrix(0,n,n)
diag(R) <- diag(R) + mu
for (j in 1:ncap) {
Z <- solve(qr(R + tcrossprod(A[j,]),LAPACK=TRUE),diag(1.0,n));
M[,,j] <- mu*Z # (5.7b)
v <- b[j]*A[j,]
if(j==1) p <- rep(0,n)
else p <- mu*E[,j-1]
w <- p + v
E[,j] <- Z \%*\% w # (5.7c)
R <- -mu*mu*Z
diag(R) <- diag(R) + 2*mu
}
# Calculate eqn (5.15) FLS estimate at ncap
Q <- -mu*M[,,ncap-1]
diag(Q) <- diag(Q) + mu
Ancap <- A[ncap,,drop=FALSE]
C <- Q + t(Ancap) \%*\% Ancap
d <- mu*E[,ncap-1,drop=FALSE] + b[ncap]*t(Ancap)
X[,ncap] <- C \%*\% d
X[,ncap] <- solve(qr(C,LAPACK=TRUE),d)
# Use eqn (5.16) to obtain smoothed FLS estimates for
# X[,1], X[,2], ..., X[,ncap-1]
for (j in 1:(ncap-1)) {
l <- ncap - j
X[,l] <- E[,l] + M[,,l] \%*\% X[,l+1]
}
X
}
}
% Add one or more standard keywords, see file 'KEYWORDS' in the
% R documentation directory.
\keyword{ ~kwd1 }
\keyword{ ~kwd2 }% __ONLY ONE__ keyword per line
88 changes: 88 additions & 0 deletions man/flskf.Rd
@@ -0,0 +1,88 @@
\name{flskf}
\Rdversion{1.1}
\alias{flskf}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
%% ~~function to do ... ~~
}
\description{
%% ~~ A concise (1-5 lines) description of what the function does. ~~
}
\usage{
flskf(A, b, mu = 1, ncap = length(b))
}
%- maybe also 'usage' for other objects documented here.
\arguments{
\item{A}{
%% ~~Describe \code{A} here~~
}
\item{b}{
%% ~~Describe \code{b} here~~
}
\item{mu}{
%% ~~Describe \code{mu} here~~
}
\item{ncap}{
%% ~~Describe \code{ncap} here~~
}
}
\details{
%% ~~ If necessary, more details than the description above ~~
}
\value{
%% ~Describe the value returned
%% If it is a LIST, use
%% \item{comp1 }{Description of 'comp1'}
%% \item{comp2 }{Description of 'comp2'}
%% ...
}
\references{
%% ~put references to the literature/web site here ~
}
\author{
%% ~~who you are~~
}
\note{
%% ~~further notes~~
}

%% ~Make other sections like Warning with \section{Warning }{....} ~

\seealso{
%% ~~objects to See Also as \code{\link{help}}, ~~~
}
\examples{
##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.

## The function is currently defined as
function (A, b, mu=1, ncap=length(b))
{
m <- nrow (A)
n <- ncol (A)
X <- array (0,c(n,ncap))
R <- array (0,c(n,n))
diag (R) <- 1/mu
a <- t(A[1,,drop=FALSE])
phi <- c((1/mu) * crossprod(a) + 1)
w <- (1/phi) * R \%*\% a
rho <- b[1]
X[,1] <- rho * w
for (j in 2:ncap) {
R = R - phi * tcrossprod(w)
diag (R) <- diag (R) + 1/mu
a <- t(A[j,,drop=FALSE])
phi <- c(t(a) \%*\% R \%*\% a)
phi <- phi + 1
rho <- c(b[j] - crossprod (a, X[,j-1]))
w <- (1/phi) * R \%*\% a
X[,j] <- X[,j-1] + rho*w
}
X
}
}
% Add one or more standard keywords, see file 'KEYWORDS' in the
% R documentation directory.
\keyword{ ~kwd1 }
\keyword{ ~kwd2 }% __ONLY ONE__ keyword per line

0 comments on commit e9f0195

Please sign in to comment.