Difference between revisions of "ExamineConvest.R"
From Organic Design wiki
(Much more efficient way of obtaining iterations) |
(version checking using limma::convest(pvalues, niter=50)) |
||
| Line 1: | Line 1: | ||
#{{R}} | #{{R}} | ||
| + | Sys.putenv("DISPLAY"=":0") | ||
| + | |||
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 17: | ||
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)) | ||
| − | + | 1 m0 <- m - m1 | |
means <- rep(0,m) | means <- rep(0,m) | ||
#c(rep(mu1,m1), rep(mu0,m0)) | #c(rep(mu1,m1), rep(mu0,m0)) | ||
| Line 27: | Line 29: | ||
return(pvalues) | return(pvalues) | ||
} | } | ||
| − | + | ||
# ---------------------------------------------------------------------------- # | # ---------------------------------------------------------------------------- # | ||
# Multtest does not allow a max of 1 (NA) going into smooth.spline | # Multtest does not allow a max of 1 (NA) going into smooth.spline | ||
| − | + | ||
#convest arguement doreport=TRUE provides iteration information | #convest arguement doreport=TRUE provides iteration information | ||
| − | pvalues <- simNorm2(m= | + | pvalues <- simNorm2(m=500) |
| − | + | ||
| − | convest(pvalues, niter= | + | convest(pvalues, niter=50, 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)]))) | ||
| + | # Checking versions | ||
convest(pvalues, niter=50) | convest(pvalues, niter=50) | ||
| + | limma::convest(pvalues, niter=50) | ||
Revision as of 11:55, 11 June 2007
Code snipits and programs written in R, S or S-PLUS Sys.putenv("DISPLAY"=":0")
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))
1 m0 <- m - m1
means <- rep(0,m) #c(rep(mu1,m1), rep(mu0,m0)) 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)
}
- ---------------------------------------------------------------------------- #
- Multtest does not allow a max of 1 (NA) going into smooth.spline
- convest arguement doreport=TRUE provides iteration information
pvalues <- simNorm2(m=500)
convest(pvalues, niter=50, 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)])))
- Checking versions
convest(pvalues, niter=50) limma::convest(pvalues, niter=50)
