Re: [R] hessian in solnp
Thanks. Yes, there is a Lagrange multiplier (though really small for the
example). However, AFAIK the hessian is the square matrix of second-order
partial derivatives and, given the simple constraint in the example, the second
(partial) derivatives of the lagrangian function should simplify to the second
(partial) derivatives of the loss function. Or do I miss something? Moreover, I
would expect smaller standard errors of the estimated parameters when
additional constraints are included, and it is not the case in the given
example.
De: "Ivan Krylov"
À: "CAMARDA Carlo Giovanni via R-help"
Cc: "Carlo Giovanni Camarda"
Envoyé: Jeudi 28 Avril 2022 06:30:00
Objet: Re: [R] hessian in solnp
В Thu, 28 Apr 2022 00:16:05 +0200 (CEST)
CAMARDA Carlo Giovanni via R-help пишет:
> when a constraint is added, hessian matrix is obviously changing, but
> in a way I don't understand.
Isn't it the point of augmented Lagrange multiplier, to solve the
constrained optimisation problem by modifying the loss function and
optimising the result in an unconstrained manner? Apologies if I
misunderstood your question.
Starting on
<https://github.com/cran/Rsolnp/blob/4b56bb5cd7c5d1096d1ba2f3946df7afa9af4201/R/subnp.R#L282>,
we can see how the functions are called: both the loss and the
constraint function are concatenated into the `obm` vector (if there's
no constraint, the function returns NULL, which is eaten by
concatenation), which form the vectors `g` (seems to be the gradient)
and `p`, which, in turn, form the matrix `hessv`.
My reading of the code could be wrong (and so could be my understanding
of augmented Lagrangian methods). Contacting maintainer('Rsolnp') could
be an option; maybe there's some documentation for the original MATLAB
version of the code at <https://web.stanford.edu/~yyye/Col.html>?
--
Best regards,
Ivan
[[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] hessian in solnp
In the commented snippet below, a behavior of the solnp() I cannot fully
explain: when a constraint is added, hessian matrix is obviously changing, but
in a way I don't understand. I know that the model below could be estimated
differently and in a much efficient way. However I wanted to present the
simplest example I could think.
In particular, I encountered this issue because I aim to analytically compute
variance-covariance matrix in a (non-linear) model with equality constraints.
Any practical hint to find a solution to this more general issue would be also
welcome.
Thanks, GC
## simulating a simple linear model
n <- 100
x <- sort(runif(n))
X <- cbind(1,x)
k <- ncol(X)
y <- X%*%c(3,4) + rnorm(n, sd=0.5)
## estimating with lm()
fitlm <- lm(y~x)
## using solnp
library(Rsolnp)
## objecting function
fn1 <- function(par){
y.hat <- X%*%par
res <- y-y.hat
RSS <- t(res)%*%res/2
return(c(RSS))
}
## using solnp without any constraint
solUR <- solnp(c(3.2, 3.5), fun = fn1)
## getting the same fitted objects as in lm()
## coefficients
cbind(solUR$pars, fitlm$coef)
## and standard errors by variance-covariance matrix computed
## with hessian internally estimated
(H.UR <- solUR$hessian)
y.UR <- X%*%solUR$pars
RSS.UR <- t(y-y.UR)%*%(y-y.UR)
sigma2.UR <- RSS.UR/(n-k)
V.UR <- c(sigma2.UR)*solve(H.UR)
se.UR <- sqrt(diag(V.UR))
cbind(se.UR,summary(fitlm)$coef[,2])
## adding equality constraint
eqn1 <- function(par) sum(par)
c <- sum(solUR$pars)
## I use the sum of the previously estimated coefficients
solR <- solnp(c(3.2, 3.5), fun = fn1, eqfun = eqn1, eqB = c)
## as expected we get the same coefficients
cbind(solR$pars,fitlm$coef)
## but the hessian is rather different
(H.R <- solR$hessian)
## which leads to completely different (much larger here)
## standard errors
y.R <- X%*%solR$pars
RSS.R <- t(y-y.R)%*%(y-y.R)
sigma2.R <- RSS.R/(n-k)
V.R <- c(sigma2.R)*solve(H.R)
se.R <- sqrt(diag(V.R))
cbind(se.R,summary(fitlm)$coef[,2])
## I noticed that differences in the hessian
## is ~constant over both dimensions
## and it depends (also) on sample size
H.R-H.UR
[[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] Reproducing lm outcomes with mixed effect model on augmented dataset
Dear R-users, let's assume a simple linear model in which to each original observation one attaches random response and covariate, but with weights that are practically zero. If one runs a linear model to both original and augmented data, one obtains the same estimates, but obviously different standard errors. Could we reproduce the lm() outcomes in a mixed effect model setting using, for instance, lme(nlme)? Is there any way to get the lm() estimates and standard deviations by accounting for correlation within each new pseudo-observation? Below a reproducible and commented instance of the problem . Thanks in advance for your help, GC library(nlme) ## simulating data m <- 1000 x <- runif(m) y <- 3 + 4*x + rnorm(m, sd=0.5) id <- 1:m dati <- cbind(id=id, x=x, y=y, w=1) ## new data: repeat n times original data n <- 10 datiA <- kronecker(dati, matrix(1,n,1)) ## trasform in dataframes with suitable names dati <- as.data.frame(dati) datiA <- as.data.frame(datiA) names(datiA) <- names(dati) ## assign ~0 weight to everyone, but the original observations fw <- 10^-30 datiA$w <- fw datiA$w[seq(1,m*n,n)] <- 1 - n*fw ## assign random numbers to all response/covariate, but for the original observations datiA$x[datiA$w==fw] <- runif(m*n-m) datiA$y[datiA$w==fw] <- runif(m*n-m) ## fits with lm() fit <- lm(y~x, weights=w, data=dati) summary(fit)$coef[,1:2] fitA <- lm(y~x, weights=w, data=datiA) summary(fitA)$coef[,1:2] ## mixed model approach, lme() ## attempting to reproduce results from "fit" fitMM <- lme(y ~ x, weights=varFixed(~1/w), random=~1|id, data=datiA) summary(fitMM)$tTable[,1:2] [[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
Dear R-users, 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] extracting the diagonal of an inverse matrix
Thanks: really nice idea.
Unfortunately the solution you suggest can not be adopted with sparse matrices:
solve command with logical is not yet implemented. (Or do I have a too-old
R-version?)
## loading package
library(Matrix)
## transform A in sparse matrix
Asp <- Diagonal(9)
pos <- which(A!=0)
Asp[pos] <- A[pos]
## attempt the solution
sapply(seq_len(nrow(Asp)), function(i)solve(Asp, seq_len(nrow(Asp))==i)[i])
## Error: not-yet-implemented method for solve(, ).
## ->> Ask the package authors to implement the missing feature.
> version
_
platform i686-pc-linux-gnu
arch i686
os linux-gnu
system i686, linux-gnu
status
major 2
minor 14.1
year 2011
month 12
day22
svn rev57956
language R
version.string R version 2.14.1 (2011-12-22)
From: William Dunlap [wdun...@tibco.com]
Sent: Monday, April 22, 2013 6:25 PM
To: Camarda, Carlo Giovanni; Greg Snow
Cc: r-h...@stat.math.ethz.ch
Subject: RE: [R] extracting the diagonal of an inverse matrix
You could save the space required to store the inverse of A by repeating
solve(A,column)
for each column of an identity matrix the size of A, and store just the
relevant element
of the result. E.g.,
> diag(solve(A))
[1] 0.39576620 0.12064108 0.20705171 0.10067892 0.03191425 0.05857497
0.10749426 0.03868911 0.08550176
> sapply(seq_len(nrow(A)), function(i)solve(A, seq_len(nrow(A))==i)[i])
[1] 0.39576620 0.12064108 0.20705171 0.10067892 0.03191425 0.05857497
0.10749426 0.03868911 0.08550176
It may take more time that you would like.
Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
> Behalf
> Of Camarda, Carlo Giovanni
> Sent: Monday, April 22, 2013 1:51 AM
> To: Greg Snow
> Cc: r-h...@stat.math.ethz.ch
> Subject: Re: [R] extracting the diagonal of an inverse matrix
>
> Dear Greg,
>
> thanks a lot for your prompt reply. I was aware of the link and the
> associated papers.
> In my previous mail I was wondering whether there was already an R-routine
> for coping
> with this issue.
>
> I add below a testable example which reproduces my problem.
>
> Thanks again,
> Giancarlo
>
> ## dimensions (my actual dimensions are much larger and I use sparse matrices)
> m <- 3
> n <- 3
> ## building A
> ## diagonal part
> d <- 1:(m*n)
> D <- diag(d)
> ## additional ~dense part (:penalty term)
> P1 <- diff(diag(m), diff=2)
> P2 <- diff(diag(n), diff=2)
> P1a <- kronecker(diag(n), t(P1)%*%P1)
> P2a <- kronecker(t(P2)%*%P2, diag(m))
> P <- 1000 * P1a + 1000 * P2a
> ## final matrix A
> A <- D + P
> ## solving A
> B <- solve(A)
> ## what I just actually need
> diag(B)
>
>
>
> From: Greg Snow [538...@gmail.com]
> Sent: Saturday, April 20, 2013 12:16 AM
> To: Camarda, Carlo Giovanni
> Cc: r-h...@stat.math.ethz.ch
> Subject: Re: [R] extracting the diagonal of an inverse matrix
>
> This link
> http://math.stackexchange.com/questions/18488/diagonal-of-an-inverse-of-a-
> sparse-matrix might help. It is about sparse matrices, but the general idea
> should be able
> to be extended to non-sparse matrices as well.
>
>
> On Fri, Apr 19, 2013 at 8:13 AM, Camarda, Carlo Giovanni
> mailto:cama...@demogr.mpg.de>> wrote:
> Dear R-users,
>
> I would like to know whether there is a way to extract a diagonal of an
> inverse matrix
> without computing the inverse of the matrix itself. The size of my matrices
> are really
> huge and, also using sparse matrix, computing the inverse leads to storage
> problems and
> low speed.
>
> In other words, given a square matrix A, I aim to know diag(B), where
> B=solve(A),
> without computing solve(A).
>
> Accidentally (I do not know whether it helps), I could write the matrix A as
> follows:
> A <- D + P
> where D is a diagonal matrix.
>
> I read there are methods around, but, before implementing one of them by
> myself, could
> you please inform whether there is already an R-routine for this issue?
>
> Thanks in advance for the help you could provide,
> Carlo Giovanni Camarda
> --
> This mail has been sent through the MPI for Demographic ...{{dropped:19}}
>
>
Re: [R] extracting the diagonal of an inverse matrix
Dear Greg,
thanks a lot for your prompt reply. I was aware of the link and the associated
papers.
In my previous mail I was wondering whether there was already an R-routine for
coping with this issue.
I add below a testable example which reproduces my problem.
Thanks again,
Giancarlo
## dimensions (my actual dimensions are much larger and I use sparse matrices)
m <- 3
n <- 3
## building A
## diagonal part
d <- 1:(m*n)
D <- diag(d)
## additional ~dense part (:penalty term)
P1 <- diff(diag(m), diff=2)
P2 <- diff(diag(n), diff=2)
P1a <- kronecker(diag(n), t(P1)%*%P1)
P2a <- kronecker(t(P2)%*%P2, diag(m))
P <- 1000 * P1a + 1000 * P2a
## final matrix A
A <- D + P
## solving A
B <- solve(A)
## what I just actually need
diag(B)
From: Greg Snow [538...@gmail.com]
Sent: Saturday, April 20, 2013 12:16 AM
To: Camarda, Carlo Giovanni
Cc: r-h...@stat.math.ethz.ch
Subject: Re: [R] extracting the diagonal of an inverse matrix
This link
http://math.stackexchange.com/questions/18488/diagonal-of-an-inverse-of-a-sparse-matrix
might help. It is about sparse matrices, but the general idea should be able
to be extended to non-sparse matrices as well.
On Fri, Apr 19, 2013 at 8:13 AM, Camarda, Carlo Giovanni
mailto:cama...@demogr.mpg.de>> wrote:
Dear R-users,
I would like to know whether there is a way to extract a diagonal of an inverse
matrix without computing the inverse of the matrix itself. The size of my
matrices are really huge and, also using sparse matrix, computing the inverse
leads to storage problems and low speed.
In other words, given a square matrix A, I aim to know diag(B), where
B=solve(A), without computing solve(A).
Accidentally (I do not know whether it helps), I could write the matrix A as
follows:
A <- D + P
where D is a diagonal matrix.
I read there are methods around, but, before implementing one of them by
myself, could you please inform whether there is already an R-routine for this
issue?
Thanks in advance for the help you could provide,
Carlo Giovanni Camarda
--
This mail has been sent through the MPI for Demographic ...{{dropped:19}}
__
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] extracting the diagonal of an inverse matrix
Dear R-users, I would like to know whether there is a way to extract a diagonal of an inverse matrix without computing the inverse of the matrix itself. The size of my matrices are really huge and, also using sparse matrix, computing the inverse leads to storage problems and low speed. In other words, given a square matrix A, I aim to know diag(B), where B=solve(A), without computing solve(A). Accidentally (I do not know whether it helps), I could write the matrix A as follows: A <- D + P where D is a diagonal matrix. I read there are methods around, but, before implementing one of them by myself, could you please inform whether there is already an R-routine for this issue? Thanks in advance for the help you could provide, Carlo Giovanni Camarda -- This mail has been sent through the MPI for Demographic Research. Should you receive a mail that is apparently from a MPI user without this text displayed, then the address has most likely been faked. If you are uncertain about the validity of this message, please check the mail header or ask your system administrator for assistance. __ 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] linear model with equality and inequality (redundant) constraints
Dear R-users, in the last days I have been trying to estimate a normal linear model with equality and inequality constraints. Please find below a simple example of my problem. Of course, one could easily see that, though the constraints are consistent, there is some redundancy in the specific constraints. Nevertheless my actual applications can get much larger and I would not like to manually check for presence of redundant constraints in the setting. First question: would it be possible to estimate a linear model with equality and inequality redundant constraints? I also understand that quadratic programming is the key for that, but I was not able to “translate” my simple linear problem with constraints for using solve.QP(quadprog) and ic.est(ic.infer). I have found pcls(mgcv) more user-friendly. And below you could find my attempt to fit the toy-data with this function with and without redundant constraints. Thanks in advance for any help you could provide, Carlo Giovanni Camarda rm(list = ls()) library(mgcv) ## estimate p: E(r) = X p, s.t. ## C p == c ## A p >= b (though pcls allows only for A p > b) ## response r <- c(4277.1265,57004.1177,7135.9204,541.6218,1455.2135) ## design matrix x <- c(6029,0,0,0,0,0,6029,0, 0,0,0,44453,0,0,0,0, 1284,0,0,0,0,3042,0, 0,0,0,10132,0,0,0,0,0, 10132,0,0,0,0,1955,0, 0,0,0,3519,0,0,0,0,0, 10132,0,0,0,0,0,44453, 0,0,0,0,1955) X <- matrix(x, 5, 12) ## equality constraints matrix C <- matrix(0,7,12) C[1,1:2] <- 1 C[2,c(3,11)] <- 1 C[3,4] <- 1 C[4,5] <- 1 C[5,c(6,7,10)] <- 1 C[6,c(8,12)] <- 1 C[7,9] <- 1 ## equality constraints vector, which is not used in pcls c <- rep(1,7) ## inequality constraints matrix: p \in [0,1] A <- rbind(diag(12), -1*diag(12)) ## inequality constraints vector b <- c(rep(0,12), rep(-1,12)) ## starting values p.st <- c(0.5,0.5,0.5,1,1,0.2,0.3,0.5,1,0.5,0.5,0.5) ## test starting values table(C%*%p.st == c) table(A%*%p.st > b) ## actually I would like to have A p >= b ## list for pcls M <- list(X=X, p=p.st, off=array(0,0), S=list(), Ain=A, bin=b, C=C, sp=array(0,0), y=r, w=r*0+1) ## estimation (p.hat <- pcls(M)) avoid redundant constraints ## response r1 <- c(4277.1265,46649.1177,3616.9204,-9590.3782,-44952.7865) ## design matrix X1 <- c(6029,-6029,0,0,0,0,44453, 0,0,-44453,0,10132,0,-10132, 0,0,0,10132,-10132,0,0,0, 1955,0,-1955) dim(X1) <- c(5,5) ## inequality constraints matrix A1 <- rbind(diag(5), -1*diag(5), c(0,0,-1,-1,0)) ## inequality constraints vector b1 <- c(rep(0,5), rep(-1,5), -1) p.st1 <- runif(5, 0.1, 0.2) ## check constraints A1%*%p.st1 > b1 M1 <- list(X=X1, p=p.st1, off=array(0,0), S=list(), Ain=A1, bin=b1, C=matrix(0,0,0), sp=array(0,0), y=r1, w=r1*0+1) (p.hat1 <- pcls(M1)) A1%*%p.hat1 > b1 -- This mail has been sent through the MPI for Demographic Research. Should you receive a mail that is apparently from a MPI user without this text displayed, then the address has most likely been faked. If you are uncertain about the validity of this message, please check the mail header or ask your system administrator for assistance. __ 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] impossible to invert a spam-object, but possible when it's a matrix-object
Dear Uwe,
thanks for the response.
I knew my matrix was almost singular, but is there any way to implement the
algorithm is solve(base) with a spam-matrix?
As workaround, I might add a ridge penalty:
ApP <- A+10^-1*diag.spam(nrow(A))
solve(ApP)
but this would completely jeopardize other parts of my model.
Thanks again,
Giancarlo
From: Uwe Ligges [lig...@statistik.tu-dortmund.de]
Sent: Saturday, February 09, 2013 8:04 PM
To: Camarda, Carlo Giovanni
Cc: r-h...@stat.math.ethz.ch
Subject: Re: [R] impossible to invert a spam-object, but possible when it's a
matrix-object
On 05.02.2013 13:29, Camarda, Carlo Giovanni wrote:
> Dear R-users,
>
> a question concerning sparse matrices in package "spam" (spam_0.29-2).
>
> On one hand I have a spam object (n X n) from which I cannot compute the
> inverse. On the other hand, if I convert this object in a plain matrix, I can
> find the inverse without any problem.
>
> Specifically I get the following error message:
>Error in chol.spam(a, ...) :
>Singularity problem when calculating the Cholesky factor.
>
> Obviously I get similar behaviour when I compute determinants of these
> objects.
>
> Please see below a toy-example which I created based on my actual problem.
>
Different algorithms are used to compute the inverse - and you are
rather close to singularity: kappa(A) results in 59638727.
Best,
Uwe Ligges
> Thanks in advance for any help,
> GC
>
>
> library(spam)
> ## creating a spam matrix A
> ent <- c(2312.12324929972,-2000,1000,-2000,1000,-2000,
> 5031.91011235955,-2000,-2000,1000,-1,1000,-2000,
> 2049.8595505618,-2000,1000,-1,-2000,5036.89635854342,
> -2000,1000,-2000,1,-2000,-2000,8119.66292134831,-2000,
> -2000,-2000,1000,-2000,5058.83426966292,-2000,-1,1000,
> -2000,2051.85434173669,-2000,1000,1,1000,-2000,-2000,
> 5043.87640449438,-2000,1,1000,-2000,1000,-2000,
> 2110.68820224719,-1,1,0,-1,1,-1,1)
> rr <- as.integer(c(1,6,12,18,24,29,35,41,47,52,55,57,59))
> cc <- as.integer(c(1,2,3,4,7,1,2,3,5,8,10,1,2,3,6,9,11,1,
> 4,5,6,7,10,2,4,5,6,8,3,4,5,6,9,12,1,4,
> 7,8,9,11,2,5,7,8,9,12,3,6,7,8,9,2,4,10,3,7,6,8))
> di <- as.integer(c(12,12))
> A <- new("spam",
> entries=ent,
> colindices=cc,
> rowpointers=rr,
> dimension=di)
> ## calculate the determinant
> det(A)
> ## computes the inverse
> solve(A)
>
> ## transform A in a plain matrix
> A1 <- matrix(c(A), di)
> ## calculate the determinant
> det(A1)
> ## computes the inverse
> solve(A1)
>
>
>> version
> _
> platform i686-pc-linux-gnu
> arch i686
> os linux-gnu
> system i686, linux-gnu
> status
> major 2
> minor 14.1
> year 2011
> month 12
> day22
> svn rev57956
> language R
> version.string R version 2.14.1 (2011-12-22)
>
> --
> This mail has been sent through the MPI for Demographic Research. Should you
> receive a mail that is apparently from a MPI user without this text
> displayed, then the address has most likely been faked. If you are uncertain
> about the validity of this message, please check the mail header or ask your
> system administrator for assistance.
>
> __
> 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-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] impossible to invert a spam-object, but possible when it's a matrix-object
Dear R-users,
a question concerning sparse matrices in package "spam" (spam_0.29-2).
On one hand I have a spam object (n X n) from which I cannot compute the
inverse. On the other hand, if I convert this object in a plain matrix, I can
find the inverse without any problem.
Specifically I get the following error message:
Error in chol.spam(a, ...) :
Singularity problem when calculating the Cholesky factor.
Obviously I get similar behaviour when I compute determinants of these objects.
Please see below a toy-example which I created based on my actual problem.
Thanks in advance for any help,
GC
library(spam)
## creating a spam matrix A
ent <- c(2312.12324929972,-2000,1000,-2000,1000,-2000,
5031.91011235955,-2000,-2000,1000,-1,1000,-2000,
2049.8595505618,-2000,1000,-1,-2000,5036.89635854342,
-2000,1000,-2000,1,-2000,-2000,8119.66292134831,-2000,
-2000,-2000,1000,-2000,5058.83426966292,-2000,-1,1000,
-2000,2051.85434173669,-2000,1000,1,1000,-2000,-2000,
5043.87640449438,-2000,1,1000,-2000,1000,-2000,
2110.68820224719,-1,1,0,-1,1,-1,1)
rr <- as.integer(c(1,6,12,18,24,29,35,41,47,52,55,57,59))
cc <- as.integer(c(1,2,3,4,7,1,2,3,5,8,10,1,2,3,6,9,11,1,
4,5,6,7,10,2,4,5,6,8,3,4,5,6,9,12,1,4,
7,8,9,11,2,5,7,8,9,12,3,6,7,8,9,2,4,10,3,7,6,8))
di <- as.integer(c(12,12))
A <- new("spam",
entries=ent,
colindices=cc,
rowpointers=rr,
dimension=di)
## calculate the determinant
det(A)
## computes the inverse
solve(A)
## transform A in a plain matrix
A1 <- matrix(c(A), di)
## calculate the determinant
det(A1)
## computes the inverse
solve(A1)
> version
_
platform i686-pc-linux-gnu
arch i686
os linux-gnu
system i686, linux-gnu
status
major 2
minor 14.1
year 2011
month 12
day22
svn rev57956
language R
version.string R version 2.14.1 (2011-12-22)
--
This mail has been sent through the MPI for Demographic Research. Should you
receive a mail that is apparently from a MPI user without this text displayed,
then the address has most likely been faked. If you are uncertain about the
validity of this message, please check the mail header or ask your system
administrator for assistance.
__
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] matrix algebra for constructing a special matrix?
Dear R-users,
would it be a better way to construct the matrix below without using any
for-loop or model.matrix? preferably with some matrix algebra?
Thanks in advance,
Carlo Giovanni Camarda
## dimensions
m <- 3
n <- 4
mn <- m*n
k <- m+n-1
## with a for-loop
X <- matrix(0, mn, k)
for(i in 1:n){
wr <- 1:m+(i-1)*m
wc <- rev(1:m+(i-1))
where <- cbind(wr,wc)
X[where] <- 1
}
## using model.matrix
ff <- factor(c(outer(m:1, seq(0,m), "+")))
X1 <- matrix(model.matrix(~-1+ff), mn, k)
--
This mail has been sent through the MPI for Demographic ...{{dropped:10}}
__
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] xyplot with different pch and col in each panel and additional line
Dear R-users,
I'd like to use an xyplot(lattice) in which in each panel I have points with
different point-character and color, and additional lines with the same color.
Please find below a toy example in which I did not manage to change such
parameters, and the associated basic plot() in which I show what I aim to
within lattice.
Thanks for any help you could provide,
Giancarlo Camarda
library(lattice)
y <- seq(0,1, length=15)
yy <- y + rnorm(15,sd=0.05)
dataset <- data.frame(time = rep(1:5,times=3),
y = y,
yy = yy,
type = rep(LETTERS[1:3], each=5))
xyplot(yy+y ~ time | type, dataset,
layout=c(3,1),
type=c("p", "l"),
col=c(1,2),
panel = function(...) {
panel.superpose.2(...)
})
par(mfrow=c(1,3))
plot(1:5, yy[1:5], pch=1, col=3, ylim=c(0,1))
lines(1:5, y[1:5], col=2)
plot(1:5, yy[1:5+5], pch=2, col=4, ylim=c(0,1))
lines(1:5, y[1:5+5], col=2)
plot(1:5, yy[1:5+10], pch=3, col=5, ylim=c(0,1))
lines(1:5, y[1:5+10], col=2)
--
This mail has been sent through the MPI for Demographic ...{{dropped:10}}
__
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] generalized additive mixed models for ordinal data
Dear R-users,
Is there any R-function for fitting generalized additive mixed
models for ordinal data? Do they actually make some sense? I can fit a
generalized linear mixed model for ordinal data using the function
clmm(ordinal) and I'm able to cope with generalized additive model for
ordinal data within the package VGAM.
But I would like to fit something like:
g(\gamma_{ij}) = \theta_{j} + x_{i1} \beta_1 + f(x_{2i}) + u_{i},
where \gamma_{ij} denote the cumulative probability that the i-th
observation falls in the j-th category or below.
Sorry for the rather out-of-R question,
Carlo Giovanni Camarda
--
This mail has been sent through the MPI for Demographic ...{{dropped:10}}
__
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] persp3d and rgl.viewpoint for rotating 3D plots
Dear Duncan,
thanks for your prompt reply.
I was tempted to use movie3d, but below I copied what I get on my PC. Then, I
have ImageMagick installed and, though I set convert equal to FALSE, I was not
able to find where the function writes the .png files.
Maybe I lack experience in managing such complex functions, but I must confess
that I have quite some hard time to grasp info from help(movie3d).
Thanks in advance for any help,
Carlo Giovanni
> library(rgl)
> # set data
> f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
> x <- seq(-10, 10, length= 30); y <- x
> z <- outer(x, y, f)
> # set working directory
> setwd("U:\\movie")
> dir() # check
character(0)
> # plot
> persp3d(x, y, z)
> movie3d(spin3d(), duration=5, convert=TRUE)
Writing movie000.png
[...]
Writing movie050.png
Error in movie3d(spin3d(), duration = 5, convert = TRUE) :
ImageMagick not found
> dir() # check again
character(0)
> movie3d(spin3d(), duration=5, convert=FALSE)
Writing movie000.png
[...]
Writing movie050.png
> dir() # check again ??
character(0)
>
> version
_
platform i386-pc-mingw32
arch i386
os mingw32
system i386, mingw32
status
major 2
minor 8.1
year 2008
month 12
day22
svn rev47281
language R
version.string R version 2.8.1 (2008-12-22)
-Original Message-
From: Duncan Murdoch [mailto:murd...@stats.uwo.ca]
Sent: Wednesday, April 08, 2009 5:13 PM
To: Camarda, Carlo Giovanni
Cc: r-h...@stat.math.ethz.ch
Subject: Re: [R] persp3d and rgl.viewpoint for rotating 3D plots
On 4/8/2009 11:01 AM, Camarda, Carlo Giovanni wrote:
> Dear R-users,
>
> within the rgl-package, I would have a question about the usage of
> persp3d in combination of rgl.viewpoint.
> I am not able to figure out how to let a 3D plot rotating around likewise the
> example in ?rgl.viewpoint. It seems that when I use persp3d(...) I see
> something on my screen, which is different from what I get when it's
> rotating. Is there any different behavior between shade3d and persp3d?
> Please find below a simple example.
>
> Actually, my final aim is to save each of the rotating graphs for creating a
> "clip". I commented below what I would use later on, but is there any way to
> do it better?
Take a look at ?movie3d. I think it does what you want.
Duncan Murdoch
>
> Thanks a lot for your help,
> Carlo Giovanni Camarda
>
>
> library(rgl)
> # simplified version from help(rgl.viewpoint)
> shade3d(oh3d())
> coo <- 1:360
> for(i in 1:length(coo)) {
> rgl.viewpoint(coo[i])
> #filename <- paste("pic", formatC(i,digits=1,flag="00"), ".png", sep="")
> #rgl.snapshot(filename)
> }
>
> # simple persp3d adaption
> f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
> x <- seq(-10, 10, length= 30); y <- x
> z <- outer(x, y, f)
> persp3d(x, y, z)
> coo <- 1:360
> for(i in 1:length(coo)) {
> rgl.viewpoint(coo[i])
> #filename <- paste("pic", formatC(i,digits=1,flag="00"), ".png", sep="")
> #rgl.snapshot(filename)
> }
>
>
>
> -
> Dr. Carlo Giovanni Camarda
> Research Scientist
> Max Planck Institute for Demographic Research
> Laboratory of Statistical Demography
> Konrad-Zuse-Straße 1
> 18057 Rostock - Germany
> Phone: +49 (0)381 2081 172
> Fax: +49 (0)381 2081 472
> cama...@demogr.mpg.de
> http://www.demogr.mpg.de/en/staff/camarda/default.htm
> -
>
>
> --
> This mail has been sent through the MPI for Demographic ...{{dropped:10}}
>
>
>
>
>
> __
> 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.
--
This mail has been sent through the MPI for Demographic Research. Should you
receive a mail that is apparently from a MPI user without this text displayed,
then the address has most likely been faked. If you are uncertain about the
validity of this message, please check the mail header or ask your system
administrator for assistance.
__
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] persp3d and rgl.viewpoint for rotating 3D plots
Dear R-users,
within the rgl-package, I would have a question about the usage of persp3d
in combination of rgl.viewpoint.
I am not able to figure out how to let a 3D plot rotating around likewise the
example in ?rgl.viewpoint. It seems that when I use persp3d(...) I see
something on my screen, which is different from what I get when it's rotating.
Is there any different behavior between shade3d and persp3d?
Please find below a simple example.
Actually, my final aim is to save each of the rotating graphs for creating a
"clip". I commented below what I would use later on, but is there any way to do
it better?
Thanks a lot for your help,
Carlo Giovanni Camarda
library(rgl)
# simplified version from help(rgl.viewpoint)
shade3d(oh3d())
coo <- 1:360
for(i in 1:length(coo)) {
rgl.viewpoint(coo[i])
#filename <- paste("pic", formatC(i,digits=1,flag="00"), ".png", sep="")
#rgl.snapshot(filename)
}
# simple persp3d adaption
f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
x <- seq(-10, 10, length= 30); y <- x
z <- outer(x, y, f)
persp3d(x, y, z)
coo <- 1:360
for(i in 1:length(coo)) {
rgl.viewpoint(coo[i])
#filename <- paste("pic", formatC(i,digits=1,flag="00"), ".png", sep="")
#rgl.snapshot(filename)
}
-
Dr. Carlo Giovanni Camarda
Research Scientist
Max Planck Institute for Demographic Research
Laboratory of Statistical Demography
Konrad-Zuse-Straße 1
18057 Rostock - Germany
Phone: +49 (0)381 2081 172
Fax: +49 (0)381 2081 472
cama...@demogr.mpg.de
http://www.demogr.mpg.de/en/staff/camarda/default.htm
-
--
This mail has been sent through the MPI for Demographic ...{{dropped:10}}
__
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] two plots side-by-side with persp3d(rgl)
Dear R-users,
I would like to place two 3D plots side-by-side in a rgl-setting. It
would nice to have something like "par(mfrow=c(1,2))" for basic plots, or an
array framework for wireframe(lattice) (see example below). I only managed to
overlap two persp3d plots.
My final idea would be to animate both surfaces using play3d(rgl).
Thanks in advance for any help.
Best,
Carlo Giovanni Camarda
# simple datasets
x <- seq(-10, 10, length= 30)
y <- x
f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
z1 <- outer(x, y, f)
z1[is.na(z1)] <- 1
z2 <- matrix(outer(1:30, 1:30)/200, 30, 30)
# simple persp
par(mfrow=c(1,2))
persp(x, y, z1)
persp(x, y, z2)
# wireframe
library(lattice)
grid. <- expand.grid(list(x=x,y=y,name=c("Z1", "Z2")))
grid.$Z <- c(c(z1), c(z2))
wireframe(Z ~ x * y | name, data = grid.)
# single plot
library(rgl)
persp3d(x, y, z1)
persp3d(x, y, z2)
# overlapping plots
persp3d(x,y,z1, color="#FF", front="lines")
persp3d(x,y,z2, col="red", add=TRUE)
play3d(spin3d(axis=c(0,0,1), rpm=10), duration=5)
-
Dr. Carlo Giovanni Camarda
Research Scientist
Max Planck Institute for Demographic Research
Laboratory of Statistical Demography
Konrad-Zuse-Straße 1
18057 Rostock - Germany
Phone: +49 (0)381 2081 172
Fax: +49 (0)381 2081 472
cama...@demogr.mpg.de
http://www.demogr.mpg.de/en/staff/camarda/default.htm
-
--
This mail has been sent through the MPI for Demographic ...{{dropped:10}}
__
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] matrix multiplication, tensor product, block-diagonal and fast computation
Dear Chuck,
thanks a lot for your reply.
Nevertheless, and though not so elegant, the for-loop solution seems faster
(see below).
Best,
Carlo Giovanni Camarda
m <- 50
n <- 40
nl <- 1000
A <- array(1:(m*m*n), dim=c(m,m,n))
M <- matrix(1:(m*n), nrow=m, ncol=n)
# for-loop solution
system.time(
for(j in 1:nl){
M1 <- matrix(nrow=m, ncol=n)
for(i in 1:n){
M1[,i] <- A[,,i] %*% M[,i]
}
})
# user system elapsed
# 5.650.035.72
# rowSums solution
system.time(
for(j in 1:nl){
M5 <- rowSums( aperm(A * as.vector( M[
rep(1:ncol(A),each=nrow(A)),]),c(1,3,2)), dims = 2 )
})
# user system elapsed
# 12.001.08 13.13
-
Dr. Carlo Giovanni Camarda
Research Scientist
Max Planck Institute for Demographic Research
Laboratory of Statistical Demography
Konrad-Zuse-Straße 1
18057 Rostock - Germany
Phone: +49 (0)381 2081 172
Fax: +49 (0)381 2081 472
cama...@demogr.mpg.de
http://www.demogr.mpg.de/en/staff/camarda/default.htm
-
-Original Message-
From: Charles C. Berry [mailto:cbe...@tajo.ucsd.edu]
Sent: Thursday, March 12, 2009 4:25 PM
To: Camarda, Carlo Giovanni
Cc: r-h...@stat.math.ethz.ch
Subject: Re: [R] matrix multiplication, tensor product, block-diagonal and fast
computation
On Wed, 11 Mar 2009, Camarda, Carlo Giovanni wrote:
> Dear R-users,
>
> I am searching to the "best" way to compute a series of n matrix
> multiplications between each matrix (mXm) in an array (mXmXn), and each
> column of a matrix (mXn).
>
> Please find below an example with four possible solutions.
> The first is a simple for-loop which one might avoid; the second
> solution employs the tensor product but then manually selects the right
> outcomes. The third solution uses too many (and too time-consuming)
> functions in order to profit of mapply. The last solution creates a
> block-diagonal from the array and multiply such matrix by the vectorized
> version of the first matrix. Here, often, the block-diagonal matrix may
> be too large and a specific list is needed (at least AFAIK).
>
> Does anyone have a further and possibly more effective way of computing
> such operation?
This is fairly quick:
rowSums( aperm(A * as.vector( M[ rep(1:ncol(A),each=nrow(A)),]),c(1,3,2)), dims
= 2 )
But my advice is to code this in C along the lines of your first solution
(using the BLAS routines to carry it out in the inner products). Your code
will be easier to read and debug and will probably be faster and easier on
memory, too.
Years ago I wrote a lot of stuff like this in native S. I 'optimized' the
heck out of it using tricks like the one above as I was debugging the
code. I had to rewrite the bulk of it in C anyway. The S code was so hard
to read that I could not migrate it to C bit by bit. I had to start over
and the effort spent debugging it in S was lost.
As an alternative you might try the 'jit' package on your first solution.
HTH,
Chuck
>
> Thanks in advance for any suggestion,
> Carlo Giovanni Camarda
>
>
> library(tensor)
> library(Matrix)
> A <- array(seq(0,1,length=48), dim=c(4,4,3))
> M <- matrix(1:12, nrow=4, ncol=3)
>
> # first solution (avoid the for-loop)
> M1 <- matrix(nrow=4, ncol=3)
> for(i in 1:3){
>M1[,i] <- A[,,i] %*% M[,i]
> }
> # second solution (direct picking of the right cols)
> A1 <- tensor(A, M, 2, 1)
> M2 <- cbind(A1[,1,1],A1[,2,2],A1[,3,3])
> # third solution (avoid as.data.frame and as.matrix)
> Adf0 <- apply(A, 3, as.data.frame)
> Adf1 <- lapply(X=Adf0, FUN=as.matrix, nrow=4, ncol=4)
> M3 <- mapply(FUN="%*%", Adf1, as.data.frame(M))
> # fourth solution (often too large block-diagonal matrix)
> Alist <- NULL
> for(i in 1:3){ # better way to create such list for bdiag
>Alist[[i]] <- A[,,i]
> }
> Abd <- bdiag(Alist)
> M4 <- matrix(as.matrix(Abd %*% c(M)), nrow=4, ncol=3)
>
> --
> This mail has been sent through the MPI for Demographic ...{{dropped:10}}
>
> __
> 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.
>
Charles C. Berry(858) 534-2098
Dept of Family/Preventive Medicine
E mailto:cbe...@tajo.ucsd.edu UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901
--
This mail has been sent through the MPI for Demographic Research. Should you
receive a mail that is apparently from a MPI user without this text
[R] matrix multiplication, tensor product, block-diagonal and fast computation
Dear R-users,
I am searching to the "best" way to compute a series of n matrix
multiplications between each matrix (mXm) in an array (mXmXn), and each
column of a matrix (mXn).
Please find below an example with four possible solutions.
The first is a simple for-loop which one might avoid; the second
solution employs the tensor product but then manually selects the right
outcomes. The third solution uses too many (and too time-consuming)
functions in order to profit of mapply. The last solution creates a
block-diagonal from the array and multiply such matrix by the vectorized
version of the first matrix. Here, often, the block-diagonal matrix may
be too large and a specific list is needed (at least AFAIK).
Does anyone have a further and possibly more effective way of computing
such operation?
Thanks in advance for any suggestion,
Carlo Giovanni Camarda
library(tensor)
library(Matrix)
A <- array(seq(0,1,length=48), dim=c(4,4,3))
M <- matrix(1:12, nrow=4, ncol=3)
# first solution (avoid the for-loop)
M1 <- matrix(nrow=4, ncol=3)
for(i in 1:3){
M1[,i] <- A[,,i] %*% M[,i]
}
# second solution (direct picking of the right cols)
A1 <- tensor(A, M, 2, 1)
M2 <- cbind(A1[,1,1],A1[,2,2],A1[,3,3])
# third solution (avoid as.data.frame and as.matrix)
Adf0 <- apply(A, 3, as.data.frame)
Adf1 <- lapply(X=Adf0, FUN=as.matrix, nrow=4, ncol=4)
M3 <- mapply(FUN="%*%", Adf1, as.data.frame(M))
# fourth solution (often too large block-diagonal matrix)
Alist <- NULL
for(i in 1:3){ # better way to create such list for bdiag
Alist[[i]] <- A[,,i]
}
Abd <- bdiag(Alist)
M4 <- matrix(as.matrix(Abd %*% c(M)), nrow=4, ncol=3)
--
This mail has been sent through the MPI for Demographic ...{{dropped:10}}
__
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] trade-off between speed and storage in matrix multiplications
Dear R-users,
I coded two equivalent ways to perform (in a simplified version)
some matrix multiplications I would like to use in a more general
framework.
In the first case I used Kronecker product and vectorization of a
certain matrix. This approach takes less time, but, as you may guess, I
run out of memory when dimensions are large.
In the second approach, I profited of sparseness and structure of the
matrices and use outer-functions for performing operations. Here it
takes more time, but I have not problem of memory.
Is there any way to enhance the second approach for speeding-up the
whole process? Or, in computing, this is a well-known trade-off between
speed and storage-spaces which I'm not aware (sorry for that).
Please have a look to the example below.
Needless to say that I'd appreciate any suggestion.
Best,
Carlo Giovanni
# dimensions
m <- 10
n <- 15
# A-matrix
rnA <- runif(m*m)
A <- matrix(rnA, m, m)
# vector
v <- runif(n)
# B-matrix
rnB <- runif(m*n)
B <- matrix(rnB, m, n)
# first solution: vectorize B + kronecker product => faster but storage
issues
system.time(
for(i in 1:100){
b <- c(B)
vKron.A <- kronecker(diag(v), A)
SOL1 <- vKron.A %*% b
})
# second solution: orig. dims + apply + mapply => slower, but w/o
storage issues
system.time(
for(i in 1:100){
Av <- outer(A, v, FUN="*")
Av.df1 <- apply(Av, 3, as.data.frame)
Av.df2 <- lapply(X=Av.df1, FUN=as.matrix, nrow=m, ncol=m)
SOL2 <- mapply(FUN="%*%", Av.df2, as.data.frame(B)) # as a matrix
})
--
This mail has been sent through the MPI for Demographic ...{{dropped:10}}
__
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] multiplication between a matrix and a block-diagonal matrix (without vectorization)
Dear R-users,
I would have a question regarding the multiplication between a
matrix (A, mXn) and a block-diagonal matrix, B, (m*n)X(m*n).
The easy solution would be to, first, vectorize A and use simple matrix
multiplication, but my aim is to avoid such vectorization. Is it
possible?
Please find below an explanatory example.
Thanks in advance for any help you can provide.
Best,
Carlo Giovanni
library(Matrix)
# matrix A
A <- matrix(1:12, 3, 4)
# elements of the block-diagonal matrix
M1 <- matrix(1:9,3,3)
M2 <- M1+1
M3 <- M2+1
M4 <- M3+1
# block-diagonal matrix
B <- bdiag(M1, M2, M3, M4)
# simple solution
a <- c(A)
B %*% a
--
This mail has been sent through the MPI for Demographic ...{{dropped:10}}
__
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] match.call in a function
Dear UseRs,
please find below a simple example which exemplifies my question.
In words, I'd like to write a function which prints formula and outcomes only
when results are not assigned or recalling the already assigned object, as e.g.
glm.
I'm sure a solution is rather simple and straightforward and it has something
to do with match.call (maybe), but I'm not able to see it.
Thanks in advance for any help
Carlo G. Camarda
# my simple function
myFUN <- function(a, b, method){
res <- ifelse(method==1, a+b, a-b)
object <- list(res=res, call=match.call())
return(object)
}
# assigning to A
A <- myFUN(2,4,1)
# recalling A or don't assign ...
A
myFUN(2,4,1)
# ... what I'd like (somehow) to have when
# I don't assign myFUN to another obj
# or when I recall the obj itself
My simple analysis
Call: myFUN(a = 2, b = 4, method = 1)
a: 2
b: 4
Outcome: 6
-
Carlo Giovanni Camarda
Research Scientist
Konrad-Zuse-Straße 1
18057 Rostock - Germany
Phone: +49 381 2081 172
Fax: +49 381 2081 472
[EMAIL PROTECTED]
http://www.demogr.mpg.de/en/staff/camarda/default.htm
-
--
This mail has been sent through the MPI for Demographic ...{{dropped:10}}
__
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] Summing up diagonals w/o for-loop
Dear R-users,
is there any way to sum up the elements of the "diagonals" of a matrix
without using a for-loop? While there is a simple way over rows and
columns, I don't see a straightforward multiplication for the diagonals,
am I too demanding? Or, more likely, I'm lack some algebra trick? Is
there any R-function that can deal with this problem w/o loop?
Actually I would need to sum up just the upper diagonals.
Here a simple, but general, example for presenting the problem.
Thanks in advance for any help,
Carlo Giovanni Camarda
m <- 7
n <- 5
mat <- matrix(1:35, n, m)
ones.r <- rep(1,n)
ones.c <- rep(1,m)
# sum over the rows
sum.r <- mat%*%ones.c
# sum over the cols
sum.c <- t(mat)%*%ones.r
# sum over the diags
sum.d <- numeric(m+n-1)
sum.d[1] <- mat[n,1]
sum.d[m+n-1] <- mat[1,m]
for(i in 2:n){
sum.d[i] <- sum(diag(mat[(n+1-i):n,1:i]))
}
for(i in 2:(m-1)){
sum.d[i+n-1] <- sum(diag(mat[,i:m]))
}
--
This mail has been sent through the MPI for Demographic ...{{dropped:10}}
__
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.
