[ESS] Emacs taking over CPU (Frank Harrell)

Rodney Sparapani rsparapa at mcw.edu
Tue Sep 17 18:10:31 CEST 2013


On 09/16/2013 05:19 PM, Vitalie Spinu wrote:
> How large your output should be to take so long?
>
> The emacs display is notoriously slow. I have fixed this back in
> March. We are now caching output in temporary buffers before flushing
> into comint buffer every 0.6 seconds. So, the output speed is pretty
> much as fast as emacs can get it. To be more concrete I would need a
> concrete example.
>
>    Vitalie

Hi Vitalie:

I can boil it down to this.  Try running long-test.R in an *R* buffer 
vs. at the command line (attached).

-- 
Rodney Sparapani, PhD
Manager of Statistical & Computational Operations
Center for Patient Care and Outcomes Research (PCOR)
Medical College of Wisconsin (MCW), Milwaukee, USA
http://www.linkedin.com/in/rodneysparapani
-------------- next part --------------
# R --no-save < long-test.R > long-test.Rt &

source("sim.bbjem.R")

    N <- 100
    p <- 0
    p1 <- 1
    q <- 1
    r <- p1+q
    s <- p1+r

delta <- 6
    gamma <- 0
    eta  <- gamma
    beta <- 0.2
    mu   <- 0 # 0.25
    rho  <- 0.6
    alpha <- 1

set.seed(66)
 
    info <- list(beta=list(init=beta, prior=list(mean=0., prec=0.001)),
             rho=list(init=rho),
             mu=list(init=rep(mu, 2)),
             theta=list(init=c(rep(gamma, p1), rep(delta, q), rep(eta, p1)),
               prior=list(mean=rep(0., s), prec=diag(0.001, s))),
             dpm=list(m=as.integer(3), alpha=alpha,
                 C=as.integer(0*(1:N)), states=as.integer(N),
               prior=list(mu0=c(0., 0.), T0=diag(0.1, 2), S0=diag(10, 2))))
    
for(i in 1:225) {
    print(sim.bbjem(N=N, p=p1, q=q,
        	gamma=info$theta$init[1:p1],
        	delta=info$theta$init[(p1+1):r],
        	eta=info$theta$init[(r+1):s],
        	beta=info$beta$init,
        	sigma=info$rho$init,
        	mu=info$mu$init))
}
-------------- next part --------------
# This function uses latent t to generate latent y then makes both binary
require(MASS)

sim.bbjem <- function(N=100,   # sample size
                    p=1,       # number of x covariates
                    q=10,      # number of instruments in z
                    eta=0,     # coefficient of x in outcome equation
                    beta=1,    # coefficient of T in outcome equation
                    delta=1,   # coefficient of instrumental variables
                    gamma=0,   # coefficient of x in treatment choice equation
                    rho=0.6,   # correlation
                    sigma=1,   # standard deviation
                    mu=1,      # intercept in the two equations
                    zbin=FALSE,# Z is continuous
                    Z=NA,      # by default, Z is randomly generated  
                    dist="normal"
                                # by default, the underlying true distribution
                                # is normal
                                # for Logistic, use "logistic"
                                # for a mixture of normal, use "mixture"       
                   )
{
  
  if (p == 0) {
    p <- 1
    x.mat <- matrix(0, N, p)
  }
  else x.mat <- matrix(runif(N*p, -0.5, 0.5), N, p)

  if ( length(gamma) == 1 ) gamma <- rep(gamma, p)
  if ( length(eta  ) == 1 ) eta   <- rep(eta  , p)
  if ( length(delta) == 1 ) delta <- rep(delta, q)

  if ( length(mu) == 1 ) mu <- rep(mu, 2)

  if(length(rho)==1){
    if(abs(rho)>=1){rho <- 0}
    Sigma <- diag(sigma^2,2,2)
    Sigma[1,2] <- rho*sigma^2
    Sigma[2,1] <- rho*sigma^2
  }

  attr(x.mat, "dimnames") <- list(NULL, c(paste("x", 1:p, sep="")))

  if (!any(is.na(Z))) z.mat <- matrix(Z, N, q)
  else if (zbin) z.mat <- matrix(rbinom(N*q, 1, 0.5)-0.5, N, q)
  else z.mat <- matrix(runif(N*q, -0.5, 0.5), N, q)
  
  attr(z.mat, "dimnames") <- list(NULL, c(paste("z", 1:q, sep="")))

  if(dist=="normal") Error <- mvrnorm(N,mu,Sigma)
  else if(dist=="mixture") Error <-1.5*((-1)^(runif(N)<0.5))+mvrnorm(N, mu, Sigma)
  else if(dist=="logistic") {
    Error <- mvrnorm(N, mu, Sigma)

    Error[ , 1] <- qlogis(pnorm(Error[ , 1], mean=mu[1], sd=sqrt(Sigma[1, 1])), 
                          location=mu[1], scale=sqrt(Sigma[1, 1]))

    Error[ , 2] <- qlogis(pnorm(Error[ , 2], mean=mu[2], sd=sqrt(Sigma[2, 2])), 
                          location=mu[2], scale=sqrt(Sigma[1, 1]))
  }
  else if(dist=="logistic.sd1") {
    Error <- mvrnorm(N, mu, Sigma)
## Logistic scale=1 ==> sd=pi/sqrt(3)
    Error[ , 1] <- qlogis(pnorm(Error[ , 1], mean=mu[1], sd=sqrt(Sigma[1, 1])), 
                          location=mu[1], scale=sqrt(3*Sigma[1, 1])/pi)

    Error[ , 2] <- qlogis(pnorm(Error[ , 2], mean=mu[2], sd=sqrt(Sigma[2, 2])), 
                          location=mu[2], scale=sqrt(3*Sigma[2, 2])/pi)
  }
  else if(dist=="logistic.ms4") {
    stopifnot(sigma==1)
    
    p.r <- c(0.106498992656952,0.458361227014536,0.374189066914829,
                  0.060950713413683)
    s.r <- 1/c(1.023135580500914,0.698774355946609,0.475127524640229,
                  0.321064655834542)

    #Error <- matrix(mu, nrow=N, ncol=2, byrow=TRUE)
    Error <- mvrnorm(N, c(0, 0), Sigma)

    for(i in 1:N) {
      #Error[i, ] <- Error[i, ]+mvrnorm(1, c(0, 0), Sigma)*
      Error[i, ] <- mu+Error[i, ]*s.r[c(c(1:4) %*% rmultinom(1, 1, p.r))]
    }
  }
  
  t.mat <- z.mat %*% delta + x.mat %*% gamma + Error[,1]
  
  tbin <- integer(N) #matrix(0,nrow=N,ncol=1)
  
  tbin[t.mat[,1]>0] <- 1
  
#  y <- beta*tbin +x.mat %*% eta +Error[,2]
  y <- beta*t.mat +x.mat %*% eta +Error[,2]
  
  ybin <- integer(N) #matrix(0,nrow=N,ncol=1)
  
  ybin[y[,1]>0] <- 1
  
data <- list(tbin=as.integer(tbin), X=x.mat, ybin=as.integer(ybin), Z=z.mat,
             dist=dist, tlat=c(t.mat), ylat=c(y), Error=Error)
  truetheta <- list(delta=delta,gamma=gamma,beta=beta,eta=eta,mu=mu,Sigma=Sigma)
  
  return(list(data=data,truetheta=truetheta))
}



More information about the ESS-help mailing list