行列の指数関数と常微分方程式について入門するためのちょっとしたメモです.独学者が行列の指数関数にいきなり出会ってしまっても怯まないように,という願いが込められています.
例によってあまり self-contained な書き方でないので文献を紹介しておきます.行列の指数関数と微分方程式については
| library(ggalluvial) | |
| library(dplyr) | |
| titanic_wide <- data.frame(Titanic) | |
| titanic1 = group_by(titanic_wide,Class, Age, Survived) %>% | |
| summarise(Freq = sum(Freq)) %>% | |
| ungroup() %>% | |
| mutate(dummy=1L) |
| library(dplyr) | |
| library(ggplot2) | |
| #see following | |
| #browseURL("https://rekihaku.pref.hyogo.lg.jp/curator/20165/") | |
| wazauta = scan("wazauta.txt", what = character()) | |
| n = nchar(wazauta) | |
| unigram = character(n) | |
| for( i in 1:n ){ | |
| unigram[i] = substr(wazauta, i, i) |
| 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)) %>% |
行列の指数関数と常微分方程式について入門するためのちょっとしたメモです.独学者が行列の指数関数にいきなり出会ってしまっても怯まないように,という願いが込められています.
例によってあまり self-contained な書き方でないので文献を紹介しておきます.行列の指数関数と微分方程式については
| 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) |
| library(BiasedUrn) | |
| softmax <- function(x){ | |
| mx = max(x) | |
| ex = exp(x-mx) | |
| ex/sum(ex) | |
| } | |
| MVNCH_Fisher <- function (x, m, weight) { |
| 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) |
| 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() |
| 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)) |
| 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) |