上図のように R で複数の軸をつけるには、axis の at と labels 引数を用いる。 df <- read.table("/tmp/test") # CTRUNCATE のログとかから切り出したもの par(mar = c(10, 4, 4, 2)) plot(df$V1, df$V2, type="l", main="Wilson plot", lwd=2, ylab="log(I…

R で、検出器上の各ピクセルの分解能を計算しようとした。最初に書いた愚直なコードは以下の通り。 lambda <- 12400 / 7000 # A distance <- 51.5 * 1000 # um pixel_size <- 50 # um buf_size <- 3000 # px buf_origin <- c(buf_size, buf_size) / 2 reso <…

バイナリファイル書き出しの練習として。 writeADSC <- function(filename, img, lambda_in_A=1.0, orgx=ncol(img)/2.0, orgy=nrow(img)/2.0, pixel_size=0.172, camera_length_in_mm=50.0) { width <- ncol(img) height <- nrow(img) conn <- file(filename…

Rcpp と CBFlib の練習も兼ねてやってみた。エラー処理は無視!cbflib-wrapper.cpp: #include <Rcpp.h> #include <cbf/cbf.h> using namespace Rcpp; using namespace std; // [[Rcpp::export]] NumericMatrix readCBF(string filename) { FILE *fh = fopen(filename.c_str(), "r</cbf/cbf.h></rcpp.h>…

画像回転のための準備。 m <- matrix(1:12, 3) nc <- ncol(m) nr <- nrow(m) m t(m)[seq(nc,1),] m[seq(nr,1),seq(nc,1)] t(m)[seq(1,nc),seq(nr,1)]

R で oversampling / density modification による位相回復のデモを作った。実験的に得られた F とランダムな位相から始める。これを実空間に戻して、マスクをかけてからフーリエ変換する。得られた位相(hopefully improved) と実験的 F を組み合わせて、反…

R の substr と substring は、どの引数がベクトル化されているかで異なる。前者は第一引数(対象の text) 、後者は位置(first, last 引数)がベクトル化されている。両方をベクトル化したものはないので、注意が必要だ。 > substr('ABCDEFG', 1:4, 4:7) [1] "…

R で、確率分布から乱数をサンプリングする runif や rnorm の引数はベクトル化されている。 > runif(10, min=c(0, 100), max=c(10, 110)) [1] 7.046310 103.526202 1.910480 103.118722 3.107876 102.799656 1.802613 109.406715 1.984414 [10] 102.594984 …

Rで一様乱数を得るには runif 関数を使う。ヘルプに、 runif will not generate either of the extreme values unless max = min or max-min is small compared to min, and in particular not for the default arguments. http://stat.ethz.ch/R-manual/R-d…

＃説明が雑!R では、関数名だけをタイプするとそのコードが表示される。例えば apply 。オブジェクト指向で多態になっているものについては、UseMethod を呼び出す総称的関数 (generic function) が表示されるだけだ。こういう時は、methods で一覧を取得で…

R の講義で、TRUE と FALSE の略として T や F を使うべきではないと習った。TRUE/FALSE は予約語であり、いわばリテラルだが、T/F は単なるグローバル変数なので、 TRUE <- FALSE # エラー T <- FALSE # 代入できてしまう! pnorm(0.3, lower.tail=T) # ??? …

R のベクトル化の例で、 a <- 1:10001 b <- a[1:10000] + a[2:10001] のように書くと for ループを回すよりも速いというのがよく紹介されているが、メモリ効率はどうなのか。内部的に右辺の2つの項に対応する a のコピーを2つ作っているとしたら、一時的とは…

実はちゃんとやったことがなかった。 a <- list() class(a) <- "TestClass" testfunc <- function(x) {print("test func")} testfunc(a) # test func testfunc.TestClass <- function(x) {print("test func for TestClass")} testfunc.TestClass2 <- functio…

Scientific Programming の講義では、Sweave ドキュメントの形で課題を提出することになっている。それに備えて早速 Sweave を使ってみたのだが、 R の出力のうち、` や ' などの記号が文字化けする。調査の結果、 \usepackage[utf8x]{inputenc} を記載すれ…

講義で「課題には R 3.x.xを使ってください」と言われたので、Mac Book の R を 2.x から 3.0.2（CRAN版)へアップグレード。そうしたら R studio が R を認識しない。"Cannot locate R" というエラーが出る。3.0.x を使うには、R studio も新しいものである…

R でバイナリファイルを読む練習をした。クロマトグラフィーのデータや結晶学のデータファイル(MTZファイル)など、バイナリファイルを直接読めると便利だ。R では、file 関数でファイルを開き、seek で目的の箇所に移動して、readBin でデータを読み込む。仕…

原子散乱因子 atomic scattering factor の異常分散項 f' と f'' は、CXRO X-Ray Interactions With Matter にある Atomic Scattering Factors から入手できる。Reference: B.L. Henke, E.M. Gullikson, and J.C. Davis. X-ray interactions: photoabsorptio…

R で ?dist のようにして使うヘルプ機能は大変便利だが、演算子について調べようとすると 演算子のヘルプ - ryamadaのコンピュータ・数学メモ(2018年9月に引越し) で記されているように不思議な挙動をする。例えば加算演算子+を調べるのに、? "+" や ? +3 は…

R でフーリエ変換を行うには、stat パッケージ(つまり、デフォルトで読み込まれる環境)に含まれる fft を使う。ついでに、いわゆる Fourier duck の実験をしてみた。 N <- 128 # 実空間の点の数 # トルエンっぽい平面分子。座標は手計算。原子距離は適当。 a…

逆格子の定義として、結晶学の教科書では外積を使った表現が出てくることが多いが、逆行列を使った表現のほうが「ちょっと実験」するにはコードしやすい。数学的な詳細はまたの機会に触れるとして、簡単に確認してみよう。以前 in-house の R-AXIS IV で取得…

R でデータフレーム等の各要素に関数を適用するには、apply の次元として c(1, 2) を指定する。apply は行ごと・列ごとに集計するイメージがあったが、第2引数としてベクトルを指定することで要素ごとに適用することも可能だったのだ。 vcf[, sample_names] …

R で plot を画像として保存する場合、 png("output%02d.png") のようにファイル名を指定すると、plot するたびに連番で保存してくれる。便利。ただし、カウンタは png コマンドを発行した時点で 1 に戻るので、既にあるファイルを上書きしないように注意が…

を台形則で数値積分してみて、分割数と誤差を評価する。 trapezoid <- function(fun, start, end, div) { x <- seq(start, end, length.out=div) step <- (end - start) / (div - 1) fx <- fun(x) return(step / 2 * (2 * sum(fx) - fx[1] - fx[div])) } div …

さっきのapply, lapply, sapply メモ - biochem_fan's noteに関連して、高階関数の型チェックつづき。グラフを書いてくれる curve 関数の引数は? fun <- function(x) {print(x); return(x * 2)} curve(fun, from=0, to=1) x の値をベクトルとして受け取って…

いつも、FUN の型がなにか分からなくなってしまうので実験してメモ。まず、行列に対して処理する場合。 m <- matrix(1:6, 2, 3) apply(m, 1, function(x) {print(x); return(1)}) apply(m, 2, function(x) {print(x); return(1)}) lapply(m, function(x) {pr…

R で整数の分割を列挙するには、partitions パッケージを用いる。 library(partitions) parts(5) # 5 の分割 (行列が返る。タテ方向に分割が入っている) P(5) # 上の個数 == ncol(parts(5)) restrictedparts(5, 3, include.zero=F) # 5を3つに分割 R(3, 5, i…

慣れの問題だと思うが、正直な話、Rのデータフレームの仕様はどうもしっくりこない。直感的に「こうだ!」と思って書くと、たいてい期待した動作にならないのでイライラする。そうやって苦労しつつ気づいた Tips をメモしておく。テーブル全体に散在する「デ…

id-val のペアがある。id は重複している場合もある。これを ID ごとに集計したい時は、tapply(val, id, sum) とすればよい。戻り値は array で、dimnames として id が入っている。問題は、dimnames は数値をセットしても文字列になってしまうこと。as.nume…

read.table で、データフォーマットは合っているのに、項目数が合わないとか文句を言われるときは、quote が悪さをしている可能性が高い。quote="" をつけるべし。

Ubuntu では、association study に用いる plink コマンドは putty-tools に含まれている plink コマンドと名前がぶつかるため、p-link に改名されていることに注意。