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)) ...) |
(Much more efficient way of obtaining iterations) |
||
| 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 23: | Line 27: | ||
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 | ||
| Line 29: | Line 34: | ||
pvalues <- simNorm2(m=50) | pvalues <- simNorm2(m=50) | ||
| − | + | convest(pvalues, niter=500, doreport=TRUE, file="conv.txt") | |
| − | + | 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(pi0 | + | plot(conv$pi0, type="l") |
| − | + | plot(diff(conv$pi0), type="l") | |
| − | plot(diff(pi0), type="l") | ||
# weighted average of convest iterations after removing burnin period | # 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)]))) |
convest(pvalues, niter=50) | convest(pvalues, niter=50) | ||
| − | |||
Revision as of 04:18, 11 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) #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=50)
convest(pvalues, niter=500, 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)])))
convest(pvalues, niter=50)
