行列の積の計算は%*%を使う。 > a <- matrix(1:4, 2, 2) > a [,1] [,2] [1,] 1 3 [2,] 2 4 > b <- matrix(0:3, 2, 2) > b [,1] [,2] [1,] 0 2 [2,] 1 3 > a+b [,1] [,2] [1,] 1 5 [2,] 3 7 > a%*%b [,1] [,2] [1,] 3 11 [2,] 4 16 > a*b #これは要素ごとの積…

RでDFAのチェック。 set.seed(1) #White noise eda.plot(DFA(rnorm(1500))) #H~0程度 #Random walk eda.plot(DFA(cumsum(rnorm(1500)))) #H~0.5程度お手軽ランダムウォークはただ単にcumsum(rnorm(1500))，でよろしい。相関をチェックしたければacf(rnorm(15…

先日読んだ論文で使われていたDFAを試してみようとRのpackageをインストール。 install.packages("fractal", dependencies = TRUE) library(fractal) #正規乱数で試して見る x <- DFA(rnorm(1024)) print(x) eda.plot(x) #正規乱数で試して見る，時間窓は1.1…

お手軽にデータサンプリングを行う． 読み込んだデータ(m行n列)から，50行だけサンプルする． > set.seed(1) > conjunction <- read.table("ave_c++.txt",header=T) > n_con <- sample(nrow(conjunction),50) > write.table(conjunction[n_con,],file="sampl…

Ubuntuのバージョンを確認 $ cat /etc/lsb-release DISTRIB_ID=Ubuntu DISTRIB_RELEASE=14.04 DISTRIB_CODENAME=trusty #trustyとうコードネーム DISTRIB_DESCRIPTION="Ubuntu 14.04.1 LTS" $ arch x86_64 #64bit次にダウンロード先の指定をしておく． $ sud…

日経225*1のデータをRからゲット． install.packages(quantmod) library(quantmod) getSymbols('^N225') head(N225) #ファイル確認 sink("day.tmp") N225$N225.Open sink() write.table(N225,file="prices.tmp",append=FALSE,quote=FALSE,sep=" ")どういうわ…

applyの引数２つ目が「１」か「２」かで行(=1)か列(=2)かが決まる． > tmp <- matrix(1:10,ncol=2,nrow=5) > tmp [,1] [,2] [1,] 1 6 [2,] 2 7 [3,] 3 8 [4,] 4 9 [5,] 5 10 > apply(tmp,1,sample) #行方向のシャッフル [,1] [,2] [,3] [,4] [,5] [1,] 1 2 8…

読み込んだファイルをお手軽ランダムシャッフル． x0 <- read.table("Forward_Number_Sherlock.txt") x <- x0[,1] set.seed(10) y <- sample(x) write.table(y,"Shuffle10_Number_Sherlock.txt",append=FALSE,row.names=FALSE,col.names=FALSE,quote=FALSE)…

単位根検定 library(tseries) a <- PP.test(x) b <- adf.test(x) ave <- round(mean(x),digits=5) std <- round(sd(x),digits=5) a2 <- PP.test(z) b2 <- adf.test(z) ave2 <- round(mean(z),digits=5) std2 <- round(sd(z),digits=5) out <- paste(commandA…

1列目と5列目のデータのみを取り出す． ついでに，中央値以上の行のみにしてから取り出している． max <- length(z1) med_z <- z1 n0 <- (1:max)[z1>med_z] up <- x[c(n0),c(1,5)]

Rstudioで日本語がプロットできない．.Rprofileというファイルを~/ 以下に作って，以下に記述． Rstudioを再起動したら，出来た．．． こちらのサイトの内容の丸写しです．ありがたやー setHook(packageEvent("grDevices", "onLoad"), function(...){ if(.Pl…

ベクトルから中央値を計算し，中央値の値だけ取り除く． x[c(-1)]で，x[1]の要素を取り除くという意味． med_x <- median(x) n0 <- (1:length(x))[x==med_x] #if(length(n0)>1) n0 <- n0[1] もし要素を1つだけ取り除く場合 x <- x[c(-n0)]

読み込んだ距離行列から，デンドログラムを作成する． リスト→マトリクス形式にawkでファイルを変換した後行う． rm(list=ls()) y2 <- read.table("Convert_Matrix_Count_v1.txt", sep=",") y <- 1/y2 #大きい方が距離近い数なので逆数に y[upper.tri(y,diag…

Rでガンマ関数を使ったFisher's exact testをしていたのだが， ガンマ関数が200!とかの階乗の計算では発散して，出来なくなった． 困ったなぁと思っていたら，Rで2×2の分割表をそのまま読み込んで， Fisher's exact testが出来ることが分かった． > c1 <- ma…

Fisher's 精密検定のため，階乗の計算をするのに，gamma関数を使うが，gamma(200)とかすると， >gamma(200) [1] Inf 警告メッセージ： 'gammafn' 中の値が範囲を超えています 範囲を超えるらしく計算してくれない．Fisher's testであれば，変数が多い場合は…

ベクトル要素の最後に数字と付け足す． w=NULL for(i in 1:10){ v=1 w = append(w, v,after=i) }他に何かいい方法がありそうなのだけれど，とりあえずこれでも処理可能．

Rで階乗の計算は，例えば5!ならば以下． gamma(5+1)普通階乗のn!におけるnは非負整数だったが，これを複素数含む一般の形にしたものがガンマ関数． n!=gamma(n+1)だから，0!=gamma(1)=1 と考える．ガンマ関数の定義式(x>0)は，18世紀(1729?)にオイラーによっ…

Rで自己相関関数を出力． namesで確認して，必要な値~$acfのみをファイルに出力する． FN <- paste("Acf_",commandArgs()[5],sep="") ax <- acf(x[,2],lag.max=60) write.table(ax$acf,FN,append=F,quote=F,col.names=F,row.names=F)

zにNAが複数含まれているので，普通にmean(z)としても平均が算出できない． そこで，na.omit(z)を使って単純にNAを取り除く． ただし，summaryだと，meanを計算してくれる模様． z1 <- na.omit(z) summary(z1) Min. 1st Qu. Median Mean 3rd Qu. Max. -1.630…

warnOnlyだと，nlsの返り値がないため，summaryするとエラーが出る． そこでtry関数を使って，とりあえず値があるかチェックしておいて， class( )で内容に"try-error"が含まれていたら，ループを飛ばすなどする． for(j in j1:jn){ y2 <- c[(n+2):j] x2 <- …

nlsで回帰をしていると，値が収束しなくて， 「繰り返し数が最大値 50 を超えました」というエラーで処理がストップすることがある． そういう場合は，controlオプションを使って，警告のみをTRUEにしておく． n <- i y1 <- c[1:n] x1 <- seq(1,n,1) dat <- …

サンプルから1つずつ選んで，戻して，ベクトルxに付け加えていく． rm(list=ls()) dat <- read.table("WordPool_Neko.txt") set.seed(1) x <- as.numeric(NULL) #random shuffle #x <- sample(dat[,1]) #word pool simulation for(i in 1:length(dat[,1])){ …

ベクトルの何番目の要素が条件を満たしているかを調査． checkのベクトル始めはすべて0で，ある要素だけ，1になっている． その要素の番号を知りたい例． N <- 100 # check flg check <- as.vector(0) for (i in 1:N) check[i] <-0 check[10] <- 1 check[20] …

組み合わせ(Combination)はchoose(n,k)を使う．[=nCk]確率pの事象が実現される確率が，n回中k回以下である確率は以下． rm(list=ls()) x <- 0 p <- 0.5 n <- 1000 k <- 500 sum_k <- c(0:k) y <- choose(n,sum_k)*p^sum_k*(1-p)^(n-sum_k) plot(y,log="y") s…

unique関数を使い，Heaps則をRで描く. uniqueは与えられたベクトルから，重複しない要素だけ抜き出すので， uniqueで抜き出した要素数をlengthで数えれば，Heaps則を描くことができる． ただし，時間がかかる． N <- 20000 DW <- 5000 doc <- sample(DW,N,re…

リストやベクトルの初期化をしておかないと，エラーが出る場合がある． [1,10]のランダムシャッフルの例．(別に必要ないけど．) N0 <- 1:10 N1 <- as.vector(NULL) N1 <- sample(N0)

行列計算を使って高速化． 要素が数字の配列a[i]を使うときには，ループの前に初期化しておくと処理が早い． a でももっと早いのはループ(for文)を使わずに，ベクトル単位で処理すること． - #Loop for(i in 1:n) a[i] #Vector unit a - 計算結果は同じでも…

igraphを使った描画． 無向，リンク数に比例したpreferential attachmentで成長するノードN=100個のBAモデルのネットワークを作成し，描画． 描画の際にはリンクにはラベルを書かないようにするには，「vertex.label=NA」が必要． N <- 100 g <- barabasi.ga…

igraphを使った解析．igraphはコマンドに重複名がある． 例えば，スケールフリーネットワークを作る，barabasi.gameとba.gameは同じ． degreeとdegree.distributionはヒストグラムと，確率密度分布だが． barabasi.gameでできるBAモデルのデフォルトはリンク…

RStudioのインストールと，仕切り直し． 非公式日本語ガイド 図の保存などが便利そう． コマンドラインからワークスペースを移動するには setwd ("~/Desktop/") getwd ()とした． せっかくなので，Rのコーディングにおけるgoogle流のルールに習って書くこと…