用Emacs 上批踢踢

December 31, 2009

其實很簡單。

在emacs下執行

M-x ansi-term

然後在新開的ansi-term的shell下執行

luit -encoding big5 telnet bbs.ptt.cc

完畢。

示範圖如下:

Advertisements

Consider $\mathbf x\sim \mathscr N(\mathbf\mu, \mathbf\Sigma)$ and $\mathbf x = (x_1, \mathbf x_2)$
we have the following R implementation for obtaining the conditional quantile of this conditional distribution:


qCondMvNorm <- function(P, x, cpar=1, meanVec, covMat) {
  x1 <- x[cpar]
  x2 <- x[-cpar]
  mu1 <- meanVec[cpar]
  mu2 <- meanVec[-cpar]
  sigma11 <- covMat[cpar, cpar]
  sigma22 <- covMat[-cpar, -cpar]
  sigma12 <- covMat[cpar, -cpar]
  sigma21 <- covMat[-cpar, cpar]
  sd11.2 <- sqrt(sigma11 - sigma12%*%solve(sigma22)%*%sigma21)
  mu1.c2 <- mu1 + sigma12%*%solve(sigma22)%*%(x2 - mu2)
  qmx <- qnorm(P, mean=mu1.c2, sd=sd11.2)
  qmx
}


The following Python is used for finding conditional quantile:

import scipy.stats.distributions as ssd

x = 0.34
a = 1.2
b = 1.8
X = ssd.norm.rvs(loc=0.5, scale=a, size=100)
Y = ssd.norm.rvs(loc=-1.5, scale=b, size=100)
w = ssd.norm.pdf(x, loc=X, scale=a)/sum(ssd.norm.pdf(x, loc=X, scale=a))

def qMixCondNorm(P, x, X, Y, a, b, tol):
    w = ssd.norm.pdf(x, loc=X, scale=a)/sum(ssd.norm.pdf(x, loc=X, scale=a))
    yk = 0.0
    P_yk = sum(w*ssd.norm.cdf(yk, loc=Y, scale=b))
    yk_ud = 2.0
    while (max(abs(P - P_yk), abs(yk -yk_ud)) > tol):
        yk = yk_ud
        noms = 2*(sum(w*ssd.norm.cdf(yk, loc=Y, scale=b))-P)*(sum(w*ssd.norm.pdf(yk, loc=Y, scale=b)))/b
        denoms = 2*sum(w*ssd.norm.pdf(yk, loc=Y, scale=b))/b - 2*(sum(w*ssd.norm.cdf(yk, loc=Y, scale=b)) - P)*(sum(w*(yk - Y)*ssd.norm.pdf(yk, loc=Y, scale=b)))/(b**2)
        yk_ud = yk - noms/denoms
        P_yk = sum(w*ssd.norm.cdf(yk_ud, loc=Y, scale=b))
    return yk


y = qMixCondNorm(0.95, x, X, Y, a, b, 2*1e-14)

sum(w*ssd.norm.cdf(y, loc=Y, scale=b))

Update: 2009 Dec 25

Using
\[ x_{t+1} = x_t – \frac{f'(x_t)}{f”(x_t)}\]

Mixture normal’s quantile

December 21, 2009

How to get the following conditional distribution’s quantile?
\[\hat{f}(y|x) = \sum_{i=1}^N\frac{\left(\frac{x – X_i}{a}\right)}{b\sum_{i=1}^N\left(\frac{x – X_i}{a}\right)}\left(\frac{y – Y_i}{b}\right)\]

The R code (simple N-R method):

qMixConF <- function(X, Y, P, x, a, b) {
   w <- pnorm(x, X, a)/sum(pnorm(x, X, a))
   yk <- 0
   yk.ud <- 100
   for ( yk.ud != yk) {
       yk.ud <- yk- b*(sum(w*pnorm(y, Y, b))-P)/(sum(w*dnorm(y, Y, b)))
       yk <- yk.ud
   }
   yk
}

戰文範例

December 13, 2009

這篇部落文真是戰文的絕佳範例, 要戰文的話最少要達到這樣的水準吧

在用Rsurvival套件中提供的coxph作cox proportional hazard model 引進時間相關解釋變數時,
\[ h(t_i) = h_0(t_i)\exp[\mathbf X_i\boldsymbol\beta_i]\]
我們會發現處理這樣的變數是有困難的, 不能直接套用。
下面的R函式可以幫我們解決這個問題:

expand.breakpoints <- function(dataset, index = "patnum", status = "status",
                               tevent = "time", Zvar = F, breakpoints = NULL) {
  ## Expand  so that each individual has a row for each
  ## unique failure time that is <= his/her/its own failure time.
  ##
  ## Original version written by Kathy Russell, for
  ## Kenneth Hess, Department of Biomathematics,
  ## University of Texas M. D. Anderson Cancer Centre.
  ##
  ## Substantially modified by JHM, Nov 23, 1998.
  ##
  ## ERROR checking
  
  onceonly <- paste("must be a character string matching exactly",
                    "one name in the first parameter, dataset")
  ## Require dataset to be of type data.frame
  if((missing(dataset) || !is.data.frame(dataset)))
    stop("\n Parameter, dataset, must be a data frame")
  varset <- names(dataset)
  covset <- varset
  lvset <- 1:length(varset) # Require dataset to have unique names
  if(any(duplicated(varset))) {
    stop(paste("\n Parameter, dataset, must have uniquely defined",
               "column names"))
  }
  ## Require index to match exactly 1 dataset name
  if(length((indexloc <- lvset[varset == index])) != 1)
    stop(paste("\n Parameter, index,", onceonly))
  covset <- covset[covset != index]
  ## Require status to match exactly 1 dataset name
  if(length((statusloc <- lvset[varset == status])) != 1)
    stop(paste("\n Parameter, status,", onceonly))
  covset <- covset[covset != status]
  ## Require tevent to match exactly 1 dataset name
  if(length((teventloc <- lvset[varset == tevent])) != 1)
    stop(paste("\n Parameter, tevent,", onceonly))
  covset <- covset[covset != tevent] # -----------------------------
  ## Form vector of breakpoints, if necessary
  if(is.null(breakpoints)) {
    times <- dataset[, tevent]
    breakpoints = 1
  if((is.null(breakpoints)) || (!(is.numeric(breakpoints))) ||
     (!(is.vector(breakpoints))) || (length(breakpoints) < 1))
    stop(paste("\n Parameter, breakpoints, must be a numeric vector",
               "with at least 1 element"))


  ##*****************
  ## Begin
  ##*****************
  n <- nrow(dataset)
  temp.stop  0)
    temp.start <- c(0, breakpoints)
  else temp.start <- c(-1, breakpoints)
  n.break <- length(breakpoints) + 1
  t.event <- dataset[, tevent]
  
  ## ---------------------------
  ## Begin calculate n.epochs

  n.epochs <- sapply(1:n, function(m, breaks, times)
                     sum(breaks < times[m]) + 1, breaks = breakpoints, times = t.event)
  ## End n.epochs
  id <- rep(1:n, n.epochs) ## Index into supplied dataset
  last.epoch <- cumsum(n.epochs)
  epoch <- unlist(sapply(n.epochs, function(x)
                         1:x)) ## Index into vectors of interval start & stop points
  if(Zvar) {
    Zmat <- diag(rep(1, n.break))[epoch, ]
    Z.name <- paste("Z", seq(1, n.break), sep = "")
  }
  Tstop <- temp.stop[epoch]
  Tstop[last.epoch] <- dataset[, tevent]
  status2 <- rep(0, length(epoch))
  status2[last.epoch] <- dataset[, status]
  new.dataset <- data.frame(dataset[id, index, drop = F], temp.start[epoch],
                            Tstop, status2, epoch, dataset[id, covset, drop = F])
  if(Zvar)
    new.dataset <- data.frame(new.dataset, Zmat)
  nam <- c(index, "Tstart", "Tstop", status, "epoch", covset)
  if(Zvar)
    nam <- c(nam, Z.name)
  dimnames(new.dataset) <- list(1:length(epoch), nam)
  return(new.dataset)
}

SAS則是非常簡單, 可參閱phreg的文件。

If the partial likelihood
\[ \prod_{i=1}^D\frac{\exp[\sum_{k=1}^p\beta_kZ_{(i)k}]}{\sum_{i\in R(t_i)}\exp[\sum_{k=1}^p\beta_k Z_{jk}]}\]
with Fisher score vector
\[ U_b(\boldsymbol\beta) = \sum_{i=1}^DZ_{(i)b} – \sum_{i=1}^D\frac{\sum_{j\in R(t_i)}Z_{jb\exp[\sum_{k=1}^p\beta_kZ_{jk}]}}{\sum_{j\in R(t_i)}\exp[\sum_{k=1}^p\beta_kZ_{jk}]}\]