abikoushi

library(janeaustenr)
library(tidytext)
library(dplyr)
library(ggplot2)

book_words <- austen_books() %>%
  unnest_tokens(word, text) %>%
  count(book, word, sort = TRUE) %>% 
  group_by(book) %>% 
  mutate(total = sum(n)) %>%

abikoushi

行列の指数関数と常微分方程式についてのイントロ

はしがき

行列の指数関数と常微分方程式について入門するためのちょっとしたメモです．独学者が行列の指数関数にいきなり出会ってしまっても怯まないように，という願いが込められています．

例によってあまり self-contained な書き方でないので文献を紹介しておきます．行列の指数関数と微分方程式については

1. 長谷川浩司 線形代数
2. Hirsch・Smale・Devaney 力学系入門

abikoushi

library(deSolve)
modLinear <- function(Time, State, A) {
  return(list(A%*%State))
}

a = 0.2
b = 0.1
A = matrix(c(-a, 0,
              a, -b), 2, 2, byrow = TRUE)
yini  <- c(X = 1, Y = 0)

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)