Difference between revisions of "GGB-correlation-snippet.R"
m (fix k param) |
(Fixing/moving plots) |
||
| Line 12: | Line 12: | ||
params <- c(2.74886, 1.36546, 4.12844) # "IPTG-a parameters" | params <- c(2.74886, 1.36546, 4.12844) # "IPTG-a parameters" | ||
| − | k <- | + | k <- 1000 |
a.shape <- params[1] | a.shape <- params[1] | ||
| Line 19: | Line 19: | ||
kDE <- rgamma(2*m1, shape=k, rate=k) # p49 2001 NEWTON paper | kDE <- rgamma(2*m1, shape=k, rate=k) # p49 2001 NEWTON paper | ||
| − | |||
| − | |||
DEscales <- rgamma(2*m1, shape=a0.shape, rate=scale) | DEscales <- rgamma(2*m1, shape=a0.shape, rate=scale) | ||
DEscales <- DEscales * kDE | DEscales <- DEscales * kDE | ||
kEE <- rgamma(m0, shape=k, rate=k) | kEE <- rgamma(m0, shape=k, rate=k) | ||
| − | |||
| − | |||
EEscales <- rgamma(m0 , shape=a0.shape, rate=scale) | EEscales <- rgamma(m0 , shape=a0.shape, rate=scale) | ||
EEscales <- EEscales / kEE | EEscales <- EEscales / kEE | ||
| Line 37: | Line 33: | ||
# Graphical check | # Graphical check | ||
| + | hist(kDE, prob=TRUE) | ||
| + | lines(density(kDE)) | ||
| + | |||
| + | hist(kEE, prob=TRUE) | ||
| + | lines(density(kEE)) | ||
| + | |||
| + | Xbar <- rowMeans(log2(X)) | ||
| + | hist(Xbar, breaks=30, prob=TRUE) | ||
| + | lines(density(Xbar)) | ||
pairs(log2(X)) | pairs(log2(X)) | ||
| − | |||
Revision as of 23:58, 28 August 2006
- p (genes) by n (slides) matrix
- m=p keeping fdr notation
m <- p <- 1000 nreps <- 5
n <- nreps * 2 pi0 <- 0.95
m1 <- round(m * (1-pi0)) m0 <- p - m1
params <- c(2.74886, 1.36546, 4.12844) # "IPTG-a parameters" k <- 1000
a.shape <- params[1] a0.shape <- params[2] scale <- params[3]
kDE <- rgamma(2*m1, shape=k, rate=k) # p49 2001 NEWTON paper DEscales <- rgamma(2*m1, shape=a0.shape, rate=scale) DEscales <- DEscales * kDE
kEE <- rgamma(m0, shape=k, rate=k) EEscales <- rgamma(m0 , shape=a0.shape, rate=scale) EEscales <- EEscales / kEE
scales <- c(rep(DEscales, each=nreps), rep(EEscales, each=2*nreps)) X <- rgamma(n* p, a.shape, rate=scales) dim(X) <- c(n,p)
X <- t(X)
- Graphical check
hist(kDE, prob=TRUE) lines(density(kDE))
hist(kEE, prob=TRUE) lines(density(kEE))
Xbar <- rowMeans(log2(X)) hist(Xbar, breaks=30, prob=TRUE) lines(density(Xbar)) pairs(log2(X))
