Skip to content

Instantly share code, notes, and snippets.

@abikoushi
abikoushi / sumouterprod.cpp
Created October 10, 2024 05:35
Sum of the all of the outer product (for Rcpp)
#include "RcppArmadillo.h"
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;
// [[Rcpp::export]]
double sumouterprod(const arma::field<arma::vec> & V){
int K = V.n_rows;
arma::vec tout = V(0);
for(int j=1; j<K; j++){
@abikoushi
abikoushi / plot_tweedie.R
Created September 28, 2024 21:38
Plot tweedie distribution
library(tweedie)
library(dplyr)
library(ggplot2)
library(ggh4x)
pow <- c(1.01,1.05,1.3,1.9)
phi <- c(0.2,0.5,1,3)
mu <- c(0.1,0.3,2)
dens <- expand.grid(y=y,power=pow,mu=mu,phi=phi) %>%
@abikoushi
abikoushi / geom.tex
Created September 5, 2024 06:39
Tikz example: relationship between binomial and geometric distribution
\documentclass{beamer}
\usepackage{tikz}
\setbeamertemplate{navigation symbols}{}
\usetikzlibrary{positioning}
\makeatletter
\usefonttheme{serif}
\usepackage{xcolor}
\def\mathunderline#1#2{\color{#1}\underline{{\color{black}#2}}\color{black}}
@abikoushi
abikoushi / Fishertest.r
Created August 31, 2024 05:44
p-value function vs. posterior probability
library(BiasedUrn)
library(exact2x2)
library(animation)
#MCMCpack::dnoncenhypergeom(x = NA, cs[1],cs[2],rs[1], 1)
saveGIF({
for(i in c(0:17,17:0)){
rs <- c(18,17)
cs <- c(17,18)
X <-matrix(c(i,cs[1]-i,rs[1]-i,cs[2]-(rs[1]-i)), nrow=2)
@abikoushi
abikoushi / cplm.R
Created August 29, 2024 08:20
linear model with "change point"
x <- seq(-1,1, length.out=100)
x2 <- 1*(x>0)
X <- cbind(x, x2, x*x2)
B <- c(0.5,2,-1)
mu <- c(X%*%B)
set.seed(2024);y <- mu + rnorm(100,0,0.1)
fit <- lm.fit(X,y)
muhat <- c(X%*%c(fit$coefficients))
@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