[prev in list] [next in list] [prev in thread] [next in thread] 

List:       r-help
Subject:    [R] running Apply takes long time
From:       Hamid Ghorbani via R-help <r-help () r-project ! org>
Date:       2019-06-30 9:28:01
Message-ID: 293554377.561656.1561886881393 () mail ! yahoo ! com
[Download RAW message or body]

Dear users,
Below I   have a   matrix, called *mysim.obs*     (548 row and (1+nsim) columns, \
which its first column is my observation and the next   columns comprises simulation \
from the   fitted   model to first column). I want to evaluate *simpsonlognormpval* \
function on each column of   *mysim.obs*. For this I have used apply function. \
Unfortunately, running apply   takes long time (i have several models, log-normal \
model in the following is just for   explanation). Many thanks in advance.
Yours,
Hamid
------------

simpsonlognormpval <- function(xx){
       # numerical integral using Simpson's rule
       # assume a < b and n is an even positive integer
                   n<-10000
a<-0;b<- 25*max(xx) #because log-normal distribution has heavy tail
meanlog = -0.216
sdlog = 1.4245521
       h <- (b-a)/n
       x <- seq(a, b, by=h)
                               y <- (plnorm(x, meanlog = meanlog0 , sdlog =sdlog0 \
)-ecdf(xx)(x))^2  if (n == 2) {
               s <- (y[1] + 4*y[2] +y[3])
       } else {
               s <- y[1] + y[n+1] + 2*sum(y[seq(2,n,by=2)]) + 4 *sum(y[seq(3,n-1, \
by=2)])  }
       s <- s*h/3
       return(s)
}
  
> meanlog0 = -0.216
> sdlog0 = 1.4245521
> nsim=100000
> my.obs<-rexp( 548,0.5*lambda ) #my.obs is acctually an observed sample, here I just \
> replaced it mysim.obs<-cbind(my.obs ,matrix(rlnorm(548*nsim, meanlog = meanlog0, \
> sdlog =sdlog0),548,nsim))

> fsimpsonlognormpval <-apply( mysim.obs, 2,simpsonlognormpval )
> fsimpsonlognormpval[1]
> lognormpval<-mean(fsimpsonlognormpval[2:(nsim+1)]>fsimpsonlognormpval[1])
> lognormpval  
	[[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.


[prev in list] [next in list] [prev in thread] [next in thread]