Commit ea1f6c7e authored by Reinhold Kainhofer's avatar Reinhold Kainhofer

Initial commit

parents
^.*\.Rproj$
^\.Rproj\.user$
.Rproj.user
.Rhistory
.RData
Package: RKTools
Type: Package
Title: Several utility functions by Reinhold Kainhofer: LogListPlot
Version: 1.0
Date: 2013-04-29
Author: Reinhold Kainhofer
Maintainer: Reinhold Kainhofer <office@open-tools.net>
Description: Utility functions for R, written by R. Kainhofer:
-) LogListPlot: Print (multiple) lists as a log-listplot
License: GPL (>=3)
exportPattern("^[[:alpha:]]+")
dbg = function(text, value, ...) {
cat(text, ": ");
if (is.matrix(value)) { cat(str(value)); }
else if (is.vector(value)) { cat(str(value)); }
else if (is.na(value)) { cat(value); }
else { cat(value); }
}
svd.composition = function(s,order=1) {
# Compose the matrix from the SVD factors, using only the first few dimensions
D = diag(s$d[1:order], nrow=order)
# need to use as.matrix, because if order=1, a vector is returned rather than a matrix
U = as.matrix(s$u[,1:order])
V = as.matrix(s$v[,1:order])
res=U %*% D %*% t(V) # X = U D V'
res
}
LeeCarter.simple = function(data, years=1:dim(data)[1], normalize=FALSE)
{
## FIXME: Add some sanity checks! ERROR messages
startyear = min(years);
endyear = max(years);
data=data[years,];
#print(data);
# Vector alpha is the mean for each age over the whole observation period
# Matrix Z is the residual containing the trend kappa_t and age-dependent strength beta_x
ax = colMeans(data,na.rm=TRUE)
dbg("ax", ax)
## calculate the residual by removing the mean from all observations:
Z = sweep(data, 2, ax)
# Consistency check: colMeans(Z)==0
dbg("colmeans(Z)", colMeans(Z));
# Normalize all ages to have roughly the same variance 1, if desired
stddev=rep(1,ncol(data))
if (normalize) {
s = svd(Z,1) # only need one dimension here
stddev=as.vector(sqrt(colMeans((Z-svd.composition(s,1))^2)))
dbg("stddev", stddev);
}
Znorm=Z;
for (i in 1:nrow(Znorm)) {
Znorm[i,] = Znorm[i,] / stddev;
}
# dbg("Z before apply", Z);
# dbg("stddev", stddev);
# Znorm=apply (Z, 1, function(q) {
# dbg("q ", q);
# dbg("q/stddev", q/stddev);
# q/stddev })
wireframe(Znorm)
print(Znorm)
# Singular Value Decomposition of matrix Z
s = svd(Znorm);
## TODO: include some analysis how much of the variation is covered by the first dimension!
# Lee-Carter vectors obtained from the SVD:
bx = s$v[,1] * stddev;
kt = s$d[1] * s$u[,1];
dbg("orig ax", ax)
dbg("orig bx", bx)
dbg("orig kt", kt)
# NORMALIZATION
# 1) Normalize to sum(bx^2)==1
# The sign(sum(bx)) is needed to ensure all bx>0 and kt is decreasing!
normalization = sqrt(sum(bx^2)) * sign(sum(bx));
print(normalization)
bx = bx / normalization;
kt = kt * normalization;
dbg("First normalization ||b||", normalization)
dbg("new bx ", bx)
dbg("new kt ", kt)
# 2) Normalize to sum(kt)==0
kmean = mean(kt)
# print(kmean);
kt = kt - kmean;
ax = ax + bx*kmean;
dbg("2nd normalization sum(kt)==00", kmean)
dbg("new bx ", bx)
dbg("new kt ", kt)
residuals=Znorm-as.matrix(kt) %*% t(bx)
result = list(ax = as.vector(ax),
bx = bx,
kt = kt,
sigmax = stddev,
normalize=normalize,
singular.values=s$d,
variance.captured = "TODO",
znorm=Znorm,
residuals=residuals)
result
}
## bx=1:5 / sqrt(55); kt=-5:5
#mx=as.matrix((-5:5)) %*% as.matrix(t(1:5)/sqrt(55))
## disturb by independent normal variates
#mx1 = mx + rnorm(n=mx,sd=(0.03*1:length(mx)))
#mx1
#mx
#wireframe(mx1)
##LeeCarter.simple(mx, normalize=TRUE)
#lc.austria.f=LeeCarter.simple(mx1, normalize=FALSE)
#lc.austria.f
#wireframe(lc.austria.f$znorm,scales = list(arrows = FALSE), colorkey = TRUE, drape=TRUE, screen = list(z = 30, x = -60))
#wireframe(lc.austria.f$residuals,scales = list(arrows = FALSE), colorkey = TRUE, drape=TRUE, screen = list(z = 30, x = -60))
#lc.austria.f.normalize=LeeCarter.simple(mx1, normalize=TRUE)
#lc.austria.f.normalize
#lc.austria.f=LeeCarter.simple(log(austria.m_fortgeschrieben[,1:96]), years=47:112)
#lc.austria.f
#lc.austria.f.normalize=LeeCarter.simple(log(austria.m_fortgeschrieben[,1:96]), years=47:112, normalize=TRUE)
#lc.austria.f.normalize
#wireframe(lc.austria.f$znorm)
#wireframe(lc.austria.f$znorm[,20:96],scales = list(arrows = FALSE), colorkey = TRUE, drape=TRUE, screen = list(z = 30, x = -60))
#wireframe(lc.austria.f$residuals[,20:96], scales = list(arrows = FALSE), colorkey = TRUE, drape=TRUE, screen = list(z = 30, x = -60))
#par(mfrow=c(1,3))
#ListPlot(list(lc.austria.f$ax, lc.austria.f.normalize$ax), xlab="Alter", ylab=bquote(alpha[x]), PlotLegend=c("ohne Varianznormierung", "mit Varianznormierung"), LegendPosition=c(10,-1.6), main="Parameter Alpha")
#ListPlot(list(lc.austria.f$bx, lc.austria.f.normalize$bx), xlab="Alter", ylab=bquote(beta[x]), PlotLegend=c("ohne Varianznormierung", "mit Varianznormierung"), LegendPosition=c(60,0.02), main="Parameter Beta")
#ListPlot(list(lc.austria.f$kt, lc.austria.f.normalize$kt), xlab="Alter", ylab=bquote(kappa[t]), PlotLegend=c("ohne Varianznormierung", "mit Varianznormierung"), LegendPosition=c(60,0.02), main="Parameter Kappa")
#par(mfrow=c(1,1))
#data=austria.f_fortgeschrieben[,1:96]
#data=austria.f_fortgeschrieben[,25:36]
#years=72:77
# full data:
#data=austria.f_fortgeschrieben[,1:96]
#years=47:112
.defaultcolors=c("blue", "red", "darkgreen", "green", "blueviolet")
.mathematicacolors=c("#3F3D99", "#993D71", "#998B3D", "#3D9956", "#3D5A99", "#993D90", "#996D3D", "#43993D", "#3D7999", "#843D99", "#994E3D", "#62993D", "#3D9799", "#653D99", "#993D4B")
.defaultcolors=.mathematicacolors;
extract.recycle = function (list, index) {
list[ (index-1) %% length(list) + 1]
}
recycle=function(list, len) {
rep(list, len/length(list)+1)[1:len]
}
LogListPlot = function(data, Joined=TRUE, PlotLegend=NULL, LegendPosition=NULL, col=NULL, ylim=NULL, xlim=NULL, xlab=NULL, ylab=NULL, main=NULL, ...) {
if (!is.list(data)) data=list(data);
if (is.null(col)) col=.defaultcolors;
colors=recycle(col, length(data));
# prepare the plot data, extract the ranges, etc.
datarange=range(unlist(data), na.rm = TRUE);
if (is.null(ylim)) {
ylim=c(10^floor(log10(datarange[1])), 10^ceiling(log10(datarange[2])));
} else if (identical(ylim,"tight")) {
ylim=datarange;
} else {
ylim;
}
plot(as.vector(data[[1]]), type="n", log="y", yaxt="n", ylim=ylim, xlim=xlim, xaxs="i", xlab=xlab, ylab=ylab, main=main, ...)
### Display all lists as lines
j=0;
for (i in data) {
j=j+1;
if (Joined) {
lines(i, col=extract.recycle(colors,j))
} else {
points(i, col=extract.recycle(colors,j))
}
}
### Plot legend from the given list, automatically use correct colors
if (!is.null(PlotLegend)) {
lpos = if (is.null(LegendPosition)) c(0, max(unlist(data), na.rm = TRUE)) else LegendPosition;
llabels=PlotLegend[1:min(length(data),length(PlotLegend),length(colors),na.rm = TRUE)]
legend(lpos[1],lpos[2], llabels, col=colors,lty=1, ...)
}
## Logarithmic ticks (large ticks at powers of 10), small ticks between
lticks=floor(log10(datarange[1])):ceiling(log10(datarange[2]))
axis(2, at=c(10^lticks), labels=sapply(lticks, function(i) as.expression(bquote(10^.(i)))),las=1,tcl=-0.5)
axis(2, at=c(outer(2:9, 10^(lticks))), tcl=-0.25, labels=FALSE)
}
ListPlot = function(data, Joined=TRUE, PlotLegend=NULL, LegendPosition=NULL, col=NULL, ylim=NULL, xlim=NULL, fill=NULL, ...) {
if (!is.list(data)) data=list(data);
if (is.null(col)) col=.defaultcolors;
colors=recycle(col, length(data));
# prepare the plot data, extract the ranges, etc.
datarange=range(unlist(data), na.rm = TRUE);
if (is.null(ylim)) {
# ylim=c(10^floor(log10(datarange[1])), 10^ceiling(log10(datarange[2])));
# } else if (ylim=="tight") {
ylim=datarange;
} else {
ylim;
}
plot(data[[1]], type="n", yaxt="n", ylim=ylim, xlim=xlim, xaxs="i", ...)
### If filling between lines is desired, create appropriate polygons before the lines:
if (!is.null(fill)) {
for (i in 1:(length(data))) {
if (is.vector(fill)) { clr = fill[i]; }
else { clr = fill; }
if (!is.null(clr) && !is.na(clr)) {
d1 = if (i==1) { rep(0, length(data[[i]])); } else { data[[i-1]]; }
d2 = data[[i]];
polygon(x=c(1:length(d1), length(d2):1), y=c(d1, rev(d2)), border=NA, col=clr);
}
}
}
### Display all lists as lines
j=0;
for (i in data) {
j=j+1;
if (Joined) {
lines(i, col=extract.recycle(colors,j))
} else {
points(i, col=extract.recycle(colors,j))
}
}
### Plot legend from the given list, automatically use correct colors
if (!is.null(PlotLegend)) {
lpos = if (is.null(LegendPosition)) c(0, max(unlist(data), na.rm = TRUE)) else LegendPosition;
llabels=PlotLegend[1:min(length(data),length(PlotLegend),length(colors),na.rm = TRUE)]
legend(lpos[1],lpos[2], llabels, col=colors,lty=1)
}
## ticks (large ticks at powers of 10), small ticks between
# lticks=floor(log10(datarange[1])):ceiling(log10(datarange[2]))
# axis(2, at=c(10^lticks), labels=sapply(lticks, function(i) as.expression(bquote(10^.(i)))),las=1,tcl=-0.5)
# axis(2, at=c(outer(2:9, 10^(lticks))), tcl=-0.25, labels=FALSE)
axis(2,las=1,tcl=-0.5)
}
MultiplePlot = function(data, PlotLegend=NULL, LegendPosition=NULL, col=NULL, ...) {
if (!is.list(data)) data=list(data);
if (is.null(col)) col=.defaultcolors;
colors=recycle(col, length(data));
# prepare the plot data, extract the ranges, etc.
# datarange=range(unlist(data), na.rm = TRUE);
# if (is.null(ylim)) {
# ylim=c(10^floor(log10(datarange[1])), 10^ceiling(log10(datarange[2])));
# } else if (ylim=="tight") {
# ylim=datarange;
# } else {
# ylim;
# }
plot(data[[1]], type="n", log="y", yaxt="n", xaxs="i", ...)
### Display all functions as lines
j=0;
for (i in data) {
j=j+1;
lines(i, col=extract.recycle(colors,j))
}
### Plot legend from the given list, automatically use correct colors
if (!is.null(PlotLegend)) {
lpos = if (is.null(LegendPosition)) c(0, 0) else LegendPosition;
llabels=PlotLegend[1:min(length(data),length(PlotLegend),length(colors),na.rm = TRUE)]
print(lpos)
legend(lpos[1],lpos[2], llabels, col=colors, lty=1, ...)
}
## Logarithmic ticks (large ticks at powers of 10), small ticks between
# lticks=floor(log10(datarange[1])):ceiling(log10(datarange[2]))
# axis(2, at=c(10^lticks), labels=sapply(lticks, function(i) as.expression(bquote(10^.(i)))),las=1,tcl=-0.5)
# axis(2, at=c(outer(2:9, 10^(lticks))), tcl=-0.25, labels=FALSE)
axis(2,las=1,tcl=-0.5)
axis(1)
}
# Usage:
# a=c(0.0001, 4, 187,0.3)
# b=c(983,1500,0.02,0.4)
# c=c(1,1,1)
# LogListPlot(list(a,b,c), PlotLegend=c("First", "Second", "Third", "Fourth"), LegendPosition=c(2, 0.1), main="LogListPlot example for multiple log-list-plots", ylab="test values", xlab="some index")
cumprod.matrix=function(x) {
y=matrix(1,nrow=dim(x)[1],ncol=dim(x)[2])
y[,1]=x[,1]
for (i in 2:dim(x)[2])
y[,i]=y[,i-1]*x[,i]
return(y)
}
cumsum.matrix=function(x) {
y=matrix(1,nrow=dim(x)[1],ncol=dim(x)[2])
y[,1]=x[,1]
for (i in 2:dim(x)[2])
y[,i]=y[,i-1]+x[,i]
return(y)
}
Version: 1.0
RestoreWorkspace: Default
SaveWorkspace: Default
AlwaysSaveHistory: Default
EnableCodeIndexing: Yes
UseSpacesForTab: Yes
NumSpacesForTab: 2
Encoding: UTF-8
RnwWeave: knitr
LaTeX: XeLaTeX
BuildType: Package
PackageInstallArgs: --no-multiarch --with-keep.source
* Edit the help file skeletons in 'man', possibly combining help files for multiple functions.
* Edit the exports in 'NAMESPACE', and add necessary imports.
* Put any C/C++/Fortran code in 'src'.
* If you have compiled code, add a useDynLib() directive to 'NAMESPACE'.
* Run R CMD build to build the package tarball.
* Run R CMD check to check the package tarball.
Read "Writing R Extensions" for more information.
\name{LogListPlot}
\alias{LogListPlot}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{Multiple logarithmic listplots in on e plot, with nice axis formatting
}
\description{
A function to print (multiple) list(s) in one logarithmic listplot. The axis ticks
are properly chosen to display all logarithmic subdivisions, a nice legend and all colors
are automatically chosen.
}
\usage{
LogListPlot(data, Joined = TRUE, PlotLegend = NULL,
LegendPosition = NULL, col = NULL, ylim = NULL,
xlim = NULL, ...)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
\item{data}{The list(s) to be plotted. For a single list-plot this can be a vector, for a multiple list plot this should be an array of the vectors to be plotted.
}
\item{Joined}{Whether the plot data is joined by lines (type="l") or displayed as separate points (type="p").}
\item{PlotLegend}{A list of legend labels. The corresponding lines and their colors will be automatically chosen from the plot colors.}
\item{LegendPosition}{2-element vector giving the position of the plot legend (x- and y-values).}
\item{col}{Vector of colors of the plot lines and legend labels.}
\item{ylim}{same as with plot}
\item{xlim}{same as with plot}
\item{\dots}{arguments passed on to plot}
}
\details{
The syntax of this call is chosen to match the Mathematica function ListLogPlot.
}
%\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{
Reinhold Kainhofer <reinhold@open-tools.net>
}
%\note{
%% ~~further notes~~
%}
%% ~Make other sections like Warning with \section{Warning }{....} ~
%\seealso{
%% ~~objects to See Also as \code{\link{help}}, ~~~
%}
\examples{
a=c(0.0001, 4, 187,0.3)
b=c(983,1500,0.02,0.4)
c=c(1,1,1)
LogListPlot(list(a,b,c), PlotLegend=c("First", "Second", "Third", "Fourth"), LegendPosition=c(2, 0.1), main="LogListPlot example for multiple log-list-plots", ylab="test values", xlab="some index")
}
% Add one or more standard keywords, see file 'KEYWORDS' in the
% R documentation directory.
\keyword{ ~hplot }
\keyword{ ~utilities }
\name{RKTools-package}
\alias{RKTools-package}
\alias{RKTools}
\docType{package}
\title{
Some utility functions (LogListPlot, etc.) by R. Kainhofer
}
\description{
Utility functions:
-) LogListPlot: (multiple) log list-plots with nice axis formatting, legend and automatic colors
}
\details{
\tabular{ll}{
Package: \tab RKTools\cr
Type: \tab Package\cr
Version: \tab 1.0\cr
Date: \tab 2013-04-29\cr
License: \tab GPL (>=3)\cr
}
~~ An overview of how to use the package, including the most important functions ~~
}
\author{
Reinhold Kainhofer <office@open-tools.net>
Maintainer: Reinhold Kainhofer <office@open-tools.net>
}
% \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 ~~
% }
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment