Skip to content

Instantly share code, notes, and snippets.

@abikoushi
abikoushi / shimarisu_ch1_repro.R
Created September 3, 2025 10:59
佐藤俊哉『宇宙怪人しまりす統計よりも重要なことを学ぶ』第1話の最初の表の数値例
########
## 佐藤俊哉『宇宙怪人しまりす統計よりも重要なことを学ぶ』(朝倉書店)
## 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
abikoushi / HDI.R
Created September 3, 2025 06:53
Highest density interval of beta distribution
# 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
abikoushi / multiple_comparisons.R
Created September 2, 2025 16:07
Plot p-value function with multiple comparisons
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
abikoushi / pvalfun_prop.R
Created August 30, 2025 08:14
animate p-value function of the binomial distribution
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
abikoushi / causalDAG.R
Created August 26, 2025 08:05
Plot Causal Diagram using DiagrammeR
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
abikoushi / qrviz.R
Created August 23, 2025 15:39
Visualize QR algorithm for covariance matrix
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)
@abikoushi
abikoushi / SliceSampling.jl
Last active September 3, 2025 08:09
Julia adaptation Mochihashi, Daichi (2020) Unbounded Slice Sampling
# Mochihashi, Daichi (2020) Unbounded Slice Sampling
# https://arxiv.org/abs/2010.01760
using Distributions
using Random
using Plots# Mochihashi, Daichi (2020) Unbounded Slice Sampling
# https://arxiv.org/abs/2010.01760
module Slice
# Unbounded slice sampling
function slice1(x, likfun; A=100.0, maxiter=1000)
# shrink/expand transformation
@abikoushi
abikoushi / entropy.md
Last active August 18, 2025 09:02
確率の対数を取る(あるいは情報エントロピー入門)

推敲して数値例も足したので以降はこっちを見てください→ 確率の対数を取る(あるいは情報エントロピー入門)

エントロピー

確率の話でよく一番基本的な例として出てくる「公平なコイン投げ」を考える.表が出る確率は $1/2$ である. $k$ 回のコイン投げで $k$ 回表が出る確率は $2^{-k}$ だ.

コイン投げの例がわかりやすいとしたら,逆に確率 $p$ をコイン投げでいうと何回表が出る程度の珍しさかで比べるという発想もありうる.

2 を $-k$ 乗したら $p$ になるような数 $k$$2^{-k} = p$ となる $k$ )は,対数の定義から次のように表せる.

@abikoushi
abikoushi / phi_and_chisq.md
Created August 12, 2025 04:11
分割表の独立性の検定と相関係数

このノートについて

分割表の独立性の検定について,「独立性というからには相関係数の話と関係あるのかな?」と思った人向けの文章です. 分割表の独立性の話と相関係数の話は関係があります.

地道手計算パート

このパートの式の部分は計算メモのようなものも含むので,もしかしたら冗長な書き方が多いかもしれない.そのように感じたら結論まで読み飛ばしてほしい.

@abikoushi
abikoushi / shimarisu_ch4_strata.R
Last active July 21, 2025 01:37
層別解析:佐藤俊哉『宇宙怪人しまりす 統計よりも重要なことを学ぶ』(朝倉書店)第4話より
library(DiagrammeR)
library(DiagrammeRsvg)
library(rsvg)
library(dplyr)
library(gt)
## 佐藤俊哉『宇宙怪人しまりす 統計よりも重要なことを学ぶ』(朝倉書店)より
g <- grViz("digraph{
graph[rankdir = TB]
node[shape = rectangle]