abikoushi

library(BiasedUrn)

softmax <- function(x){
  mx = max(x)
  ex = exp(x-mx)
  ex/sum(ex)
}


MVNCH_Fisher <-  function (x, m, weight) {

abikoushi

library(BiasedUrn)
softmax <- function(x){
  mx = max(x)
  ex = exp(x-mx)
  ex/sum(ex)
}

MVNCH_Fisher <- function(K,M,OR){
  CB = combn(K,M)
  nr = ncol(CB)

abikoushi

library(BiasedUrn)
library(dplyr)
K=10
set.seed(1234); OR = rexp(K)
M=3
ressim = replicate(100000, sort(sample.int(K, size = M, prob = OR, replace = FALSE))) %>% 
  t() %>% 
  as.data.frame() %>% 
  group_by_all() %>% 
  tally()

abikoushi

library(BiasedUrn)
library(dplyr)
K=10
set.seed(1234); OR = rexp(K)
ressim = sample.int(K, size = 10000, prob = OR, replace = TRUE) |>
  table() |>
  as.data.frame()

ressim = mutate(ressim, prob = Freq/sum(Freq))

abikoushi

library(brunnermunzel)
#https://cran.r-project.org/web/packages/brunnermunzel/vignettes/usage.html

## median of the difference: point a s.t. P(X-Y<a)+.5*P(X-Y=a) = 0.5
conv <- function(x){
  #numerical convolution
  integrate(function(y){dnorm(x+y,mu)*dexp(y)},-Inf,0)$value + 
    integrate(function(y){dnorm(x+y,mu)*dexp(y)},0,Inf)$value
}
conv <- Vectorize(conv)

abikoushi

library(DiagrammeR)
library(DiagrammeRsvg)
library(rsvg)
library(dplyr)
library(dqrng)
library(ggplot2)

g <- grViz("digraph{
                     graph[rankdir = TB]
                     node[shape = none]

abikoushi

library(dagitty)
library(ggdag)
library(ggplot2)
library(gridExtra)

g1 <- dagitty("dag{U -> X; U ->Y; X -> Y}")
exposures(g1) <- "X"
outcomes(g1) <- "Y"

plot(g1)

abikoushi

library(ggplot2)

NB_p <- function(p,x,r){
  pnbinom(x, size = r, prob = p,lower.tail = TRUE)  
}

B_p <- function(p,x,r){
  pbinom(r-1, size = r+x, prob = p, lower.tail = FALSE)  
}

abikoushi

library(rootSolve)
findall_q = function(q_ini, a, b, ...){
  dp = dbeta(q_ini, shape1 = a, shape2 = b)
  sol = rootSolve::uniroot.all(function(X)dbeta(X, shape1 = a, shape2 = b)-dp, interval = c(0, 1), ...)
  sort( sol )
}

prob_diff <- function(q, a, b){
  p = pbeta(q, shape1 = a, shape2 = b)
  diff( p )

abikoushi

#install.packages("ggalluvial")
library("ggalluvial")
library("dplyr")
V_row=matrix(rexp(100*10),100,10)
V_col=matrix(rexp(100*10),100,10)

W_row = matrix(rexp(10*5),10,5)
cl_row = apply(V_row%*%W_row,1,which.max)

W_col = matrix(rexp(10*7),10,7)