サイトトップ | Software Tips

統計処理ソフトウェアRについてのTips

最終更新:2016年 3月 17日 (木曜日)

このページでは,国際共同研究のオープンソースなプロジェクトで開発され,GNU GPLに従って公開,配布されている高機能な統計ソフトであるRについてのTipsを扱う。

TOPICS:
参照情報(リンク) | なぜR? | インストール | 基本操作 | 入出力 | グラフィック | 人口ピラミッド | 計算テクニック他 | 数学 | 続・数学 | 文字列

News/更新情報

保管庫1(2004年1月まで) | 保管庫2(2010年6月まで) | 保管庫3(2014年5月まで) | 保管庫4(2015年6月まで) | 保管庫5(2015年7月から)

保管庫内の主なトピック:
平方和(SS)「Rによる統計解析の基礎」(保管庫外だがサポート掲示板正誤表)|オッズ比レーダーチャート | EZR on Rcmdrの紹介文 | マルチレベル分析 | 栄養関連
1.6.01.6.11.6.21.7.01.7.11.8.0 | 1.8.11.9.01.9.1
2.0.02.0.12.1.02.1.12.2.02.2.12.3.02.3.12.4.02.4.1
2.5.02.5.12.6.02.6.12.6.22.7.02.7.12.7.22.8.02.8.12.9.02.9.12.9.22.10.02.10.12.11.02.11.1
2.12.02.12.12.12.22.13.02.13.12.14.02.14.12.14.22.15.02.15.12.15.2 | 2.15.3
3.0.0 | 3.0.1 | 3.0.2 | 3.0.3 | 3.1.0 | 3.1.1 | 3.1.2 | 3.1.3 | 3.2.0 | 3.2.1 | 3.2.2 | 3.2.3 | 3.2.4 | 3.2.4 Rev

3.2.4改訂版リリース(2016年3月17日)
R 3.2.4-revised is releasedというアナウンスが流れた。3.2.4は3.2系最後の超安定板であるべきなので、異例だが2つの不具合を修正したバージョンで置き換えたということらしい。
■Makefileの不具合1つと,format.POSIXlt()のバグフィックスとのこと。
■例えば,format(as.POSIXlt(paste0(1940:2000,"-01-01"), tz="CET"), usetz=TRUE)と打つと,最後の2つが"CET"でなく"CEST"という表示になっていたが,それを修正した。
mapdataパッケージのmap()で太平洋を描く工夫(2016年3月16日)
■mapdataパッケージでmap("worldHires")の地図は日付変更線が左端のX座標-180(西経180度)と右端のX座標180(東経180度)になっているので,太平洋を中心にサモアとハワイを1枚の図に入れようと思うと工夫が必要。
■layout(t(c(1,1,1,2,2)))としてグラフィックウィンドウを左対右が3:2になる二画面に分割し,まずpar(mar=c(1,1,1,0))で左画面の右マージンをゼロにしてからmap()のオプションでxlim=c(120,180)と指定して60度分を描画し,次にpar(mar=c(1,0,1,1))で右画面の左マージンをゼロにしてからmap()のオプションでxlim=c(-180,-140)と指定して40度分を描画すると,それらしい地図になった。普通はどうするんだろうか?
順位相関係数の信頼区間についてのメモとR-3.2.4インストール(2016年3月14日,鵯記,15日に追記)
■奥村さんが昨日tweetしていたKryvasheyeu Y et al. (2016) Rapid assessment of disaster damage using social media activity. Science Advances, 2(3), e1500779, DOI: 10.1126/sciadv.1500779は,確かにρやτのような順位相関係数とともに,明らかに相関があるのにp値を出しているのはあまりinformativeではない(この場合,検定の帰無仮説は「無相関」だから,p値がどれほど小さかろうと,「無相関ではない」確信の度合いが増すだけ)。当然,信頼区間を出す方が良いのだが,Rのcor.test()ではmethod="pearson"の場合しか信頼区間は表示されない。cor.test()でmethod="spearman"のとき,二つの変数をxとyとするとスピアマンの順位相関係数ρそのものはcor(rank(x), rank(y))で求めている(getS3method("cor.test", "default")と打って内部コードを見るとそうしている)が,ρをz変換したものが近似的に正規分布に従うと仮定しないとcor.test(rank(x), rank(y))で出てくる信頼区間がρの信頼区間にはならないので,そう言い切ってしまうのは抵抗を感じる。たぶんbootstrapで求めるのが原理的には正しいやり方で,この記事によると,ρの信頼区間を出してくれるのはパッケージRVAideMemoirespearman.ci()関数であり(bootstrapの回数はnrep=オプションで指定するが,まあ1000もあれば良かろう),ケンドールのτの信頼区間を出してくれるのはNSM3パッケージのkendall.ci()関数である(bootstrapするときは,bootstrap=TRUEオプションを付け,B=オプションで回数を指定する)。せっかくだからエビデンスベーストヘルスケア特講のテキストの順位相関係数のところに付記しておこう。なお,当該論文のFig.2のパネルAみたいなU字型の関係は順位相関ではダメなので,かつてメモしたMICを使うしかないだろう。
■(2016年3月15日追記)上記のいくつかの方法で順位相関係数の信頼区間を求めるためのRコード
■R-3.2.4がミラーにも来たのでWindowsにもインストール中。mran (Microsoft R Application Network)から,Microsoftが買収するまでRevolution R Openだったものが(つまり,マルチコア対応のRが),Microsoft R Openとしてダウンロードできるのだが,そっちはまだ3.2.3だな。
R-3.2.4リリース(2016年3月10日,鵯記
R-3.2.4が予定通りリリースされた。よく使う機能で目に見える変化は,cor.test()がきわめて小さいp値をより正確に計算するようになったという点くらいか(試しに,まだ3.2.3のWindows版でcor.test(1:10, c(2:10, 11.2))と実行するとp-value = 2.665e-15となったが,同じコードを3.2.4にアップデートされたubuntu上のRで実行したらp-value = 2.54e-15となった。95%信頼区間はどちらでも同じ値が表示されていた)。今回,[Rd]メーリングリストに,リリース直後にOracle Enterprise Linuxでビルドしようとした人から,たぶん@であるべきところが$になっているためにエラーが出てビルドできないというバグレポートがあり,すぐに,R-patchedで修正済みでR-3.3.xではなくなる問題だという返事があった。もちろん$を@に直すだけでもビルドできるとのこと。
cranパッケージ8,000到達(2016年2月29日,鵯記
■[Rd]メーリングリストにcranパッケージ数が8,000に到達という投稿があったので推移をグラフにしてみた(tweet)。コードは140字に入らないので,スラドに投稿しておいたが,ここにも貼っておく。

Npkg <- data.frame(
dates = as.Date(c("2003-04-01","2004-10-01","2007-04-12","2009-10-04",
"2011-05-12","2012-08-23","2013-11-08","2014-10-29","2015-08-12",
"2016-02-29")),
N = c(0.25, 0.5, 1:8)*1000)
plot(Npkg, type="b",main="CRANパッケージ数の推移")
# Source: https://stat.ethz.ch/pipermail/r-devel/2016-February/072388.html

■このコードで描かれるグラフを見ると,cranに登録されるパッケージ数が加速度的に増えていることがわかる。cran以外にもGitHubだけにあるパッケージとかもたくさんあるので,今後どこまで増えていくのか,まだ見当もつかない。
R-3.2.4とR-3.3.0のリリース予定について,及び多項ロジスティック回帰(2016年2月23日,鵯記
■2月17日にアナウンスメールが流れていたのをメモし忘れていたが,Rの次のリリースは3.2系最終版のR-3.2.4(コード名:"Very Secure Dishes")で3月10日予定。その次はいよいよ3.3系に入ってR-3.3.0が4月14日リリース予定なわけだが,公表されたコードネーム"Supposedly Educational"をこれまで通りPeanutsから探しても,これくらいしか近いものが見つからなかった。しかしこれはSally Brownの言葉だが,コードネームとは若干違って"Supposed to be educational"となっている。Peanuts以外からそのものずばりを探すと,文の一部として使われている例(とくに「おそらく教育的な」と"Slave Tetris"を形容するもの)が多かったが,イカレたプログラマの短編コミックみたいなものがあった。なかなか面白いので他のも見てみたいところだが,ひょっとするとこの方のオリジナルなんだろうか?
■多項ロジスティック回帰分析をしたいという相談を受け(厳密に言うと,アウトカムが3水準のカテゴリで,3水準以上のカテゴリからなる要因について,リファレンスカテゴリに対しての各アウトカムのオッズ比を,他の変数の影響を調整した上で得たいという相談だったが,たぶん素直に考えたらこれだろう),共立出版の「Rで学ぶデータサイエンス」シリーズに入っている,藤井良宣『カテゴリカルデータ解析』を貸して,第6章を読んで真似してみたら? とアドヴァイスしたのだが,6.5にあったのはnnetパッケージのmultinom()関数を使うやり方で,オッズ比を求めるのが若干大変かもしれない。他のパッケージと比べたページを見たところ,VGAMパッケージのvglm()関数を使うのがわかりやすいと思った。とくにこのページは参考になると思う(翌日,エビデンスベーストヘルスケア特講Iのテキストに追記した)。
RjpWikiの旧URL廃止(2016年2月16日,鵯記
■RjpWikiの旧URLが廃止になったということでリンクを修正していて気づいたが,R-jpメーリングリストが消滅してしまったようだ。こうして時代は移ってゆく。
『R言語徹底解説』の書評情報(2016年2月12日,鵯記
■Qiitaに「R言語徹底解説」入門と題した同書のレビューがあることを@kazutanさんのtweetで知った。実に的確なレビューだと思った。
世界の年平均気温の推移について(2016年2月2日,鵯記
気象庁が発表した世界の年平均気温の推移について,Rでweb上の表からデータをスクレイピングして図示するコードはXMLパッケージを使えばとても簡単で,数行で済む。たぶん数値が更新されてもコードは書き換えずに新しいグラフが作れるはず(なお,グラフの縦軸は基準値からの偏差なので注意。実際の値には意味が無いそうだ)。
世界の年平均気温の推移
世界全体についての回帰直線を追記したいなら,abline(lm(x[,2]~as.integer(strtrim(x[,1], 4))), lty=1, col=1) でいいし,x[,2]x[,3]にすれば北半球,x[,4]に変えれば南半球の回帰直線が描ける(重ね描きしていくなら,もちろんltyとcolもそれぞれ2と3にする)。もう少し真面目に時系列データにするならこうか? こうしておけば,library(forecast); plot(forecast(auto.arima(z[,1]), h=50))で世界全体の今後50年の幅付き予測とかもできる。
相対値と順位の関連を見るための散布図を描く方法(2016年2月2日,鵯記
■十数種類のモノについて数十人の対象者にある軸での順位と,別の軸での0~100%の相対値をつけて貰い,それぞれのモノについて順位と相対値の関係を見るための図を描こうとすると,どうしても散布図では重なりが出てしまって様子がわからない。さりとて普通にstripchart()で描こうと思うと,無い順位が詰まってしまって,ほぼ全員が高い順位を付けたモノと,ほぼ全員が低い順位を付けたモノの違いが図の上で一目でわからない。つまり別々のモノについて,共通の順位軸を与えたい。こういう場合は,順位をそのまま数値としてstripchart()に与えるのではなく,factor(levels=)を使って要因にしてしまえばよい。例えば,

x1 <- sample(1:6, 100, replace=TRUE)
x2 <- sample(10:15, 100, replace=TRUE)
y <- sample(1:10, 100, replace=TRUE)
layout(1:2)
stripchart(y ~ factor(x1, levels=1:15), vert=TRUE, method="jitter")
stripchart(y ~ factor(x2, levels=1:15), vert=TRUE, method="jitter")

とすれば良い。

Rの参照情報

なぜRを使うべきなのか?

Rの最大の利点は,オープンソースなフリーソフトで,かつ拡張性が高い点である。世界中の研究者がGISを含む空間統計解析やゲノム解析などに至るまでさまざまな追加ライブラリを公開しているしCRAN(http://cran.r-project.org)というサイト(筑波大ミラー)に集積される仕組みもあるので,CRAN内で検索すれば大抵の処理は見つけることができる),自分で新しい拡張関数を付け加えることもできる(群馬大学社会情報学部の青木繁伸教授のように,検索するより作ってしまう方が早いと言われる方も多いが)。

オープンソースということは,誰でもその気になれば計算の中身をインプリメンテーションレベルでチェックできることを意味する。これは,商用ソフトにはありえない利点である。商用ソフトでは,利用している計算式はわかっても,コードそのものは公開されないために,インプリメンテーションにバグがないかどうかは,サンプルデータ解析を実行してクロスチェックすることでしか確認できない場合が多い。Rの場合は,世界中の人が使いながらコードチェックもしているので,計算の信頼性もかなり高い。

また,R自体の開発は,メーリングリストで連絡をとりながら行われているが,コアグループが受け入れないとバージョンアップに取り入れられない。コアグループは,かなりはっきりした思想をもってコード開発を行っているので,コードが入り乱れてしまうことはない(これは,裏を返せば,コードをチェックしておかしいと思われる点を報告しても,開発コアグループを説得できないと「仕様だ」ということになって,バージョンアップに取り入れられない場合もあることを意味する)。

プログラムコードが洗練されているため,ソフトウェアのサイズが小さく,動作が軽快なのも利点である。大学等での統計教育を考えると,完全に無料で利用できるため,卒業後も覚えた技術が無駄にならない(SASやSPSSのような高価なソフトウェアで実習を受けた場合,卒業しても使えるような環境に就職できる人はほんの一握りなので,多くの学生はその技術を活用することができなくなってしまう)。教育面では,多くの有用なサンプルデータが含まれているため,いちいちデータ入力をしなくても,分析手法に合わせて適切なデータを利用できる点も利点といえる。Fisherのirisデータのような有名なものは当然として,白血病患者の生存時間データとしてよく使われるGehanのデータも,MASSライブラリ(RecommendedなライブラリなのでWindows版バイナリディストリビューションには標準で含まれている)にgehanとして入っている。

また,国際協力などの場面でもライセンスを気にすることなく共用することができる。英文のみならず,仏文,西文などのマニュアルも公開されている。Windowsだけでなく,MacintoshでもLinuxでもFreeBSDでも動作するので,さまざまな環境で同じ統計解析を行なうことができる。R以上に各国語対応している統計解析のフリーソフトウェアとして,米国CDCが提供しているEPIINFOがあるが(ただしEPIINFOはWindowsのみ),利用できる統計解析手法の種類はRの方がずっと多い。RにはSPSSでさえ実装されていないような新しい分析手法が多く含まれている。結果の信頼性も高く,最近では多くの学術論文が統計解析にRを使っている(例えば,2004年2月にはNatureにもRを使って分析された論文が掲載されている。Morris RJ, Lewis OT, Godfray HCJ: Experimental evidence for apparent competition in a tropical forest food web. Nature 428: 310-313, 2004.)。

もちろん,使いやすさもかなり高い水準にある。例えば,Excelで独立2標本の平均値の差の検定をするには,2群の標本データをワークシート入力し,メニューのツールから分析ツールを選んで(パッケージのインストール時にアドインとして分析ツールをインストールしておかないとメニューに出てこない),等分散でないときの2群の平均値の差の検定を選んで2つの標本の範囲を選んで実行しなくてはならない。結果は別々のシートに出力されるが,余分な統計量が雑然と並んでいて,表の体をなしていない。

Rで同じことをする場合,サンプルサイズが小さければ,変数xy(変数名は何でも構わない)に直接2つの標本データを付値(代入)してから,t.test(x,y)一発でWelchの検定が完了する。

独立2標本の平均値の差の検定は,古典的には,まずF検定をして2群の分散に差がないかを調べ,差がない場合は通常のt検定,有意差があったらWelchの検定と2段階でなされてきたが,奥村先生青木先生のシミュレーション結果を見ると常にWelchの方法が最良の結果が得られることが明白なので,この辺りの記述を改めた。なお,古典的な2段階式の場合,Excelではいくつもの手順が必要だが,Rでは関数として

a.t.test <- function(x,y) {
if(var.test(x,y)$p.value < 0.05) {t.test(x,y)} else {t.test(x,y,var.equal=T)}}

と定義しておけば,次からはa.t.test(x,y)とするだけで,自動的に等分散性の結果で条件判定して適切な検定をさせることができる。もっといえば,

t.test(x, y, var.equal=(var.test(x,y)$p.value>=0.05))

とすれば,関数定義さえ必要とせず2段階式の検定ができる(頭の固い雑誌のレビューアがどうしてもやれと言ってこない限り,やる必要はないが)。もちろん,表形式のデータを読み込んでから,分析する関数だけ指定することもできる。

ヘビーユーザーにとっての利点は,RがS言語にほぼ互換な言語のインタープリタであり,それが関数型言語だということから生まれる。つまり,結果を変数に代入して保存したり加工したりできる。xtableというライブラリをインストールしておけば,結果をxtable( )の括弧内に入れるだけで,HTML形式やLaTeX形式に変換でき,きれいな表が得られる。

さらに素晴らしいことに,そうやって実行したすべてのプロセスを,テキストファイルとして記録し,保存しておけるので,後になって,どういう分析をしたかをチェックすることができる。しかも,保存しておいたファイルは(例えばtest.Rというファイル名だとすると),source("test.R")とすると再実行できる。どんなに複雑な作業をしても,それを何度でも簡単に再現できるということである。逆に考えれば,適当なテキストエディタでプログラムとしてRのコマンドを書き連ねておいたものを読み込ませれば,複雑な分析手続きでも1回の操作で終わらせることができる。R Commander(Rcmdr)というライブラリを使うと,メニュー形式で操作することもある程度可能になるが,その場合でも,ログはきちんと関数リストとして保存される。

美しい図を作るのも実に簡単で,しかもその図をpdfとかpostscriptとかpngとかjpegとかWindows拡張メタファイルの形式(emf)で保存でき(ただしemfで保存できるのはWindows版のRのみ),他のソフトに容易に取り込める。例えばwin.metafile()関数を使ってemf形式で保存すれば,Microsoft PowerPointやOpenOffice.orgのDrawなどの中で,ベクトルグラフィックスとして再編集することが可能である。

以前は,多くの日本人にとって最大の難点は,日本語が使えないことだったが(データとしては入っていても大丈夫だったが変数名に使えなかったしグラフ内で使えなかった),中間栄治さんが日本語も扱えるようにするパッチを開発して公開されたのでこの問題は解決した。バージョン2からは本体が国際化対応したので,日本のRユーザ有志の手によってメッセージまで日本語化されたものも使えるようになった。2010年10月16日現在の最新版は2.12.0である。

かつては日本語による解説書があまりなかったが,2003年10月に出版された拙著『Rによる統計解析の基礎』(ピアソン・エデュケーション)を皮切りに,現在では多数出版されているので問題ない。ウェブ上の情報も,群馬大学社会情報学部の青木繁伸教授のサイト内の「Rによる統計処理」や岡田昌史氏による「RjpWiki」を始めとして充実したので,環境としてはRを使う上で支障はなくなったといえよう。

インストール関連(2016年3月16日更新)

@ITで実践!Rで学ぶ統計解析の基礎を2010年から2011年まで連載されていた柏野雄太さんが自社サイトで提供されていたサポートページの動画は大変わかりやすかったが,いつの間にか消滅してしまったようだ。2016年現在の動画としては,

Windows
山形大学統計数理研究所のCRANミラーサイトからR-3.2.4のインストール用バイナリファイル(R-3.2.4-win.exe)をダウンロードし,ダブルクリックして実行する。
Japaneseを選んでリターンキーを押すと,[セットアップ - R for Windows 3.2.4]というウィンドウが起動するので,[次へ]というボタンをクリックし,ライセンス表示にも[次へ]をクリックし,インストール先は特別な理由がなければデフォルトのままで[次へ]をクリックする。次の,インストールするコンポーネントを選ぶウィンドウもそのまま[次へ]ボタンをクリックすると,スタートアップオプションをカスタマイズするか尋ねるウィンドウが表示されるので,ここはデフォルトの[いいえ]でなく,[はい]の方をマークして[次へ]をクリックすることをお薦めする(スタートアップオプションがデフォルトでは,Rを起動した後のすべてのウィンドウが,1つの大きなウィンドウの中に表示されるMDIモードになってしまうのだが,それだとR Commanderが非常に使いにくくなるからである)。次に表示されるウィンドウでSDIにチェックを入れて[次へ]をクリックする。次のHelpの形式はどちらを選んでも良いが,筆者はテキスト形式の方が好みである。[次へ]をクリックすると,インターネット接続を標準設定にするか,Internet Explorerのプロキシ設定に合わせる(Internet2)かを聞いてくるが,通常はどちらでも問題ないはずである(プロキシを介している場合はInternet2を選択)。[次へ]をクリックするとスタートメニューフォルダ名を聞いてくるが,ここはデフォルトのままで問題ない。次のウィンドウではデスクトップアイコンを作り,Rのバージョン番号をレジストリに記録し,.RDataという拡張子をもつファイルをRに関連付けするというオプションがデフォルト指定されているが,通常はそのままで問題ない。
[次へ]をクリックするとインストールが始まる。暫く待つとセットアップウィザードが完了したという意味のウィンドウが表示されるので,[終了]をクリックする。以上の操作でR本体のインストールは終わりである。
Macintosh
群馬大学社会情報学部・青木繁伸教授のサイトに詳細な解説記事があるので参照されたい。
Linux
メジャーなディストリビューションについては有志がコンパイルしたバイナリがCRANにアップロードされているので,それを利用すればインストールは容易であろう。ubuntuの場合なら,Ubuntuソフトウェアセンターの「科学&工学」の「数学」からRを選んで「インストール」ボタンをクリックするだけでインストールできる。ただし,rglパッケージなど,X11を要するパッケージを使う場合は,端末からsudo apt-get install r-cran-rglもやっておく必要がある。また,gmpパッケージを使う場合は,端末からsudo apt-get install libgmp3-devもやっておく必要がある。マイナーな環境の場合や,高速な数値演算ライブラリを使うなど自分のマシンに最適化したビルドをしたい場合は,CRANミラーからソースR-3.2.4.tar.gzをダウンロードして展開して自力でコンパイルする。最新の環境であれば,./configuremakeでビルドしてから,スーパーユーザになってmake installで済むことが多いが,場合によっては多少のパッチを当てる必要があるかもしれない。
追加パッケージ
追加パッケージのインストールは,ネットワークにつながった環境なら,Rを起動して,install.packages("パッケージ名", dep=TRUE)とすればいい(CRAN=というオプションを指定しない限り,1度のRの起動ごとに,どのミラーからダウンロードするか尋ねるウィンドウが表示されるので,そこから選ぶ)。UNIX系OSではRのディレクトリへの書き込み権限のあるユーザで実行しないと失敗する。Windows環境では,Windows版バイナリがあるパッケージなら,CRANには会津大学のミラーサーバを指定してもいい。例えばWindows版のRにsemパッケージをインストールする場合は,install.packages("sem", dep=TRUE)とすればよい。パッケージ名はCRANのパッケージソースのページで調べられる。また,Windowsでハードディスクなどに保存しておいたzipファイルからインストールするには,例えばinstall.packages("./survival.zip", .libPaths()[1], CRAN = NULL)のようにすればよい。
■既にインストール済のパッケージのバージョンなどを確認するには,packageDescription("パッケージ名")またはlibrary(help=パッケージ名)とする。前者は二重引用符が必要で,後者は二重引用符を付けない。
Windows環境の場合のパッケージ管理のヒント(2013年1月16日修正;2013年5月20日追記;2016年3月16日に大幅書き換え)
●Windows 7以降の場合,追加パッケージはとくに指定しなければシステムディレクトリでなくユーザごとのディレクトリ(例えば,C:\Users\Minato\Documents\R\win-library\3.2)にインストールされる。R本体がメジャーバージョンアップされると,パッケージのインストール先ディレクトリ名が変わるので,新しいバージョンのR(例えばR-3.3.0)をインストールした後で,新しいバージョンのRのライブラリフォルダ(例えばC:\Users\Minato\Documents\R\win-library\3.3)にそれらを全部移動(ただし既に同名のファイルやフォルダがある場合は,更新されたバージョンがシステムと一緒にインストール済みということを意味するので移動しない)してから,新しいバージョンのRを起動し,プロンプトにupdate.packages(checkBuilt=TRUE, ask=FALSE)と入力する。
●この方式は,とくにネット回線が遅い場合に時間の節約にはなる。しかし,バージョンアップを繰り返すと,古くなってdeprecateされたパッケージや使わなくなったパッケージが残り続け,だんだん重くなってくる欠点がある。使うパッケージがある程度決まっているなら,ホワイトリスト方式の方が良いと思う。つまり,常用しているパッケージを予めリストしておき,下記リンク先のように,Rのコードとしてinstall.packages()を並べておいて,Rの新バージョンをインストールした後はいつも,このコードを実行するようにすれば良い。というわけで,中澤はこの方式をお薦めしたい。

source("http://minato.sip21c.org/swtips/instmyusuallibs.R", encoding="UTF-8")

Windows環境における起動アイコンのカスタマイズ(2013年1月16日修正)
●Windows環境では,Rの作業ディレクトリは,Rを起動するショートカット("%Program Files%R\R-3.2.4\bin\x64\Rgui.exe"[64ビット版の場合]または"%Program Files%R\R-3.2.4\bin\i386\Rgui.exe"[32ビット版の場合]を直接ダブルクリックしているのでない限り,Rはショートカットから呼び出されているはず)の「プロパティ」の「作業フォルダ」に指定したものになる(そうすると,Rを起動してからFileのChange Dirで作業フォルダを指定するのと同等のことが起動時からできる)。
●起動アイコンを複製し,片方を"R x64 3.2.4"(ショートカットなので拡張子は.lnkである),もう片方を"R x64 3.2.4-en"などと名前をつけて,"R x64 3.2.4"のアイコンで右クリックして「プロパティ」を選び,「リンク先(T)」の右側のテキストボックスに"C:\Program Files\R\R-3.2.4\bin\x64\Rgui.exe"となっているので,その後に半角スペースを1つ空けてから--sdi LANG="ja"と文字列を追加してOKする。次に"R x64 3.2.4-en"のアイコンで右クリックし「プロパティ」の「リンク先(T)」の方では,同じく"C:\Program Files\R\R-3.2.4\bin\x64\Rgui.exe"となっているので,その後に半角スペースを1つ空けてから--sdi LANG="en"と打ってOKする。すると,"R x64 3.2.4"アイコンからRを起動すれば日本語メニュー・日本語メッセージでRが起動し,"R x64 3.2.4-en"アイコンからRを起動すれば英語メニュー・英語メッセージでRが起動する環境が構築できる。Rcmdrパッケージは多言語対応なので,日本語メッセージのRから呼び出せば日本語メッセージになるし,英語メッセージのRから呼び出せば英語メッセージになる。これにより,Windows自体が日本語版であってもRは英語メッセージで使うという環境が実現できる。RcmdrPluginとしてEZRを使うときも言語環境は引き継がれる。

最も基本的な操作(2010年1月5日追補)

データの入出力関連

グラフィックテクニック

グラフィック上での日本語を正しく扱うための注意事項(2010年5月19日再掲)
●グラフィックデバイスによって異なる。
postscript()やpdf()では,中間さんのAI_UCS2.Rをsource()を使うなどして先に実行してからグラフィックデバイスを開き,par(family="Japan1GothicBBB")をしてグラフ出力すべき。はしご高のように,UTF-8では表示できるがEUC-JPでは表示できないような文字も表示できるようになる(参考:2008年6月23日のメモにいただいたコメント)。
win.metafile()では
windowsFonts(JP1=windowsFont("MS Gothic"),JP2=windowsFont("MS Mincho"))
par(family="JP1")
を実行してからplot()などのグラフ出力を実行すべき。こうしておけば,拡張メタファイル(emf)をPowerPointやOpenOffice.orgのImpressなどに読み込んでから編集するためにオブジェクト変換しても漢字が文字化けしない。
●いずれの場合でも,描画が終わったらdev.off()を忘れずに。
日本語をベクトルグラフィックとしてプロット(注:歴史的意味しかない技)
●RにはHersheyEUCとして約700の日本語フォントが含まれていて,text関数から使えるのだが,フォントの指定が面倒なので使いにくかった。これを簡単に使えるようにする関数jtext(実行例を含むソース)を書いてみた。ただし内部で漢字文字列をJISコード番号文字列に変換するPerlプログラムjiscode.plを呼び出していて,その中でcgiなどで良く使われているjcode.pl(リンク先は一次配布元FTPサイト内のver 2.13のファイル)を呼び出しているので,それらが実行時に呼び出せるディレクトリに置かれていて,かつActive Perlにパスが通っている必要がある。もっとも,LinuxやFreeBSD環境ならばjiscode.plの先頭にPerlのパスを書けばいいはずである。なお,Perl5.8以降では標準のEncodeライブラリを使う方がスマートな気がするし,それ以上に,単にjisx0208にできればいいのだからCで書いてライブラリ化した方が便利だと思うが,とりあえずPerlで用が足りるので今後の課題ということにする。
実行例を下に示す(上でリンクしているRプログラムを実行してできるpngファイル)。
日本語文字列プロットの例
●なお,pLaTeXと組み合わせて使えば,上のような面倒なことをしなくても,日本語を含むプロットが簡単にできる(ベクトルグラフィックとしてではなく,文字としてのプロットになる)。textで普通に日本語文字列をプロットした後で,dev.copy(pictex,"jtext3.tex"); dev.off()としてpictex形式のファイルjtext3.texを作り,これを
\documentclass{jarticle}
\usepackage{pictex}
\usepackage[dvipdfm]{graphicx}
\begin{document}
\begin{figure}[h]
  \centerline{\input{jtext3.tex}}
  \caption{3都市の緯度と経度の2次元表示}
\end{figure}
\end{document}
のようにして読み込むファイルtest.texを作って,platex testdvipdfm testとすれば,test.pdf (16,518 bytes)ができあがる。こちらの方が使えるフォントが多いし,美しくて実用的である(時折そのままでは通らない文字があるのだが,出力されるpictexファイルを直接エディタで編集すれば問題ない。textだけではなく,xlabとかmainでも日本語が使えた)。
群ごとにマークを変えたプロット
●簡単にやるには,plot関数のpch=に値を与えればよい。例えば,X,Y,Zという3つの村があって,それぞれ身長と体重のデータがあって,その関係を村ごとにマークを変えてプロットしたいとする。
●データは1行目がタブ区切り変数名で2行目以降がタブ区切りデータのテキストファイルからread.delim関数で読む(上述)のが普通だが,プログラム内で与えるには,
hx <- c(161.5, 167.0, 157.5, 160.5, 156.0, 165.0, 157.5, 169.5, 168.5, 178.5); wx <- c(49.2, 72.8, 61.6, 62.6, 73.6, 57.8, 60.8, 92.4, 58.2, 71.8); hy <- c(159.5, 161.5, 148.0, 164.0, 169.5, 164.5, 161.0, 172.0, 168.0, 163.0, 153.5, 165.5, 168.0, 161.0, 160.5, 158.0, 156.0, 155.5); wy <- c(43.6, 65.2, 49.6, 57.4, 67.0, 61.4, 56.2, 66.4, 67.0, 60.2, 44.6, 63.2, 62.4, 53.8, 63.6, 59.0, 59.0, 53.6); hz <- c(171.5, 176.5, 166.0, 175.0, 166.0, 169.0, 175.5, 167.0, 163.5); wz <- c(73.2, 77.8, 58.0, 57.2, 67.0, 73.6, 97.8, 63.0, 67.4);として村ごとのデータを与え,
x <- data.frame( VG=c(rep('X',NROW(hx)),rep('Y',NROW(hy)),rep('Z',NROW(hz))), HEIGHT=c(hx,hy,hz), WEIGHT=c(wx,wy,wz) )としてデータフレームにすればよい。
●村名をそれぞれ違う色で身長と体重の座標位置にプロットするには,
plot(x$WEIGHT ~ x$HEIGHT,pch=as.character(x$VG),col=as.single(x$VG),main="Relationship between HEIGHT(cm) and WEIGHT(kg)\n in Solomon Adult Males",xlab="HEIGHT(cm)",ylab="WEIGHT(kg)")とすればよい。下図が得られる。なお,村ごとに細かく指定したい場合は,plot(x$WEIGHT[x$VG=='X'] ~ x$HEIGHT[x$VG=='X'], pch='X', col="blue", xlim=c(140,180), ylim=c(40,100))などとしてから,points関数を使って重ね書きすればよい。なお,この例のようにプログラム内で付値したのであれば,もちろん,plot(hx, wx, pch='X', col="blue", xlim=c(140,180), ylim=c(40,100)); points(hy, wy, pch='Y', col="red"); points(hz, wz, pch='Z', col="green")でよい。
Scatter plot by village as a sample of R techniques
●もっとも,3種類くらいならいいが,120種類のマークを変えるのは難しい。マークとして使えるのは1文字だけなので,アルファベットや数字では不足してしまう。解決策は3つあって,
  1. 色とマークを組み合わせる: pchに与える値として,1から25まではプロットとして適切なマークが既に定義されている(26から32は空白で,33以上は文字や記号)ので,col="red"とかcol="blue"とか指定して(または剰余を使うなどして周期的変数を生成して)色を変えれば120くらいは何とかなりそうである。欠点は,モノクロプリンタに印字すると色の区別がつかないことである。
  2. text( )関数を使って文字列を重ね打ちする: plot(x,y)の後に(たぶんpch='.'とかpch=20とかで小さい点にしておいた方がいいと思う),text(x, y, paste(string), pos=4, offset=0.5)などとすれば,stringを(x,y)の点の右側に印字してくれる。(plotでlog="x"により横軸を対数軸にしておいても,text( )による文字追加は自動追従してくれる)。この場合文字を小さくしておきたいので,予めpar(ps=6)などとしておく。それでも重なる場合はあると思う。(注:「日本語をベクトルグラフィックとしてプロット」の項に書いたように,文字列をベクトルグラフィックとしてプロットすることもできる)
  3. identify( )関数を使う: すべてのデータ点を特定する必要はないので,必要な点についてだけ情報を表示できるのがベストであろう。plot(x,y)の後にidentify(x,y,labels=string)としておくと,プロットの後に十字型のマウスカーソルが出現するので,画面のstringを表示したい点の上でクリックすればstringが出現する。描画ウィンドウのメニューのstopからか,右クリックメニューからstopするまで,複数の点をクリックできる。
度数ゼロを含む度数分布図を描く
●連続した自然数をカテゴリとする度数分布を描きたいけれども,度数がゼロのカテゴリを飛ばしたくない場合がある。例えば,子ども数の分布を描きたいときは,たまたま子ども数7人の母親がいなかったとしても,8人の人がいたら,子ども数7人の母親の度数がゼロであることも,図には表したいだろう。とくに,複数の図を比較する場合は尚更である。ところが,普通に子ども数を代入(付値)した変数をxとすると,barplot(table(x))では,度数がゼロのカテゴリがtableの出力に出てこないので描画されないのである。この解決のために汎用の関数を作ってみたので紹介する。
●関数定義は以下の通り。

sbarplot <- function(freq,xnames,xlabel,title) {
maxc <- max(as.single(dimnames(freq)[[1]]))
z <- rep(0,NROW(xnames))
for (i in c(1:NROW(freq))) {
z[as.single(dimnames(freq)[[1]][i])+1] <- freq[[i]]/sum(freq) }
barplot(z,names.arg=xnames,main=title,xlab=xlabel)}

●使い方は,上の関数をRに読み込んだ後で,例えば,y <- c(3,4,2,1,3,4,2,2,3,2,1,2,6,6,7,4,2)というデータがあったときに,0や5も含めた度数分布図を書きたければ,x <- c(0:7)として横軸を定義しておいて,sbarplot(table(y),x,"label for x axis","title of the figure")とするだけで,下図が得られる(なお,解像度を指定してpngファイルを作るには,画面にグラフが描かれている状態で,例えば,dev.copy(png,"sample.png",width=400,height=480); dev.off()とすれば,画面に出ているグラフがsample.pngというファイルに保存される)。縦軸を割合でなくて実数にしたい場合は,sbarplot関数の中から/sum(freq)を削除すればよい。
frequency bar plot including zero category
●forなど使っていてR的ではないコードだし,おそらくもっと簡単な方法があると思うので,改善のご意見がいただければ幸いである。
●群馬大学の青木先生からいただいた2通りの関数(実行の仕方を含む)を掲載する。例えばx<-c(162, 159, 163, 157, 152, 168, 153, 156, 167, 161, 154, 162, 160, 157, 169, 160, 162, 158, 161, 160, 163, 160, 163, 153, 164, 163, 163, 153, 155, 155, 162, 163, 168, 160, 158, 168, 163, 163, 158, 153, 161, 153, 168, 156, 155, 159, 158, 161, 157, 155, 161, 156, 167, 156, 158, 152, 160, 160, 155, 157, 158, 160, 157, 156, 164, 157, 161, 158, 161, 153, 163, 161, 160, 162, 159, 162, 161, 158, 160, 177)としてから試されたい。
  1. mbarplot.R。使い方はxをデータベクトルとして,mbarplot(x,5)とすれば階級幅5の度数分布が計算され,同時にグラフが描かれる。
  2. dosuubunpu.R。使い方は同じくxをデータベクトルとして,dosuu.bunpu(x,2.5)とすれば階級幅2.5の度数分布が計算され,同時にグラフが描かれる。
◎この項,長らく放置していたが,裏RjpWikiさんの記事を見て目からウロコが落ちた。なるほど,tableをとる前にfactor化してlevelsを強制指定すれば,それだけで良かったのだな。つまり,上述のyについて0と5も含んだ度数分布図を描くためには,barplot(table(factor(y, levels=0:7)))とすれば良かった。

人口ピラミッドを作る(2014年9月4日 ver 1.4)

■pyramid()関数の開発
●かつては,軸を書かせないでbarplotしてからaxis関数でカスタマイズした軸を出力するなど,いろいろなテクニックを使って作っていたが,関数内でpar()を操作して画面を2分割して描画するなど不自然であり,間に合わせとしか言いようのない関数であった。当時の関数はpyramid.Rとして残しておく。当時の描画例も下に残しておく。
Population pyramid -- old version
◆なお,人口学のページで紹介している人口ピラミッドの作り方は,まだこの旧バージョンによる説明が残っているが,下記の新しいバージョンでも試してみたので参照されたい。年齢階級数が多い場合は,まだ微調整が必要かもしれない。
●その後数年が経つうちに,低レベル関数の扱いに慣れてきたし,レーダーチャートを作る関数なども作ってみたりしたので,自分で長方形の座標を計算して描画するように全面的に書きかえることにした。ちまちまと計算式を書いただけなので,何も難しいことはしていないし,じっくりソースを読めば誰でもわかるコードだと思う。ポイントは,最初に何も描画せずに左下が(-1.1-Cgap/2,-0.1),右上が(1.1+Cgap/2,1.1)となる架空座標系を設定した上で,ラベルや軸要素は最初から架空座標系で描き,データは架空座標系に変換してプロットさせている点である。たぶんまだ改良の余地があると思うが,とりあえず以前のバージョンより使いものになると思う。
●圧縮したファイルは,pyramid_1.1.zip(rev 1.1; 13,031 bytes),またはpyramid_1.1.tar.gz (rev 1.1; 3,513 bytes)なので,ダウンロードしてからパッケージインストールして,library(pyramid)とすれば,pyramid()とpyramids()が使えるようになる。example(pyramids)で描ける図を下に示す。
Population pyramid -- new version
■パッケージ化
●Windows XP環境での作業を書いておくが,LinuxやFreeBSDでも大差ないであろう。なお,ここではRの標準関数を組み合わせて作った関数について説明するが,srcサブディレクトリとinstサブディレクトリを作ってCやC++やFORTRANで書いた外部プログラムを組み合わせることも可能である。
●まず,パッケージ用に適当なディレクトリを用意する。仮に,"d:\work\pyramid"としよう。ここに,テキストファイルDESCRIPTIONを用意する。DESCRIPTIONファイルは,例えば次のような内容にする(継続行は行頭を空白文字またはタブコードにすることに注意)。
Package: pyramid
Version: 1.1
Date: 2010/1/6
Title: Functions to draw population pyramid
Author: Minato Nakazawa <minato-nakazawa@umin.net>
Maintainer: Minato Nakazawa <minato-nakazawa@umin.net>
Depends: R (>= 2.2.0)
Description: Drawing population pyramid using (1) data.frame or (2) vectors.
	 The former is named as pyramid() and the latter pyramids(), as lapper
	 function of pyramid().
License: GPL (>= 2)
URL: http://minato.sip21c.org/swtips/R.html#PYRAMID
○次に,"d:\work\pyramid\R"と"d:\work\pyramid\man"という2つのサブディレクトリを作る。Rサブディレクトリには,このパッケージに入れたい関数のコードをファイルとして置く。拡張子は.Rとすることが推奨されている。ファイル名の先頭はアルファベットでなければいけない。ここでは関数定義部分だけをpyramid.Rとして置く。
○実行例はmanサブディレクトリに入れるドキュメントファイル(拡張子は.Rdにすることが推奨されている)に含める。ドキュメントファイルは英文できちんと書かねばならないので面倒である(マークアップ言語はTeX風に\で始まっているのでとっつきやすいが)。例えば,pyramid.Rdは次のようになる(manのファイルは関数単位で用意する--そうしないとexample(pyramids)ができない--ので,pyramids.Rdも作る)。
\name{pyramid}
\alias{pyramid}
\title{Drawing population pyramid using data.frame}
\description{
  Drawing population pyramid using data.frame.
  Detailed explanation is given in Japanese
  at http://minato.sip21c.org/swtips/R.html#PYRAMID.
}
\usage{ pyramid(data, Laxis=NULL, Raxis=NULL, Cgap=0.3, Cstep=1, 
 Llab="Males", Rlab="Females", Clab="Ages", GL=TRUE, Cadj=-0.03 
 Lcol="Cyan", Rcol="Pink", Ldens=-1, Rdens=-1, main="", ...)
}
\arguments{
 \item{data}{A data.frame including left pyramid numbers in the 1st column and
 and right pyramid numbers in the 2nd column, where the numbers of males in 
 each age-class are usually given to left numbers and those of females are to
 right numbers.  If the data.frame includes 3rd column, it is used as age-class
 labels, otherwise the row.names(data) is used as age-class labels.}
 \item{Laxis}{A vector of axis for left pyramid.  If missing, automatically given.}
 \item{Raxis}{A vector of axis for right pyramid.  If missing, Laxis is used.}
 \item{Cgap}{The width of center gap (as ratio to each panel) to draw age-class. Default is 0.3}
 \item{Cstep}{The interval to write the labels of age classes. Default is 1}
 \item{Cadj}{The vertical adjustment factor for the labels of age classes. Default is -0.03}
 \item{Llab}{The label of the left pyramid.  Default is "Males".}
 \item{Rlab}{The label of the right pyramid.  Default is "Females".}
 \item{Clab}{The label of the center age-class.  Default is "Ages".}
 \item{GL}{Logical value to draw the vertical dotted lines.  Default is TRUE.}
 \item{Lcol}{The color of the left pyramid.  Default is "Cyan".}
 \item{Ldens}{The density of hatching lines (/inch) for left pyramid.
 Default is -1, when the pyramid will be filled.}
 \item{Rcol}{The color of the right pyramid.  Default is "Pink".}
 \item{Rdens}{The density of hatching lines (/inch) for right pyramid.
 Default is -1, when the pyramid will be filled.}
 \item{main}{The main title of the pyramid.}
 \item{...}{Other options.}
}
\author{Minato Nakazawa \email{minato-nakazawa@umin.net} \url{http://minato.sip21c.org/}}
\examples{
ages <- c('0-9','10-19','20-29','30-39','40-49','50-59','60-')
males <- c(34,19,11,11,8,7,5)
females <- c(26,25,16,11,7,5,1)
data <- data.frame(males,females,ages)
pyramid(data)
# another example
py.Males <- c(80,40,30,20,10)
names(py.Males) <- c('0-9','10-19','20-29','30-39','40-')
py.Females <- c(60,50,40,30,5)
names(py.Females) <- names(py.Males)
py.df <- data.frame(py.Females,py.Males)
pyramid(py.df,Llab="Females",Rlab="Males",Lcol="navy", Ldens=5, Rcol="red", 
 Rdens=10, GL=FALSE,main="An example of population pyramid\n with auto-axis")
# The census 2005 for Gunma prefecture data obtained from
# http://www.e-stat.go.jp/SG1/estat/List.do?bid=000001005042&cycode=0
# as "a00411.xls"
GunmaPop2005 <- data.frame(
 Males=c(8872, 9144, 9528, 9812, 9817, 10049, 10234, 10047, 10222, 
 10187, 10319, 10420, 10135, 10473, 10219, 10492, 10943, 11243, 10579, 
 9653, 9941, 10252, 10396, 10649, 11343, 12031, 12420, 13015, 13354, 
 14132, 14624, 15873, 16013, 15809, 15338, 14876, 14568, 14536, 14580, 
 10713, 13844, 12911, 12303, 12158, 12086, 12117, 12620, 12331, 12376, 
 13138, 14144, 13553, 14353, 15194, 16064, 17286, 18561, 18347, 18230, 
 12001, 11938, 14363, 13684, 14122, 13475, 12241, 10279, 10866, 10939, 
 10853, 10093, 9996, 9575, 9334, 9147, 8643, 8249, 7760, 7353, 6812, 
 6170, 5041, 4099, 3206, 2711, 2778, 2113, 1895, 1630, 1431, 1146, 966, 
 675, 559, 384, 263, 212, 133, 84, 38, 24, 15, 11, 6, 1, 2, 1, 0, 0),
Females=c(8323, 8750, 8964, 9359, 9559, 9605, 9511, 9800, 9790, 9848, 
 9939, 9755, 9696, 9884, 9734, 9767, 10293, 10616, 10040, 9383, 9731, 
 9748, 10211, 10204, 10500, 11086, 11758, 12248, 12548, 13524, 13907, 
 14930, 15129, 15057, 14344, 13982, 13942, 13587, 13972, 10359, 13212, 
 12435, 11934, 11673, 11668, 11583, 12100, 11867, 11917, 12746, 13364, 
 13170, 13968, 15318, 16251, 17125, 18253, 18042, 17927, 11981, 11773, 
 14450, 14124, 14438, 13502, 12960, 10729, 11710, 11697, 11884, 11413, 
 11442, 11087, 11035, 11209, 10646, 10482, 9784, 9777, 9491, 8891, 8188, 
 7636, 7034, 6123, 6103, 4577, 4415, 3861, 3426, 3035, 2571, 2064, 1683, 
 1209, 878, 739, 536, 333, 193, 134, 86, 56, 28, 13, 8, 5, 3, 1),
Ages=0:108)
pyramid(GunmaPop2005,Llab="Males",Rlab="Females",Clab="",Laxis=seq(0,20000,len=5),
 Cstep=10,main="Population pyramid of Gunma Prefecture\n (Data: Census 2005, total by gender)")
}
\keyword{hplot}
●以上のファイル群ができたらチェックしてからビルドする。pdflatexやGNUの各種開発ツールが必要だが,Rtoolsにまとまっているようである。2010年1月5日現在,Rtools210.exeをダウンロードしてインストールすれば,gccコンパイラを含むMinGWやPerlやlibpngや各種DLLなどが利用可能になるようである(2011年12月28日追記:2011年12月28日現在はRtools215.exe)。加えて,一応Microsoft HTML Help Workshopもインストールしておいた。しかし,pdfLaTeXに関しては当該ページにMikTeXへのリンクがあるが,普通にWin32版のpLaTeXをインストールしてパスが通った環境になっていれば(日本語の情報としては奥村さんのTeX Wikiが充実している),それでも大丈夫なようである。以前はMinGWのMsysという名前のbashからでないと操作できなかったように記憶しているが,2010年1月5日現在,Windows XPの場合だと,普通のコマンドプロンプトから操作可能である。
●具体的には,d:\work\pyramidに上記ファイル群を作ってある場合,d:\workをカレントにしておき,"C:\Program Files\R\R-2.10.0\bin\Rcmd" check pyramidする。コード中,TRUEを省略形でTとか,FALSEを省略形でFとか書いているところがあるとエラーを出して止まってしまうし,\item{}{}で説明を書いていない引数があると警告が出るが,それらを直せばcheckを通すのは難しくなかった。checkの結果はd:\work\pyramid.Rcheck内にlogなどが出来上がる。
●チェックを通ったらビルドである。チェックの時と同様に"C:\Program Files\R\R-2.10.0\bin\Rcmd" build pyramidとすると,d:\work\pyramid_1.1.tar.gzが出来上がった。Windows用のzipは"C:\Program Files\R\R-2.10.0\bin\Rcmd" build --binary pyramidとすればd:\work\pyramid_1.1.zipができる。なお,手動で"d:\work\pyramid.Rcheck\pyramid"をzip圧縮(+Lhacaなどでやれば簡単だろう)したものも,ローカルでのパッケージインストールには使えたが,--binaryオプションを使う方が簡単だしサイズも小さくなったのでお薦めである(2011年12月28日追記:R-2.13.0以降は,Windows環境でのzip作成はRCMD INSTALL --build pyramidとする)。
●後はCRANに登録してもらうだけなので,トライしてみた。helpを読むと,ftp://cran.r-project.org/incoming/に(sftpでは不可だそうだ),ユーザanonymous,パスワードを自分のメールアドレスにしてログインしてアップロードしてから,cran@r-project.orgにメールで知らせると,審査した上で無問題なら晴れてCRANに登録されることになるようだ。そうなれば,"Packages"の"Update packages from CRAN"でバージョン管理したアップデートができるようになるので,やってみる価値はある。FFFTPにcran.r-project.orgを設定し,PASVモードで接続して/incoming/にアップロードし,次のメールをcran@r-project.orgに送ってみた。さてどうなるだろうか。
Dear Sirs,

I have uploaded the package named `pyramid' including two functions 
to draw population pyramid.  It's very simple and small.
But, as far as I know, there is no such package ever, so that I submitted.
I hope that this package can pass the review process and go to main archive.
Thank you very much.

With Best Wishes,
Minato Nakazawa, Ph.D.
Associate Professor, Gunma University Graduate School of Medicine
●2010年1月6日の夜に,Kurt Hornikさんからの「投稿してくれてありがとう,今CRANに載ったよ」(意訳)というメールと,Uwe Liggesさん名義で自動的に生成されたものらしい,「投稿されたtar.gzからWindows用にbuildしたので24時間以内にCRANの/bin/windos/contrib/2.10/に載るよ」(意訳)というメールが相次いで届き,無事にCRANへの登録が完了した。素晴らしいことだ。
●ちなみに,example(pyramid)で表示される3つの人口ピラミッドの図では,図の下が妙に広く空いている。通常のグラフでは軸目盛や目盛ラベル,軸ラベルなどは描画領域外になるので,共通描画パラメータでpar()$marのデフォルト値がc(5.1,4.1,4.1,2.1)と広めのボトムマージンが取られているわけだが,この人口ピラミッド描画では,人数軸の軸目盛や目盛ラベルもすべて(軸ラベルは無いが)描画領域内であるため,par()$marの値がデフォルトのままだとボトムマージンが広すぎるのだ。この問題を解決するには,関数定義内部ではparはまったく弄っていないのだから,図を描かせる前に,例えばpar(mar=c(2,4,4,2))などとやっておけばいい。その状態で描かせた2005年の群馬県の人口ピラミッドが下図のようになる。
Age structure of Gunma prefecture by gender, in 2005, based on the national census record.
●別の手としては,関数内部でxpdをTRUEにしておき(これもpar()のパラメータだが,これくらいはたぶん,関数内部で最初にop <- par(xpd=TRUE)として最後にpar(op)するという常套手段を使えば問題ないだろう),軸目盛や軸ラベルを最初に設定した架空座標系外に描かせてしまうという手もあったかもしれない。しかし"The simple is best."という観点から,現在のような仕組みにした。
●(2010年3月11日更新)1月末にPhilippe Guillet, MD.という方からメールをいただいて,年齢階級の文字列だけをサイズを変える手段と,人口の目盛りラベルの書式を指定する手段を提供しようと考えていたのだが,漸く実装できたので,バージョン1.2とした。前者は年齢階級文字列をプロットする際のtext()のオプションでcex=で指定すればいいので,それをpyramid()のオプションCsizeで与えることにした。デフォルトを1にしておけば,指定しない場合はこれまで通りの動作になる。後者は目盛りラベルをプロットするtext()の中で,ただpaste(数字)としていたのを,paste(formatC(数字,format=AxisFM,big.mark=AxisBM,big.interval=3))とすることで可能であった。AxisFMを"g"としておけばこれまで通りだし,"d"としてAxisBMを","にすれば3桁ごとにコンマが入った書式になるし,AxisBMはデフォルトでは""なので,AxisFMを"d"とするだけなら桁の多い数字でも浮動小数点表示ではなく普通の整数としてプロットすることができる。このオプション追加に伴って,群馬県の人口ピラミッドをプロットする例で使い方を示すようにした。
●(2012年3月15日更新)サイト移転に伴い,ドキュメントのURLを書き換えた他,GunmaPop2005をドキュメントのexampleからコードの方にうつしたので,pyramidパッケージをロードするだけでGunmaPop2005が使えるようになった。CRANからもinstall.packages("pyramid")でインストールできるが,pyramid_1.3.zip,またはpyramid_1.3.tar.gzとして,このサイトからダウンロードすることもできる。
●(2014年9月4日更新)鵯記にメモした通り,フレーム内に描画するバージョンとして関数pyramidf()を追加した。オプションとして,描画フレームの左下と右上の座標を指定するframe=が増えただけで,他はpyramid()と同じである。frame=のデフォルトはc(-1.15, 1.15, -0.05, 1.1)であり,これはデフォルトのCgap値0.3でpyramid()を実行したときの座標と同じである。バージョン1.4をインストールすれば,library(pyramid); example(pyramidf)とするだけで,2005年と2010年の群馬県の人口ピラミッドを重ね描きする実行例が表示される。CRANからもinstall.packages("pyramid")でインストールできるが,pyramid_1.4.zip,またはpyramid_1.4.tar.gzとして,このサイトからダウンロードすることもできる。
●(2014年9月6日追記)ついでに,神戸市の地図上に区別に人口ピラミッドを描くサンプルコードも載せておく。国土地理院のウェブサイトから兵庫県のshapeファイルを含むN03-140401_28_GML.zipをダウンロードし展開した後,EpiInfoに含まれているEpiMapのAddLayerPartial機能を使って神戸市だけの地図情報をkobe-city.shpとして保存し,同じディレクトリに神戸市の区別の人口データ(神戸市が独自に国勢調査結果をまとめているページから,平成22年国勢調査結果の「全市」をダウンロードして加工したもの)を置いて,Rでそのディレクトリを作業ディレクトリにしてこのコードを実行すると下図ができる(低解像度ラスタデータなので,国土地理院のシェイプファイル地理情報を使って描いたが個別許可は必要ないはず)。
神戸市区別人口ピラミッド

計算テクニック,他

追加パッケージを使うには?
●プルダウンメニューから選んでもいいが,library(パッケージ名),あるいはrequire(パッケージ名)とする方が簡単である。library()とrequire()の最大の違いは返り値。
●未インストールのパッケージをCRANからインストールするには,install.packages("パッケージ名". dep=TRUE)とするのが基本。既にインストールしてあるパッケージを最新版に更新するには,update.packages(ask=FALSE, checkBuilt=TRUE)とするのが簡便である。
●CRAN以外に存在するパッケージをインストールするには,レポジトリを適切に指定すること。
●既にインストール済のパッケージのバージョンなどを確認するには,packageDescription("パッケージ名")またはlibrary(help=パッケージ名)とする。
標準化偏回帰係数を得る
●ベクトルとして計算させると楽である。例えば,従属変数をy,独立変数をx1, x2, x3とすれば,
res <- lm(y~x1+x2+x3)で線型回帰を行った後,
sdd <- c(0,sd(x1),sd(x2),sd(x3))として各独立変数の不偏標準偏差ベクトルを作り(0は切片用),
stb <- coef(res)*sdd/sd(y)として,偏回帰係数ベクトルに不偏標準偏差ベクトルを掛けて,従属変数の不偏標準偏差で割ってやれば,stbに標準化偏回帰係数のベクトルが得られる。
変数選択
●Kendallも書いているように,線型モデルにおける機械的なステップワイズの変数選択は,あまり薦められる方法ではなく,できれば,それ以前に先見的知識から独立変数群に入れられる変数を選んで作ったモデルをそのまま解釈するか,またはSASでいうMAXR法のように,K個の独立変数群候補の中からJ個(J<K)の,R2を最大にする変数の組を総当りで計算し,尤度比検定によって尤度が有意に変わらない最小の独立変数群を選ぶとかすべきと思う。
●しかし,いろいろな事情によって,変数選択をしなければならないこともあるだろう。そういう場合,Rでは,標準で組み込まれているstep()関数を使うことで,AICによる変数選択をしてくれる。step()に渡すモデルはlm()でもglm()でもよい。directionというオプションで,増加法(direction="forward"),減少法(direction="backward"),増減法(direction="both")を指定することができるが,scopeとして1つの式(upper扱いされる),あるいはupper(独立変数候補のすべて)やlower(常にモデルに入れる独立変数)からなるlistを指定すれば変数増減法(あるいは減増法)になるし,scopeを指定しなければ(その場合,モデルとして与えた式そのものをupperのscopeと解釈してくれる)変数減少法になる。MASSライブラリに入っているstepAIC関数は,正確なAICを計算してくれるとのこと。
分類変数(因子,factor)について
●身長heightを150cmから180cmまで5cm刻みのカテゴリ変数hcに代入したいとき,
hc <- cut(height,seq(150,180,by=5))とすればよい。例えば最小区間が(150,155]となる。ただし,このままでは分類変数(factor)に使えないので,
hcf <- factor(hc)とすれば,hcfは分類変数になる。また,hco <- ordered(hc)とすれば,hcoは順序のある水準になる。
●この場合,lm(y~hc)lm(y~hcf)lm(y~hco)はすべて結果が違うし,plot(y~hc)ではドットがプロットされるが,plot(y~hcf)plot(y~hco)では層別箱ヒゲ図になる。
●2つの要因で層別した同時散布図をcoplot(y~x|a*b)によって実行する場合,aやbが分類変数でないと層別数はaについてもbについてもデフォルトでは6である(numbers=で変更可)。aやbがfactorまたはorderedなら,カテゴリごとにplot(y~x)がなされる。
バッチ実行(Windows2000環境の場合)
●Rterm.exeコマンドにリダイレクト入力することが必要。つまり,例えば,"C:\Program Files\R\rw1081\bin\Rterm.exe --no-save < %1"という内容のバッチファイルrun_R.batを作って*.Rに関連付けしておいて*.Rをエクスプローラでクリックするとか,このバッチファイルをパスの通ったディレクトリにおいてrun_R foo.Rとすれば,Windows2000の場合,コマンドラインでfoo.Rが実行される。しかし,バッチ処理が終わると同時にRが終了してしまうので,関連付けからの実行でも画面への出力を残しておくには,バッチファイルの2行目に"pause"を入れておく必要がある。あるいは,画面に残さなくてもいいなら,"C:\Program Files\R\rw1081\bin\Rterm.exe --no-save < %1 > _Rresult.txt"としておけば,結果の出力先は画面でなく,カレントディレクトリの_Rresult.txtというファイルになる。グラフィック出力はそれでも残らないが,そもそもwindows()をせずにグラフィック出力をバッチ実行させるとカレントディレクトリのRplots.psというポストスクリプトファイルに書き込まれるので,それを加工してもいいし,画面でみつつpdfで残したければwindows()してからplot()などの描画コマンドを実行し,描画後にdev.copy(pdf,"_Rresult.pdf"); dev.off()するようにバッチ実行するRスクリプト中に書いておけばいい。画面で見なくてもよければ,例えばpdf("_Result.pdf",horizontal=F)とdev.off()で描画コマンドを挟めばよい。
データの再コーディング
●例えば,9つの地域区分がxというデータフレームのAREAという数値変数(値は1~9)に入っているのだが,それを地域名(A~I)がついた分類変数に変換し,かつ3種類の街区(市街地=A,C,G,農村部=B,F,H,工業地区=D,E,I)に区分し直した新しい分類変数REGも使いたいとする。次のコードで可能である([2004年2月18日訂正]……と書いてずっと公開していたが,実はうまく動かなかったことをお詫びする。AREAは数値変数だから分類変数を1つ付値した段階でそれ以降が欠損になってしまうのだった。別の変数を作る限りは大丈夫なので,新しい変数REGを作るだけにしてAREAはそのままにすれば問題なく動く)。
NREG <- c('工','農','市','農','市','市','市','農','農')
x$REG <- rep(0,length(x$AREA))
for (i in 1:length(x$AREA)) {
 if (!is.na(x$AREA[i])&(x$AREA[i]>0)&(x$AREA[i]<10)){x$REG[i] <- NREG[x$AREA[i]]}
 else {x$REG[i] <- NA}
}
x$REG <- as.factor(x$REG)
●(2004年2月18日追記)実はここに載せていたコードではうまく動かないことは群馬大学の青木先生から御指摘いただいたのだが,ついでにもっとシンプルでうまく動くコードをお教えいただいたので許可をいただいて下に掲載する(1行目はデータフレームxに1から9までの整数値をもつ変数AREAを付値している)。ifelse()関数をこんな風に使えるとは驚いた。
x <- data.frame(AREA=1:9)
NAREA <- c('A','B','C','D','E','F','G','H','I')
NREG <- c('工','農','市','農','市','市','市','農','農')
ifelse(!is.na(x$AREA) & (x$AREA>0) & (x$AREA<10), {x$REG <- NREG[x$AREA]; x$AREA <- NAREA[x$AREA]}, x$REG <- x$AREA <- NA)
●(2004年3月14日追記)上のAREAの書き換えのように,水準数が変わらずに名前だけ変えるときは,下記のようにlevels()関数を使うと楽である。
x$AREA <- as.factor(x$AREA)
levels(x$AREA) <- c('A','B','C','D','E','F','G','H','I')
Rの文字列操作について(2014年6月1日)
●R本体が最も苦手とする処理の1つが文字列操作である。
●read.delim()関数などで,文字列をファイルから読み込むとき,通常は自動的に因子型になる。この自動変換をさせないグローバルオプションがoptions(stringsAsFactor=FALSE)である。ちなみに,データフレームでは列ごとに変数の型が異なっていてもいいが,matrix()ではすべての要素が同じ型でなければならない。
●既にデータフレームにしてしまった後で因子型の変数を一括で文字列型に変えたい場合は,Stack Overflowのスレッドにあった(bobがデータフレームだとして),

i <- sapply(bob, is.factor)
bob[i] <- lapply(bob[i], as.character)

というコードが便利だと思う。
●数値を書式付きで文字列にするのは,C言語と同様の仕様で,sprintf()という関数が使える。表示桁長を見やすく揃えるときも便利。例えば,sprintf("%09d", 4)の結果は,
[1] "000000004"と表示される。
●文字列処理関数としてよく使われるのは,paste(),substr(),strsplit()などであるが,あまり機能は充実していない。stringr()パッケージを使うと,例えば,ある文字列に含まれる別の文字列の個数を返すstr_count()関数などが使える。例えば,str_count("abc1234def5432", "4")は,第二引数の文字列が第一引数に2回出現するので2を返す。
Rで数学をやってみる
■川端裕人『算数宇宙の冒険 アリスメトリック!』実業之日本社,ISBN 978-4-408-53563-0(Amazon | bk1)を読み,解説を書かれている柏野雄太さんのサポートページを見て触発されたので,Rでリーマン・ゼータ関数まわりのことがどれくらいできるのか試してみた。
■とりあえずゼータ関数を探してみたら,VGAMライブラリ[manual pdf]とgslライブラリ[日本語版マニュアル……このgslというライブラリはGNU Scientific LibraryをRに実装したもののようだ]に入っていた。どちらも関数名はzeta()であった。
■しかし,VGAMでもgslでもリーマン・ゼータ関数は整数あるいは実数に対してしか定義されていないので,解析接続(注:この概念が今一つ良く理解できていないのがいけないんだけれども)ができず,柏野さんのページで示されているような図は作れない。Rでも複素数(complex)の四則演算などはできるのだが(例えば,x <- 1+1iと打ってから,x*2と打てば[1] 2+2iと表示されるし,x^2と打てば[1] 0+2iと表示される),zeta()は複素数を引数にとれず,gslでは虚部が捨てられて計算されてしまうし,VGAMではNAになってしまうのだった。
■実数範囲でできることは,例えばζ(2)=π2/6[式1]を確かめてみるとか,ζ(2),ζ(3),...,ζ(n)とnを増やしていったときの挙動を図示してみるといったことである。これは,Rで普通に試せる。
■最初の例は,ζ(2)=1+1/(2^2)+1/(3^2)+...という無限級数だから,とりあえず一万項くらいやってみるには,
sum(1/(1:10000)^2)
とすればいい。出力結果は,[1] 1.644834となる。一方,円周率はpiとして組み込まれているので,[式1]の右辺は,
pi^2/6
でいい。出力結果は,[1] 1.644934となる。つまり,ちょうど一万分の一の誤差があるわけである。誤差が見えなくなるまで項数を1桁ずつ増やしていったら,千万項が必要なことがわかった。つまり,
sum(1/(1:10000000)^2)
の結果は,[1] 1.644934となった。なお,後で青木先生からご指摘を受けた通り,コンピュータでの数値計算の常識として,小さい方から足さないと計算誤差が出るので,足し算は小さい方からやるべきと考えたら,
sum(1/(10000000:1)^2)
と書く方がいい。けれども,実はRではどちらでも小さい方から計算されているようだ。つまり,明示的に繰り返しループを書くと,
a <- 0; for (i in 10000000:1) a <- a+1/i^2; a
とした方がいいということだ。この場合はoptions(digits=20)くらいにしないと違いがわからないのだけれども,確かに
a <- 0; for (i in 1:10000000) a <- a+1/i^2; a
とすると,結果に微妙な誤差が出る。
■(2009年12月11日追記)黒川信重編著『数学セミナー増刊:リーマン予想がわかる』日本評論社の巻頭,桜井進「高校生からわかる――超入門:リーマン予想」に書かれていた,オイラーが求めた収束が速い関数というのは,本当に収束が速かった。つまり,
sum(1/((1:50)^2*2^(0:49)))+log(2)^2
の出力結果は,options(digits=20)のとき[1] 1.644934066848226となって,完璧にpi^2/6と一致した。
approxzetatwo <- function(x) { sum(1/((1:x)^2*2^(1:x-1)))+log(2)^2 }
と関数定義してから,sapply(1:50,approxzetatwo)とすると,実は41項目で同じ結果が得られていることがわかる。
■計算する項数を増やすにつれて級数がpi^2/6という極限に近付いていく様子を図示するには,例えば,
plot(cumsum(1/(1:1000)^2),type="l",log="x")
lines(c(1,1000),c(pi^2/6,pi^2/6),col="red")
とすれば下図が描ける。ここでは,見やすくするため項数を対数軸表示にした。なお,pi^2/6は,library(gsl)の後でなら,zeta(2)と置き換え可能である。参考までに書いておくと,ここまでの内容は,いつの間にかMaximaと融合し,Gnuplotも内蔵した高機能なものになったEulerというソフトでも,同じような感じで実行できる。
ζ(2)への漸近
■続いて,ゼータ関数の引数を2以上の整数で徐々に大きくしていくとどうなるかを試してみる。もちろん,賢いgslライブラリを使えばplot(2:1000,zeta(2:1000),type="l",log="x")などということも簡単にできるが,定義式通りにx乗分の1を計算させてzeta()並みの精度を出すために一千万項まで計算させてしまうと,ζ(100)まででもハングアップしたかと思うくらいの途方もない時間がかかるのでやめておく方が無難である。
■そこで,計算させる項数を1万項くらいにとどめてζ(2)からζ(100)までの近似をプロットさせるには,べき乗関数powを定義してから,sapplyの結果の行方向の和をとればいいので,Rのコードは次の通りである。
pow<-function(x,y) x^y
plot(2:100,rowSums(sapply(1:10000,pow,-(2:100))),type="l",log="x",xlab="ζ(n)の引数n",ylab="ζ(n)")
require(gsl)
lines(2:100,zeta(2:100),col="red",lty=2)
legend("topright",lty=c(1,2),col=c("black","red"),c("近似","ζ"))
ζ(2)~ζ(100)
■上図のように横軸に2から100まで,縦軸にそれに対応するゼータ関数の近似値をとったグラフが描けるわけである。この状態でlines(2:100,zeta(2:100),col="red",lty=2)とすると,gslライブラリのゼータ関数で計算した値が重ね描きされる。パソコン画面くらいの解像度では,2つのグラフのズレはほとんどわからないくらいである。
■さていよいよ,複素数を引数にとれるゼータ関数だが,webを検索していたら,Javaで書かれたものが見つかった。コンパイル済みのもの(*.jar)しか公開されていないようだが,このソフトのソースを何とか入手してC++に移植すればRから呼べるかもしれない。しかし公開されていないわけだから,たぶんそれは難しいだろう。それで,このページに書かれていた式をそのままRに置き換えてZeta(s)を関数定義したら(但しもちろん無限和はできないので柏野さんに倣って500項の和にして),複素数を受け付けるゼータ関数ができた。もっとも,Zeta(-1)が-0.0833313となってしまうし(本当は-1/12だから-0.08333333となるはずだし,不思議なことにgslライブラリのzeta(-1)はそうなる。なぜ不思議かといえば,ζ(-1)=-1/12は解析接続をしたことで得られる解だからである),Zeta(2)は1.644930074848となってしまうので近似精度は低い。しかも,Rのガンマ関数gamma()が複素数を引数にとれないので実部が負の場合はまだ計算できない。
□それでも,実部が正の場合,例えばZeta(1/2+1i)は計算できて, 0.144659-0.6974934iとなる。もし,これが正しければ,柏野さんのサポートサイトに載っていたグラフが作れるはずである。
□複素ガンマ関数は,わざわざ作らなくても,gslライブラリにlngamma_complex()という対数をとった形で入っていることがわかった。exp()は複素数を通すので,これで複素数を通すζ関数の近似ができるはずである。また,大浦さんという方が開発されたライブラリにはCとFortranで複素ガンマ関数のコードが提供されているので,これをRに移植すればgslライブラリを使う必要がなくなるように思う。ぱっと見た感じではわかりやすいソースコードなので,そう苦労しなくても移植できそうな気がした。
■大浦さんのコードをRに移植し,複素数を受け付けるゼータ関数改良してみた。500項の和だったところを1万項の和にしたら,Zeta(2)は1.644934となってくれたのだが,Zeta(-1)は実部は-0.0833333となってくれるものの,なぜか虚部がついてしまう。もう少し検討が必要か。
□試しにgslライブラリを使って,Zeta()の関数定義中の複素変数のガンマ関数をexp(lngamma_complex(t))とした関数Zeta2()を定義したら,Zeta2(-1)は[1] -0.08333333+0iとなったので,移植したコードで虚部がついてしまったのは,計算精度の問題だったようだ。
■柏野さんのサイトで言及していただいたのにできていないのは悔しいので,とりあえずgslのlngamma_complex(t)を使ったZeta2()で図示してみる。
source("http://minato.sip21c.org/swtips/Zeta.R")
AbsZeta <- function(x) Mod(Zeta2(1/2+x*1i))
curve(AbsZeta(x),0,50)
lines(c(0,50),c(0,0))
とすると,条件判定のところで警告メッセージが出る部分があるが,下図ができる。残念ながら計算精度が低いのか,非自明なゼロ点がきちんとゼロになってくれないのだが,だいたい近い形のグラフにはなる。
Abs(zeta(1/2+t*1i))のカーブ
■(2010年5月19日追記)大きな勘違いをしていたのだが,上記欠点は,計算精度が低いというよりも,curve()だとステップが粗いため,非自明なゼロ点を与えるxが計算されていないということが原因だった。また,条件判定の警告は,Zeta2()関数がベクトル化されていないためであった。一方,Zeta2()を定義する部分は,100000項でなくて10000項でもZeta2(-1)が[1] -0.08333333+0iとなってくれたし,1000項でも[1] -0.08333328+0iとなった。以上を踏まえ,以下のように書き換えたコード(Zeta2.R)を試してみた。
require(gsl)
Zeta2 <- function(s) {
 if (s==1) { ret <- Inf } else if (Re(s)>0) {
  a <- 1/1^s
  for (i in 1000:2) { a <- a + ((-1)^(i-1))/(i^s) }
  ret <- 1/(1-2^(1-s))*a
 }  else {
  t <- 1-s
  ret <- 2*(2*pi)^(-t)*cos(t*pi/2)*exp(lngamma_complex(t))*Zeta2(t)
 }
ret
}

AbsZeta <- function(x) Mod(Zeta2(1/2+x*1i))
x <- 0:50000/1000
y <- sapply(x,AbsZeta)
plot(x,y,type="l")
lines(c(0,50),c(0,0),col="red")
かなりの計算時間がかかったが,下図の通り,ある程度満足のいくグラフになったように思う。速いコンピュータがあれば,Zeta2()の計算の項数を増やしたりAbsZeta()を描画するステップを細かくすれば,もっときれいなグラフが描けるだろう。
Abs(zeta(1/2+t*1i))のカーブ・改良版
続・Rで数学を(2014年6月1日)
■滅多に使わないため,自分でもよく忘れるので,このページのアーカイヴからここにコピペしておく。
(多倍長演算)Rで多倍長演算をサポートするパッケージは,gmpが便利である。gmpパッケージをロードすると,as.bigz()で多倍長整数,as.bigq()で多倍長実数を扱える(過去のメモ参照)。しかし,高精度計算と表示には,Rmpfr(R Multiple Precision Floating-Point Reliable)というパッケージの方が向いている。πなどの定数も高精度で含まれているからだ(pi,gamma,catalan,log2については,Const("pi",100)などとすれば高精度で得られる)。mpfr()という関数で高精度オブジェクトができ,そのまま計算ができる。なお,Const("pi",100)の結果はasin(mpfr(1,100))*2と一致するが,mpfr(pi,100)とは一致しないことには注意が必要である。mpfrになった値に関数を適用すると自動的に高精度計算になるが,逆では単なる型変換になってしまうので下位の桁は適当に補われるだけになる。100桁のexp(1)が欲しい時は,mpfr(exp(1),100)ではなくて,exp(mpfr(1,100))でなくてはいけない。

リンクと引用について

Correspondence to: minato-nakazawa[at]umin.net.