Difference between revisions of "ExamineConvest.R"
From Organic Design wiki
(Much more efficient way of obtaining iterations) |
(tol tests) |
||
| (One intermediate revision by the same user not shown) | |||
| Line 1: | Line 1: | ||
#{{R}} | #{{R}} | ||
library(limma) | library(limma) | ||
| − | + | ||
# Load local copy of convest | # Load local copy of convest | ||
source("convest.R") | source("convest.R") | ||
| − | + | ||
# ------------- Generating a mixture of normally distributed data ------------- # | # ------------- 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) { | ||
| Line 15: | 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 27: | 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) | ||
| + | |||
| − | convest(pvalues, niter= | + | # 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) | conv <- read.table("conv.txt", header=TRUE, as.is=TRUE) | ||
# Interestingly the convergence rate varies if m is small | # Interestingly the convergence rate varies if m is small | ||
| − | |||
plot(conv$pi0, type="l") | plot(conv$pi0, type="l") | ||
plot(diff(conv$pi0), type="l") | plot(diff(conv$pi0), type="l") | ||
| − | # | + | # Weighted average of convest iterations after removing burnin period |
# Raw | # Raw | ||
mean(conv$pi0[-(1:10)]) | mean(conv$pi0[-(1:10)]) | ||
| Line 48: | Line 51: | ||
weighted.mean(conv$pi0[-(1:10)], w=seq(0,1, length=length(conv$pi0[-(1:10)]))) | weighted.mean(conv$pi0[-(1:10)], w=seq(0,1, length=length(conv$pi0[-(1:10)]))) | ||
| − | convest(pvalues, niter= | + | # 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") | ||
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")
