忍者ブログ
物理学者(ポスドク)による日々の研究生活のメモ書きです ( python/emacs/html/Japascript/シェルスクリプト/TeX/Mac/C言語/Linux/git/tmux/R/ポケモンGO)
dir_r_result = sprintf("r_%s", header)
dir.create(dir_r_result)

library(fs)
dir_create(dir_r_result)
みたいな感じ

2種類あるけど、どちらもディレクトリを作れる
dir.create()はすでにディレクトリがあると、エラーで止まる

fsというパッケージの方のdir_create()はすでにあると何もエラーが出ない
後者の方が便利そう





PR

■ あるベクトルに他のベクトルを追記したいとき
x <- 1:10="" br="">x = 1:10
y = 11:20
x
[1] 1 2 3 4 5 6 7 8 9 10
y
[1] 11 12 13 14 15 16 17 18 19 20
x = c(x, y)
x
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20





■ for文で回した計算をベクトル(リストではない)に追加していく
# まず空のベクトルを作る
list_sim = c()

# n_simがシミュレーション回数
for (i_sim in 1:n_sim){
 
 # ここでなんか計算する
 ....
 # 計算結果をベクトルに追加したい
 list_sim = c(list_sim, result)

}



タイトルにある通り、ここではベクトルを使っている

最初ハマったのはこういうときはリストを使うんだろうなぁ〜
と思って「R リスト 追記」とかで調べてしまった・・・

そしたら
list_sim = list()

# n_simがシミュレーション回数
for (i_sim in 1:n_sim){
 
 # ここでなんか計算する
 ....
 # 計算結果をベクトルに追加したい
 list_sim = append(list(list_sim), result)

 とか

 list_sim = c(list(list_sim), result)

}

とか出てくるが、なんか結果が自分のほしいものとは違う・・・という罠にハマった










#fitdiftrを使うためのパッケージ
library(MASS)

# tableでデータ読み込み、1列目のみ抜き出しているが必要かはわからん
data1 <- read.table("hoge.txt")
data2 <- data1[, 1]

# chi square分布でフィットする
# chi_kには推定された自由度を代入しておく
chi_df <- fitdistr(data2, "chi-squared", start=list(df=2), method="BFGS")
chi_k <- chi_df[[1]][1]

# histではなくtruehistで確率密度関数にする
chi_hist <- truehist(data2, h=0.1, prob = TRUE, col="gray", main="hogehoge")
curve(dchisq(x, df=chi_k),add=TRUE,col="red", lwd=3)







■ 2つのデータに対してKS検定をする場合
% x=rnorm(10, mean = 1)
% y=rnorm(10, mean = 1)
% x
[1] 0.91430534 1.59262261 1.52565219 1.68065990 2.35755910 1.91386531
[7] -2.26359817 0.08465784 0.06163455 -0.10943752
% y
[1] 2.2200235 2.2241460 1.0733582 0.8002340 1.6228854 0.2652880
[7] 1.4286084 -0.6869875 0.1918104 2.2083442
% ks.test(x, y)

Two-sample Kolmogorov-Smirnov test

data: x and y
D = 0.3, p-value = 0.7869
alternative hypothesis: two-sided


■ 1つのデータに対した正規分布のKS検定をする場合
% x=rnorm(10, mean = 1)
% ks.test(x, pnorm)

One-sample Kolmogorov-Smirnov test

data: x
D = 0.43645, p-value = 0.02996
alternative hypothesis: two-sided

2つの引数は、比較したい分布を指定する文字列、正規分布なら pnorm、ポアソン分布なら ppois、一様分布の場合 punif を指定できる

■ 正規分布の場合だと、平均と分散も指定可能
% ks.test(x=vx,y="pnorm",mean=mean(vx),sd=sd(vx))

One-sample Kolmogorov-Smirnov test

data: vx
D = 0.13793, p-value = 0.9196
alternative hypothesis: two-sided


■ 片側検定/両側検定の切り替えはalternativaで行う
デフォルトでは両側検定になっている
% ks.test(x, y, , alternative="t")

Two-sample Kolmogorov-Smirnov test

data: x and y
D = 0.3, p-value = 0.7869
alternative hypothesis: two-sided

% ks.test(x, y, , alternative="l")

Two-sample Kolmogorov-Smirnov test

data: x and y
D^- = 0.1, p-value = 0.9048
alternative hypothesis: the CDF of x lies below that of y

% ks.test(x, y, , alternative="g")

Two-sample Kolmogorov-Smirnov test

data: x and y
D^+ = 0.3, p-value = 0.4066
alternative hypothesis: the CDF of x lies above that of y



■ 検定結果からp値を取り出したいときは、
結果をなにかに代入して、namesで中身を調べる→ p値のみを取り出す
% out <- ke.test(x, y)
% names(out)
[1] "statistic" "p.value" "alternative" "method" "data.name"
% out$p.value
[1] 0.7869298







rm(list=ls())

■ 参考 : Rで使わない変数を消す




% .libPaths()
でわかる

自分の場合は
/Library/Frameworks/R.framework/Versions/3.5/Resources/library
だった


自分で追加したい場合は
% .libPaths("./hoge")
のようにする
これを ~/.Rprofile に書いておくと勝手に起動時に読み込まれる
(パッケージのインストール先の優先順位はどうなる? )





■ plotの基本
plot(x, y)
■ 参考 : グラフィックスパラメータを設定する、確認する


■ 対数表示、常用対数のみ使える
plot(x, y, log="x")
plot(x, y, log="y")
plot(x, y, log="xy")


■ 範囲指定
plot(x, y, xlim=c(0, 1))
plot(x, y, ylim=c(0, 1))


■ タイトルとかラベル
plot(x, y, main="上につくタイトル")
plot(x, y, xlab="x軸のラベル")
plot(x, y, ylab="y軸のラベル")
ラベルとかを大きくするにはtmag=1.5とかするらしいけど、現在のRではもはやサポートされていないらしい
今使えるのは↓
mainを大きくしたいときは cex.main=1.5
labelを大きくしたいときは cex.lab=1.5
数字を大きくすると大きくなる


■ 重ね書き
plot(x, y)
par(new=T)
plot(a, b)
とか、linesを使う場合はpar(new=T)はいらないっぽい
ググった感じ、そっちを奨励してる人もいた
あとpar(new=T)を使う場合はグラフの範囲指定を意識的にしないといけない
plot(x, y)
lines(a, b)


■ 線の色
col = 1
col="blue"
col = rgb(1, 0, 0)


■ 線の太さ
lwd=2





(2019/03/06 追記)
# rhoをギリシャ文字で入れたい場合は
main = expression(rho)

# 通常の文章にギリシャ文字を混ぜたい場合はexpressionの中でpasteを使う
main = expression(paste(lambda, "=5 のポアソン分布 "))

■ 参考 : ggplot2でギリシャ文字や数式を表示したい

■ 参考 : Rコード:プロットのタイトルにギリシャ文字を利用する。







ここがわかりやすかった
といってもヒストグラムを1から書いたわけではなくて、どういうオプションが使えるのかをちょっと知りたかっただけ
■ 参考 : Rによるヒストグラムの描画



■ グラフの中に最尤推定法で得られたパラメーターを書き込みたい

■ 参考 : Rの数式描画の中で変数を展開したいときはbquoteを使う
hist(dat, freq = FALSE, ylim = c(0, max(max(h$density),max(y))), xlab = "z", ylab = "f(z)",
main = paste("Density Plot\nmu = ", round(a[1], 3), "\ntheta = ", round(a[2], 3), "\ngamma = ", round(a[3], 3) ), cex.main=1 )
cex.mainでmainの文字サイズを指定できる

pasteで文字列と変数を繋げた文字列を作ることができる
これで、グラフの中に最尤推定法で得られたパラメーターを書き込める

pasteでつなげたときに間に空白を入れたくないときは
paste("AAA", "BBB", sep="")
round(x, 2)
は桁数を短くするための関数
これで出力する小数点以下の数を2桁にできる







(2019/08/02 追記)

とりあえずRで困ったらhelp(hist)って打てば大体わかるから

ヒストグラムのbin数を変えたくてあれこれ調べた
ビン数を30にしたかったら、break=30

y軸をログにするには、result = hist()で結果を一度格納してから、そのresultのlogを取るしかないっぽい


↓もっと詳細なヘルプは実際にRで見てください
hist package:graphics R Documentation

Histograms

Description:

The generic function ‘hist’ computes a histogram of the given
data values. If ‘plot = TRUE’, the resulting object of class
‘"histogram"’ is plotted by ‘plot.histogram’, before it is
returned.

Usage:

hist(x, ...)

## Default S3 method:
hist(x, breaks = "Sturges",
freq = NULL, probability = !freq,
include.lowest = TRUE, right = TRUE,
density = NULL, angle = 45, col = NULL, border = NULL,
main = paste("Histogram of" , xname),
xlim = range(breaks), ylim = NULL,
xlab = xname, ylab,
axes = TRUE, plot = TRUE, labels = FALSE,
nclass = NULL, warn.unused = TRUE, ...)






Rを使うことになったので新しくカテゴリを作った
ヘビーユーザーにはならないと信じてるので、凝ったことはしない




■ インストール
Mac用のRは https://cran.r-project.org の「Download R for (Mac) OS X」からインストールできる

追加のパッケージをインストールしたいときは起動後に
パッケージとデータ → パッケージマネージャ
からインストールできる


R起動時に毎回パッケージを読み込むようにするのがめんどくさいので設定ファイルを用意する
~/.Rprofile を作成する
読み込みたいパッケージ名が hoge だとすると、Rprofileの中に
library(hoge)
と書いておけばOK


RをGUIで実行したときもCUIで実行したときもその履歴を残したいときは、.Rprofileに次の文章を追記しておく
.Last = function() { if (interactive()) try(savehistory()) }





■ はじめに読む

■ 参考 : Rの初歩





■ コマンドメモ
このあと追記していきます
getwd()
今いる場所を調べる

setwd("/Users/hoge/foo")
ワーキングディレクトリを変える





■ データを読み込みたい

こんなデータを読み込みたい
1 5.2
2 2.3
3 1.1
4 6.9
5 10.2
6 1.1
7 2.4
みたいな感じで、複数行・複数列のデータ

↓ ここには情報はなかった
■ 参考 : Rでデータ読み込みから前処理までのTips


Rのwikiがあったので見てみたが、↑のようなフォーマットのデータの読み込み方法は書いてなかった
■ 参考 : ファイルを読み込む tips 集(暫定版)


ググってるとdata.frameというのを使うのが良さそうだった
data.frameについてよくまとめられていた
■ 参考 : データフレームTips大全
# テキストデータをその構造を維持したまま読み込み
# 列(変数名)ラベルが自動的につけられる、これがいらない場合は read.table("hoge.txt", header=TRUE)とすればいい?
data1 <- read.table("data/hoge.txt")


# 1列目のみ欲しいとき
data2 <- data1[, 1]


# data2のうち範囲指定してデータを取り出したいとき
# この場合は 4.6 < x < 5.0 を満たすデータのみ取り出せる
data3 <- subset(data2, 4.6<data2, data2<5.0)



■ スクリプト的にRを使いたいとき

やることを hoge.R というファイルに書いておく
Rを起動して
source("hoge.R")
で書いてあることを順次実行してくれる



■ データ書き出し

■ 参考 : テキストデータのファイル出力に関係する関数の例

writeを使うのが一番簡単そう
writeでは1つの変数しか書き出すことができないので、pasteを使って複数の変数を1つの文字列に収める必要がある
それを何度も書き出せばいいっぽい
Rだとヘビーな計算はさせないやろうし



■ whileループ
i <- 1
while (i < 6) {
print(i)
i = i+1
}
■ forループ,
for (i in 1:5) {
 print(i)
}
■ グラフを画像として保存したい
png("hoge.png")
のあとでプロット
画像の書き込みを終了するときは
dev.off()

pdfにしたいときは
pdf("hoge.pdf")
このあとでplotを何個も書くと、その分だけページができて複数の画像を同時に保存できる

epsは
postscript("hoge.eps")
tex用のepsは
postscript("hoge.eps", horizontal = F, onefile = F, paper = "special")


gnuplotでいうところの
set terminal something
set output "somthing"
を同時にできる感じなので便利





(2019/03/01 追記)
■ 指数表示にしたいとき

「R 指数表示にしたい」でぐぐったら指数表示にしない方法しか出てこなかった
とりあえず sprintf("%e", a)って感じで対応した
sprintf("%.3e", a)とかC言語と同じフォーマットはそのまま使える
便利




(2019/03/07 追記)

xという配列にデータ点を追加したいとき
先頭に1を追加したいときは
new_x < c(1, x)

後ろに1を追加したいときは
new_x < c(x, 1)




(2020/09/08 追記)

■ シェルスクリプト上でR

単純にシェルスクリプトの上でRを呼び出したら以下のようなエラーが出た
Fatal error: you must specify '--save', '--no-save' or '--vanilla'
zsh: exit 2

以下のようにオプションを付けたらエラーが出なかった
#/bin/sh

r --quiet --vanilla <<EOF
source("./r/hoge.r")
EOF










HOME |
プロフィール
HN:coffee
職業:物理屋(自称)
趣味:映画鑑賞、登山
出身:大阪府の南の田舎
自己紹介:
import MyProfile
import coffee_pote from TWITTER
import amazonのほしい物リスト from WISH_LIST

print "先月子供が産まれました!"

# 最終更新 2022/10/25
カウンター
カウンター カウンター
ブログ内検索
ツイートするボタン
リンク
相互リンク募集中です (Twitterにてお知らせください)

Demo scripts for gnuplot version 5
(gnuplotのさまざまなデモ画像と作り方がまとめられている、眺めているだけでできるようになった気分になれる)

gnuplotスクリプトの解説
(米澤進吾さんの個人ページ、gnuplotと言えばこのかた)

gnuplot のページ
(Takeno Lab、うちのブログがリンクされていたのでリンク返し)

Twitterから映画の評価が分かる & 映画の鑑賞記録が残せる coco
(映画の感想をまとめられるサイト、いつもお世話になっています)

Astronomy Picture of the Day Archive
(天文や宇宙関連の最新の話題について画像とともにNASAが説明しているページ)

今日のほしぞら
(任意の時刻の空で見える星を表示してくれる、国立天文台が管理している)

GNUPLOTとアニメーション
(応用の項目の「見せてあげよう!ラピュタの雷を!!」あたりからすごすぎる)

読書メーター
(読んだ本をリストできる便利なサイト)

flickr難民の写真置き場
(20XX年、flickrは有料化の炎に包まれた。あらゆるflickr無料ユーザーは絶滅したかに見えた。 しかし、tumblr移住民は死に絶えてはいなかった。)

教授でもできるMac OS X へのLaTeX, X11, gccのインストレーションと環境設定
(阪大の山中卓さんのwebページ、タイトルにセンスが溢れている、内容は超充実してる、特にTeX関連、学振DCとかPDの申請書類作成時にはお世話になっております)

英語論文執筆用の例文検索サービス
(とんでもないものを見つけてしまった・・・・ arXivに収録されている 811,761報の 英語論文から,例文を検索するための検索エンジン)


Template "simple02" by Emile*Emilie
忍者ブログ [PR]