[R] numerical approximation of survival function
Dear R-users,
my question concerns numerical approximation and, somehow, survival
analysis.
Let’s assume we have a density function and we aim to numerically
compute the hazard, which is, in theory, the ratio between density and
survival.
In the following example, I take a pdf from a log-normal and proceed as
I would have known only this function. I numerically approximate the
survival function by splinefun() and integrate().Finally I compute the
hazard. It’s easy to see that (1) true hazard is different from the
numerically approximated one and (2) discrepancy depends upon the
approximation carried out to compute the survival. Is there any way to
overcome or attenuate this issue? Or do I just miss something obvious?
In my real study I do not deal with analytically defined pdf, so I would
like to be sure that approximation is done at its best.
Thanks in advance for your help.
Giancarlo Camarda
## x-values
delta <- 0.1
x <- seq(0, 200, by=delta)
m <- length(x)
## true pdf
fxT <- dlnorm(x, meanlog=3.5, sdlog=0.5)
fxT <- fxT/sum(fxT)
## true survival
SxT <- 1 - plnorm(x, meanlog=3.5, sdlog=0.5)
## true hazard
hxT <- fxT/SxT
## numerical approximation of the pdf
fx <- spline(x, fxT, n=m)$y
## numerical approximation of the survival
Sx <- numeric(m)
fxfun <- splinefun(x, fx/delta)
for(j in 1:m){
Sx[j] <- 1 - integrate(fxfun, x[1], x[j])$value
}
## hazard: pdf divided by survival
hx <- fx/Sx ## both numerical approx
hx1 <- fx/SxT ## only the survival numerically approx
hx2 <- fxT/Sx ## only the pdf numerically approx
## plotting hazards
plot(x, hxT)
lines(x, hx, col=2, lwd=2)
lines(x, hx1, col=3, lwd=2)
lines(x, hx2, col=4, lty=2, lwd=2)
## it seems that we have got something
## at the survival approx level:
plot(x, Sx/SxT)
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] matrix built by diagonal matrices with a given structure (2nd trial)
Dear R-users, three weeks ago I sent the mail below, but I didn't receive any response. Maybe it was overlooked. Thanks anyway for all the help you gave us by this mailing-list, Giancarlo Camarda [...] I have a matrix with a series of not-overlapping in a row dimension vectors in a given structure. Something like: |a1, 0, 0, 0| | 0, a2, a3, 0| |a4, 0, 0, a5| where ai are column-vectors of the equal length, m. My aim is to construct a new matrix formed by diagonal matrices built by the mentioned vectors and placed following the original structure. Something like: |diag(a1),0,0,0| | 0, diag(a2), diag(a3),0| |diag(a4),0,0, diag(a5)| Of course the zeros are vectors of length m, and empty (m times m) matrices in the first and second scheme, respectively. I found a way to obtain what I need by selecting an augmented version of the original matrix which I have constructed using the kronecker product. I was wondering whether there is a more elegant and straightforward procedure. See below a simple reproducible example of my challenge in which the length of the vectors is 4. Thanks in advance for your help, Giancarlo Camarda ## size of the diagonal matrices ## or length of the vectors m - 4 ## the original matrix ze - rep(0,m) A - cbind(c(1,2,3,4,ze,13,14,15,16), c(ze,5,6,7,8,ze), c(ze,9,10,11,12,ze), c(ze,ze,17,18,19,20)) ## augmenting the original matrix A1 - kronecker(A, diag(m)) ## which rows to select w1 - seq(1, m^2, length=m) w2 - seq(0, 2*m^2, by=m^2) w0 - outer(w1, w2, FUN=+) w - c(w0) ## final matrix A2 - A1[w,] [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Simplifying matrix computation
Dear David, and R-users, thanks for the response. I'll do my best to describe the context. My data consists in two matrices, D and T. The former provides simple values, the latter is a Boolean matrix. The number of rows in D is equal to the number of columns in T. My aim is to construct a larger matrix which is made up of sum(T) diagonal matrices constructed from the rows of D, and placed according to a structure given by a transformation of T. Whereas the positions of the columns in which T is equal to 1 informs about the placement of the diagonal matrices, the positions of the rows in which T is equal to 1 informs about which row to take from D in order to create these diagonal matrices. In the code below (corrected in the last line with respect to the previous version, sorry), I first created the matrices D and T. Then I built a 0/1 matrix (M0) that presents the structure of the final block-matrix made up of diagonal matrices, and, by Kronecker product, a matrix with the final 0/1 structure is constructed, M. Using the rows in which T is equal to 1, I select the values from D and place in the matrix M, where such matrix is equal to 1. It's a long explanation and I could not find a better way to program (and explain) it. Thanks for you help, Giancarlo ## data in matrices D - matrix(1:15, 3, 5) T - matrix(0, 3, 3) T[c(2,4,6,8)] - 1 ## the col of T equal to 1 gives the position wr - which(T==1, arr.ind=TRUE)[,2] ## we aim to sum(T) diagonal matrices wc - 1:sum(T) ## structure: how to place the diagonal matrices M0 - matrix(0, nrow(T), sum(T)) M0[cbind(wr,wc)] - 1 ## number of columns m - ncol(D) ## final 0/1 matrix M - kronecker(M0, diag(m)) ## the row of T equal to 1 gives which rows to take from D pos - which(T==1, arr.ind=TRUE)[,1] ## filling up with data M[M!=0] - t(D[pos,]) On 29/01/2014 01:49, David Winsemius wrote: On Jan 27, 2014, at 8:04 AM, Carlo Giovanni Camarda wrote: Dear R-users, I would like to know whether you know some trick for skipping some of the steps in the example below (especially the last step in a way that would make easier to be written succinctly in a text). I could try to explain in words the whole process, but I'm sure the code below would be clearer. After looking at the code and output, I must disagree. The lack of any response from the rest of the readership suggests to me that I am not the only one who thinks a natural language description of the context and goals for this effort would help. Thanks in advance for your help, Giancarlo ## data in matrices D - matrix(1:15, 3, 5) T - matrix(0, 3, 3) T[c(2,4,6,8)] - 1 ## how to place the diag matrices of each row M0 - matrix(0, nrow(T), sum(T)) wr - which(T==1, arr.ind=TRUE)[,2] wc - 1:ncol(M0) M0[cbind(wr,wc)] - 1 ## number of columns m - ncol(D) ## non-zero positions M - kronecker(M0, diag(m)) ## which rows to take pos - which(T==1, arr.ind=TRUE)[,1] ## filling up with data M[M!=0] - t(D[wr,]) [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius Alameda, CA, USA [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Simplifying matrix computation
Dear R-users, I would like to know whether you know some trick for skipping some of the steps in the example below (especially the last step in a way that would make easier to be written succinctly in a text). I could try to explain in words the whole process, but I'm sure the code below would be clearer. Thanks in advance for your help, Giancarlo ## data in matrices D - matrix(1:15, 3, 5) T - matrix(0, 3, 3) T[c(2,4,6,8)] - 1 ## how to place the diag matrices of each row M0 - matrix(0, nrow(T), sum(T)) wr - which(T==1, arr.ind=TRUE)[,2] wc - 1:ncol(M0) M0[cbind(wr,wc)] - 1 ## number of columns m - ncol(D) ## non-zero positions M - kronecker(M0, diag(m)) ## which rows to take pos - which(T==1, arr.ind=TRUE)[,1] ## filling up with data M[M!=0] - t(D[wr,]) [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
