Skip to content

Instantly share code, notes, and snippets.

@abikoushi
abikoushi / pois_anim.r
Created August 27, 2024 16:24
A Comparison between confidence interval vs. credible interval
library(animation)
pvfun <- Vectorize(function(r, y){
poisson.test(y, T = tau, r = r)$p.value
})
postprob <- function(r, x, tau){
pl <- pgamma(r, shape = x+0.5, rate = tau, lower.tail = FALSE)
pu <- pgamma(r, shape = x+0.5, rate = tau, lower.tail = TRUE)
2*pmin(pl,pu)
}
@abikoushi
abikoushi / wilcox_compositional.r
Created August 26, 2024 09:25
An simulation of wilcox.test for the compositional data
library(gtools)
library(ggplot2)
library(ggbrick)
set.seed(1)
wp <- numeric(10000)
ksp <- numeric(10000)
tp <- numeric(10000)
for(i in 1:10000){
y1 <- rdirichlet(50,rep(1,10))
y2 <- rdirichlet(50,rep(10,10))
library(KFAS)
set.seed(123456)
alpha <- beta <- numeric(100)
alpha[1:2] <- rnorm(2)
for(i in 3:100){
alpha[i] <- rnorm(1,2*alpha[i-1]-alpha[i-2],0.2)
}
beta <- rnorm(100,alpha,0.2)
y1<-rnbinom(100,size=10,mu=exp(alpha))
y2<-rnbinom(100,size=10,mu=exp(beta))
\documentclass[tikz,border=14pt, dvipsnames]{standalone}
\usepackage{tikz}
\usetikzlibrary{matrix,shapes,decorations.pathreplacing, backgrounds, positioning}
\usepackage{amssymb}
\begin{document}
\begin{tikzpicture}[
%Styles
Matrix/.style={
matrix of nodes,
@abikoushi
abikoushi / poistest.R
Created August 18, 2024 02:31
Check `poisson.test`
library(rootSolve)
pvfun0 <- function(r, x, tau){
pu <- ppois(x-1, r*tau, lower.tail = FALSE)
pl <- ppois(x, r*tau, lower.tail = TRUE)
2*pmin(pl, pu)
}
sol0 <-uniroot.all(f = function(p){pvfun0(p,5,10)-0.05},
interval = c(0.1,3))
sol0
res_p
@abikoushi
abikoushi / binomtest.R
Created August 17, 2024 15:59
Check `binom.test`
library(rootSolve)
pvfun0 <- function(p, x, n){
pu <- pbinom(x-1, n, p, lower.tail = FALSE)
pl <- pbinom(x, n, p, lower.tail = TRUE)
2*pmin(pl, pu)
}
sol0 <-uniroot.all(f = function(p){pvfun0(p,5,10)-0.05},
interval = c(0.1,0.9))
res_b <- binom.test(5, 10)
print(res_b$conf.int)
@abikoushi
abikoushi / confint_lr.R
Created August 7, 2024 07:33
Example: Getting confidence interval numerically
library(rootSolve)
LRtest <- function(p,x,n){
phat <- x/n
lp <- 2*(dbinom(x, n, phat ,log = TRUE)-
dbinom(x, n, p ,log = TRUE))
pchisq(lp, 1, lower.tail = FALSE)
}
set.seed(123)
n <- 20
x <- rbinom(1, n, 0.5)
@abikoushi
abikoushi / pvalfun_prop.R
Last active August 6, 2024 14:19
P-value function of prop.test
set.seed(1234)
n <- 20
x <- rbinom(1,n,0.5)
print(x)
probbeta <- function(p,x,n){
ahat <- x + .5
bhat <- n - x + .5
pl <- pbeta(p, ahat, bhat, lower.tail=TRUE)
pu <- pbeta(p, ahat, bhat, lower.tail=FALSE)
@abikoushi
abikoushi / check_prop.R
Created August 6, 2024 03:13
Checking prop.test (confidence interval of the proportion)
set.seed(1234)
n <- 12
x <- rbinom(1,n,0.5)
res_p <- prop.test(x, n, p = 0.5, correct = FALSE)
CI_score <- function(x, n, level){
z <- qnorm(0.5*(1-level), lower.tail = FALSE)
phat <- x/n
t_1 <- phat+(z^2)/(2*n)
t_2 <- sqrt(z^2+4*n*phat*(1-phat))*(z/(2*n))
@abikoushi
abikoushi / miniLineTable.R
Created July 19, 2024 09:20
sparkline with dendrogram
plot.miniLineTable <- function(tab,...){
oldpar <- graphics::par(no.readonly = TRUE)
N <- length(tab$y)
graphics::par(mar = c(4, 10, 0, 2), oma = rep(1, 4))
graphics::layout(mat = cbind(N+1,1:N,1:N), respect = FALSE)
for (i in 1:N) {
tmpy <- tab$y[[i]]
matplot(tmpy, type = "l", lty=1, xaxt="n", yaxt="n",
xlab = "", ylab = "", frame.plot = FALSE,
col=1)