Difference between revisions of "ExamineConvest.R"
From Organic Design wiki
(New page: library(limma) # ============= Generating a mixture of normally distributed data ============== simNorm <- function(m=1000, pi0=0.8, mu0=0, mu1=3, sigma=1) { m1 <- round(m * (1-pi0)) ...) |
(tol tests) |
||
| (2 intermediate revisions by the same user not shown) | |||
| Line 1: | Line 1: | ||
| + | #{{R}} | ||
library(limma) | library(limma) | ||
| − | + | ||
| − | # | + | # Load local copy of convest |
| + | source("convest.R") | ||
| + | |||
| + | # ------------- Generating a mixture of normally distributed data ------------- # | ||
simNorm <- function(m=1000, pi0=0.8, mu0=0, mu1=3, sigma=1) { | simNorm <- function(m=1000, pi0=0.8, mu0=0, mu1=3, sigma=1) { | ||
m1 <- round(m * (1-pi0)) | m1 <- round(m * (1-pi0)) | ||
| Line 11: | Line 15: | ||
return(pvalues) | return(pvalues) | ||
} | } | ||
| − | + | ||
simNorm2 <- function(m=1000, pi0=0.8, mu0=0, s1=1, s2=3) { | simNorm2 <- function(m=1000, pi0=0.8, mu0=0, s1=1, s2=3) { | ||
m1 <- round(m * (1-pi0)) | m1 <- round(m * (1-pi0)) | ||
m0 <- m - m1 | m0 <- m - m1 | ||
means <- rep(0,m) | means <- rep(0,m) | ||
| − | |||
sigma <- c(rep(s2,m1), rep(s1,m0)) | sigma <- c(rep(s2,m1), rep(s1,m0)) | ||
DEflag <- as.logical(sigma - min(sigma)) | DEflag <- as.logical(sigma - min(sigma)) | ||
| Line 23: | Line 26: | ||
return(pvalues) | return(pvalues) | ||
} | } | ||
| − | # | + | # ---------------------------------------------------------------------------- # |
| − | # | ||
#convest arguement doreport=TRUE provides iteration information | #convest arguement doreport=TRUE provides iteration information | ||
| − | pvalues <- | + | pvalues <- simNorm(m=500, pi0=0.25) |
| + | |||
| − | + | # 1) Illustrating that modified convest function is identical to limma version | |
| − | + | limma::convest(pvalues, niter=50, doreport=TRUE) | |
| − | + | convest(pvalues, niter=50, doreport=TRUE) | |
| − | |||
| − | |||
| − | + | # 2) Writing report information to file | |
| + | convest(pvalues, niter=100, doreport=TRUE, file="conv.txt") | ||
| + | conv <- read.table("conv.txt", header=TRUE, as.is=TRUE) | ||
| − | plot(diff(pi0), type="l") | + | # Interestingly the convergence rate varies if m is small |
| + | plot(conv$pi0, type="l") | ||
| + | plot(diff(conv$pi0), type="l") | ||
| − | # | + | # Weighted average of convest iterations after removing burnin period |
# Raw | # Raw | ||
| − | mean(pi0[-(1:10)]) | + | mean(conv$pi0[-(1:10)]) |
# Weighted mean | # Weighted mean | ||
| − | weighted.mean(pi0[-(1:10)], w=seq(0,1, length=length(pi0[-(1:10)]))) | + | weighted.mean(conv$pi0[-(1:10)], w=seq(0,1, length=length(conv$pi0[-(1:10)]))) |
| + | |||
| + | # 3) Tolerance adjustment | ||
| + | convest(pvalues, niter=300, doreport=TRUE, file="conv.txt", tol=1e-6) | ||
| + | conv <- read.table("conv.txt", header=TRUE, as.is=TRUE) | ||
| − | + | X11() | |
| − | pi0 | + | plot(conv$pi0, type="l") |
| + | plot(diff(conv$pi0), type="l") | ||
Latest revision as of 23:55, 24 June 2007
Code snipits and programs written in R, S or S-PLUS library(limma)
- Load local copy of convest
source("convest.R")
- ------------- Generating a mixture of normally distributed data ------------- #
simNorm <- function(m=1000, pi0=0.8, mu0=0, mu1=3, sigma=1) {
m1 <- round(m * (1-pi0)) m0 <- m - m1 means <- c(rep(mu1,m1), rep(mu0,m0)) DEflag <- as.logical(means) X <- rnorm(m, means, sigma) pvalues <- 2*(1-pnorm(abs(X),0,sigma)) return(pvalues)
}
simNorm2 <- function(m=1000, pi0=0.8, mu0=0, s1=1, s2=3) {
m1 <- round(m * (1-pi0)) m0 <- m - m1 means <- rep(0,m) sigma <- c(rep(s2,m1), rep(s1,m0)) DEflag <- as.logical(sigma - min(sigma)) X <- rnorm(m, means, sigma) pvalues <- 2*(1-pnorm(abs(X),0,s1)) return(pvalues)
}
- ---------------------------------------------------------------------------- #
- convest arguement doreport=TRUE provides iteration information
pvalues <- simNorm(m=500, pi0=0.25)
- 1) Illustrating that modified convest function is identical to limma version
limma::convest(pvalues, niter=50, doreport=TRUE) convest(pvalues, niter=50, doreport=TRUE)
- 2) Writing report information to file
convest(pvalues, niter=100, doreport=TRUE, file="conv.txt") conv <- read.table("conv.txt", header=TRUE, as.is=TRUE)
- Interestingly the convergence rate varies if m is small
plot(conv$pi0, type="l") plot(diff(conv$pi0), type="l")
- Weighted average of convest iterations after removing burnin period
- Raw
mean(conv$pi0[-(1:10)])
- Weighted mean
weighted.mean(conv$pi0[-(1:10)], w=seq(0,1, length=length(conv$pi0[-(1:10)])))
- 3) Tolerance adjustment
convest(pvalues, niter=300, doreport=TRUE, file="conv.txt", tol=1e-6) conv <- read.table("conv.txt", header=TRUE, as.is=TRUE)
X11() plot(conv$pi0, type="l") plot(diff(conv$pi0), type="l")
