abikoushi

# https://en.wikipedia.org/wiki/Wallenius%27_noncentral_hypergeometric_distribution#Multivariate_distribution
# https://markheckmann.github.io/notes/wallenius-distribution.html
library(dplyr)
library(tibble)
library(BiasedUrn)
set.seed(1234); R = rexp(10)

ressim = t(replicate(1e5L,sort(sample.int(10, size = 3, prob = R))))
pmf = data.frame(ressim) %>% 
  group_by_all() %>% 

abikoushi

\documentclass[tikz, border=14pt]{standalone}
\usepackage{tikz}
\usetikzlibrary{bayesnet}
\begin{document}
\begin{tikzpicture}
    \node[obs] (x) {$X$};
    \node[obs, above right = of x] (z) {$Z$};
    \node[obs, below right = of z] (y) {$Y$};
    \edge{x}{y};
    \edge{z}{y};

abikoushi

#browseURL("https://en.wikipedia.org/wiki/Weibull_distribution")

median_weibull <- function(shape, scale=1){
  ln2 = log(2)
  scale * (ln2^(1/shape))  
}

relativepoverty <- function(shape, scale=1){
  pweibull(0.5*median_weibull(shape, scale), shape, scale)
}

abikoushi

########
## 佐藤俊哉『宇宙怪人しまりす統計よりも重要なことを学ぶ』（朝倉書店）
## Chapter 1: 再現性
########
library(dplyr)
library(ggplot2)
TPratio = function(true_ratio,
                           beta = 0.8,
                           alpha = 0.05){
  study1 =cbind(true_ratio*c(beta,1-beta), 

abikoushi

########
## 佐藤俊哉『宇宙怪人しまりす統計よりも重要なことを学ぶ』（朝倉書店）
## Chapter 1: 再現性
########
library(dplyr)
library(ggplot2)
reproducibility = function(true_ratio,
                           beta = 0.8,
                           alpha = 0.05){
  study1 =cbind(true_ratio*c(beta,1-beta), 

abikoushi

# the following are used as reference 
# https://stats.stackexchange.com/questions/381520/how-can-i-estimate-the-highest-posterior-density-interval-from-a-set-of-x-y-valu

hdi_beta = function(shape1, shape2, coverage, n){
  x = seq(0,1, length.out=n)
  best = 0.0
  for (ai in 1 : (length(x) - 1)){
    for (bi in (ai + 1) : length(x)){
      mass = pbeta(x[bi], shape1, shape2) - pbeta(x[ai], shape1, shape2) 
      if (mass >= coverage && (mass / (x[bi] - x[ai])) > best){

abikoushi

library(ggplot2)
library(tidyr)
library(dplyr)

X = matrix(0,10,10)
set.seed(1234)
X[,1:5] = rnorm(50, 0)
X[,6:10] = rnorm(50, 2)

colwise_ttest <- function(X, mu, method=NULL){

abikoushi

library(ggplot2)
library(dplyr)
library(animation)
library(rootSolve)

probbeta <- function(p, x, n, a=0.5, b=0.5){
  ahat <- x + a
  bhat <- n - x + b
  pl <- pbeta(p, ahat, bhat, lower.tail=TRUE)
  pu <- pbeta(p, ahat, bhat, lower.tail=FALSE)

abikoushi

library(DiagrammeR)
library(DiagrammeRsvg)
library(rsvg)
#https://elevanth.org/blog/2021/06/21/regression-fire-and-dangerous-things-2-3/
g <- grViz("digraph{
                     node[shape = plaintext]
                     U[label = 'U', shape=circle]
                     B1[label = 'B1']
                     B2[label = 'B2']
                     M[label = 'M', fontcolor='red']

abikoushi

library(mvtnorm)
library(reshape2)
library(ggplot2)
library(gganimate)
library(dplyr)

Sig = matrix(c(1,0.8,0.8,1),2,2)
scatters = vector("list", 5)
heats = vector("list", 5)
X = rmvnorm(5000, sigma=Sig)