2007.4.24

紙に書かれたマニュアルを見たい場合には、 参考文献[1]の第I分冊の前半に目を通すと概要が分かります。 第II分冊は関数のマニュアルとなっています。
2. Splusとの対話
Unix 環境で次の最初の1行を入力して下さい。
% Splus
S-PLUS : Copyright (c) 1988, 2002 Insightful Corp.
Copyright (c) 2002 Mathematical Systems, Inc.
S : Copyright Lucent Technologies, Inc.
Version 6.0.4 Release 1 for Sun SPARC, SunOS 5.7 : 2002
Working data will be in /home/cs/seki/MySwork
>
となって、Splusの環境に入れるはずです。最後の1行は Splus のプロンプトで、
ここに、何かを入力すると、入力した文字列が S言語として解釈できれば、
解釈した結果を次の行以下に表示して、プロンプト表示状態になります。
次の1行("2+3")を入力して下さい。
> 2+3 [1] 5 >このように普通の四則演算(
+,-,*,/,^)は素直に実行してくれます。
また、%/%, %%で整数除算の商と余りを計算します
( (7%/%3) は 2、 (11%%4)は 3となる)。
2行目行頭の[1] については、次節で説明します。
仕事を終えて、Splus から出て、UNIX環境に戻りたい時には、
> q()と入力します。以上、Splusを起動して終了するまでを、1つのセッション(session)と 呼びます。
S言語は怠惰な表現の許される言語で、適当に思いつくやり方で表現しても、 何とかしてくれます。何でもやって見て下さい。
3. オブジェクト
原始オブジェクト
解析対象となるデータの基本構成要素(原始オブジェクト)には、
論理値(TURE, FALSE [T,Fと省略可])
、整数、実数、(半角)英数字などがある
。
これらは、この順に後ろのものほど一般表現となる。
ベクトルと行列
同じ種類のものは、c( ) と言う関数でベクトルにまとめられる。
> c(1,2,3) [1] 1 2 3S言語の演算子に"
:" がある。初項と末項から公差1の等差数列を作って
くれる
(関数 seq()も参照)。
> 101:121 [1] 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 [20] 120 121S が返した2行目の
[20] は、先頭要素が表示ベクトルの 20番目の要素である
ことを示す。
ベクトルの一般形として、行列もオブジェクトとして扱える。
matrix() と言う関数で作れる。
> matrix(1:6,nrow=2) # 1から6の数列を行数2の行列とせよ。
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
<-''
(2字の間に空白を含んではならない)、
または ''_''または ''=''を使う。
> x <- 3:10 > x [1] 3 4 5 6 7 8 9 10Sは この例のようにオブジェクト名だけを入力すれば、その内容を表示する。
なお、Sのオブジェクトを表す変数 (上の例の"x"などオブジェクト名) に(原則として)変数宣言/型宣言は必要ない。 任意の変数は、付値されたときに、付値された内容と型を持つ。
注意 : あなたが作成したオブジェクト(特に関数)が既成の関数と同じ名称であると、
その関数は参照できなくなる(使えなくなる)場合がある。
関数が一つ使用不能になると、その関数を使っている関数も使えなくなることがある。
特に、Splusで重要な関数 "c"や "t"を
上書きしないように注意のこと。
上書きしてしまった時には、次の"rm()"で
自分のオブジェクトを消去すれば良い。
自分が定義したオブジェクトは"ls()"で、
そのリストを見ることができる。
また、"ls()"でリストされるオブジェクトは
不要ならば、"rm("その名前")"とすることで、消去できる。
ls()"で見られるオブジェクトのしまい場所は
session frame と呼ばれる。その他にも、
各種のライブラリが用意されている。また、関数呼出しによって生成される
frame もある。
> x [1] 3 4 5 6 7 8 9 10 > x+2 [1] 5 6 7 8 9 10 11 12 # x の各要素に2を足す。 > x^2 [1] 9 16 25 36 49 64 81 100 # x の各要素を巾乗する。 > y <- 1:4 > y [1] 1 2 3 4 > x+y [1] 4 6 8 10 8 10 12 14 # y の 長さが足らないので y を2回複製して、xに足し込んだ。 # つまり、3+1 4+2 5+3 6+4 7+1 8+2 9+3 10+4 を実施した。 > z <- 1:3 > z [1] 1 2 3 > x+z Problem in x + z: length of longer operand (8) should be a multiple of length of shorter (3) Use traceback() to see the call stack # z の 長さが足らないので z を複製して、xに足し込もうとしたが、 # 整数回の複製でピッタシにならなかったので、意図不明でエラーになった。なお、version 4.x 以前では、以下のように警告を出しながらも実行していた。
[1] 4 6 8 7 9 11 10 12
Warning messages:
Length of longer object is not a multiple of the length of the shorter object
in: y + x
# つまり、3+1 4+2 5+3 6+1 7+2 8+3 9+1 10+2 を実施した。
# 最後に複製した方が余ったので、本当にこの計算をしたかったのか
# Sは貴方の意図が不安になって、Warning messages を出してくれている。
行列についても、同様のルールが適用される。 なお、次元数が異なるときには 1次元に直した順で演算される。
> xx <- matrix(1:6,nrow=2)
> xx
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> z
[1] 1 2 3
> xx+z
[,1] [,2] [,3]
[1,] 2 6 7
[2,] 4 5 9
# z の 個数が足らないので z を2回複製して、xに足し込んだ。
# S では配列を1次元に解釈する時、前の添字が忙しく回るように
# 並べ直す原則を採用しているので、たとえば[2,1]要素は 2+2
# (x[2,1]+z[2])をしている。
論理演算子と参照
論理演算子は
& | ! && ||
を
(前二つは、論理値のベクトルを扱える。
後ろ二つは、C言語同様に先頭からの条件つき比較をし、単一の論理値を返す。
関数 ifelse() も参照)、
比較演算は
< > <= >= == !=
を使う。ベクトルなどに関しても、これらの比較演算ができ、たとえば、
先ほどの xについて、以下となる。
> x %% 2 ==1 # 奇数かどうかの論理値ベクトル [1] T F T F T F T Fこの表現では、ベクトル(
x)と
単一値( 2や 1)が
混在していて論理的にはおかしいが、S言語のベクトル演算では
必要に応じて短い方のベクトルを適当に反復して複製して、
長さを揃えてから演算してくれる原則がここでも適用されている。
ベクトル、行列を表すオブジェクトの要素は次のように参照できる。
> x[2] # 2番目の要素の参照 [1] 4 > x[2:4] # 2から4番目の要素の参照 [1] 4 5 6 > x[ x %% 2 == 1] # 値が 奇数であるものの参照 [1] 3 5 7 9
> c( T, 1, 2.4, "a", matrix(1:6,nr=2)) [1] "TRUE" "1" "2.4" "a" "1" "2" "3" "4" "5" "6"となって、もっとも一般的な原始オブジェクト(この場合、文字型)に強制変換され、 配列などの構造を潰して、一つのベクトルにしてしまう。
このようなタイプの異なるものを一つにして扱いたいときには
list() と言う関数を使い、リストにすると、
元の構造を保存したまま、ひとつのオブジェクトとできる。
> L <- list( T, 1, 2.4, "a", matrix(1:6,nr=2))
> L
[[1]]: #一番目のリストの要素の
[1] T # ベクトルの一番目の要素は 論理値 T
[[2]]: #二番目のリストの要素の
[1] 1 # ベクトルの一番目の要素は 整数 1
[[3]]: #三番目のリストの要素の
[1] 2.4 # ベクトルの一番目の要素は 実数 2.4
[[4]]: #四番目のリストの要素の
[1] "a" # ベクトルの一番目の要素は 文字 a
[[5]]: #五番目のリストの要素は
[,1] [,2] [,3] # 整数の配列
[1,] 1 3 5
[2,] 2 4 6
ベクトル同様に、部分参照が出来るが、リストの場合には
二重カギカッコを使う。
> L[[2]] # リストの2番目の要素は 整数 1
[1] 1
> L[[5]] # リストの5番目の要素は 配列
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
リストに対して一重カギカッコも使えるが、この場合部分リストの抽出を
することになる。結果は単一の要素の指定でもリストである。
> L[2:3] # リストの2、3番目要素の部分リストは [[1]]: [1] 1 [[2]]: [1] 2.4 > L[2] # リストの2番目要素の部分リストは [[1]]: [1] 1また、リストの要素には名前を付けることができる。
> L <- list( L=T, N=1, C="a") > L $L: [1] T $N: [1] 1 $C: [1] "a"この場合は、以下のように
$を使って、名前でもリストの要素を参照できる。
> L$C [1] "a"
4. 関数
既成の関数
関数は関数名の後に ( )で括った引数を並べて参照する。
> rnorm(10, 1, 2) # 平均 1, 標準偏差 2 の正規乱数 10個 [1] 1.8253494 -2.0670295 0.8855791 4.6236200 0.9862706 -1.1668234 [7] -3.1592610 6.6921252 -0.5139193 -0.6347702引数がないときにも
( )は省略できない。省略すると、
関数自身の定義を表示する (関数定義も一種のオブジェクトです)。
> ls() # S 環境にあるオブジェクトのリストを表示せよ。 > ls # 関数オブジェクト ls の内容を表示せよ。
引数には名前がついているので、参照時に名前を記述すれは順番を変更してよい。 参照時に名前をつけずに引数並びに並べられた引数は、 名前指定で対応づけられた残りの関数定義中の引数と、引数並びでの順番を頼りに、 対応づけられる。 また、引数の名前は区別できる限り省略してよい。以下は、皆同じこととなる。
rnorm(10, mean=1, sd=2) rnorm(10, sd=2, mean=1) rnorm(10, s=2, m=1) rnorm(10, s=2, 1) rnorm(10, 1, 2)
引数には省略時の値(デフォルト)が指定されているものがある。 たとえば、先ほどの関数 rnorm は次の1行目のように定義されている。 よって、 rnorm(10)と参照しても、上述の例と同じ結果になる。 なお、デフォルトの指定されていない引数は省略できない。
> rnorm # rnormという関数は
function(n, mean = 0, sd = 1) # 左のようにn,mean,sd という引数を持つ関数で
# mean,sd はデフォルトとして0,1が設定されている。
.Internal(rnorm(n, par1 = mean, par2 = sd), "S_rfuns", T, 9)
# 実体は内部関数の呼び出しとなっている。
既存の関数についての疑問は 関数名が分かるものについてはオンラインマニュアルで調べるとともに、 その実体をみてみると良い。 たとえば、mean (平均を求める関数) の定義を表示させ、 ?mean で表示されるマニュアルの内容がどう実現されているのか、読取ってみよう。
> help(mean) # mean という関数のヘルプをみせよ。
# unix のコマンド more とおなじ環境、キー"q"で出られる。
> ? mean # (上の行に同じ)
> mean # 関数のオブジェクトとしての実体をみる。
最後の例のようにして、関数オブジェクトの実体をみると、
先頭行のfunctionの後ろから、
その関数がどんな引数を用い得るかがわかる。
また、その下の関数の中身をみることで、手続きとして何をしてくれるのかがわかる。
> x <- rnorm(100) # 正規乱数を 100 個 生成 > mean(x) # その平均値をもとめる。という調子です。上の結果を簡潔に求めるには、
> mean(rnorm(100))としても結構です。ただし、途中経過に自信がないときには、 前者のように一つずつ付値して、途中の内容を確認してゆく方が、 結局早道であることもあります。
my.log <-
# 自家版対数関数(関数名 my.log。負の数の場合0を返す対数関数)
function(x=1) # function(引数名1=引数1の省略値, 引数名2=引数2の省略値, ... )
{ # 引数を使った関数本体始まり
if( x>0 ) y <- log(x)
else y <- 0 # y の計算
y # この関数の返し値
} # 関数本体の定義の終了
上記の関数を 1行ずつコマンドラインから入力してみよ。
( 各行を入力すると、
2行目以降はプロンプトが" > "から" + "
に変わる。これは、まだ、その前の行での入力での
構文が閉じていないので、S が続きを要求している印である。
また、S 言語の各行のうち、# 以降の内容はコメントであると見なされて
構文解析されない。)
その後、ただfun.exと入力すると、
入力した通りの関数オブジェクトが実現したかどうかを確認出来る。
> my.log
# 自家版対数関数(関数名 my.log。負の数の場合0を返す対数関数)
function(x = 1)
{
# 引数を使った関数本体始まり
if(x > 0) y <- log(x) else y <- 0
# y の計算
y
}
実行するには、下のように引数を具体的に括弧内に指定して実行してみれば良い。
> my.log(2) # 関数の引数を 2としての参照 [1] 0.6931472 > my.log(-1) # 関数の引数を -1としての参照 [1] 0 > my.log() # 関数の引数を 省略しての参照 [1] 0 # x はデファオルトの1として計算される
長い関数を定義したいときには、適当なエディターのウィンドーを新らたに開いて、 その上で、関数をエディットして、マウスで Sのウィンドーへコピーするのが早い。 ( コピー内容が沢山の行にわたる場合、コピーに失敗する可能性が高くなるので注意する。 本当にとても長い関数を定義したいときには、 source("ファイル名") コマンドなどで、 定義を記述したファイルを S環境に直接読み込む。 ) または、emacs で M-x shell を実行して、 emacs 上で Splus を実行すると便利。
関数の返し値は、その関数本体の最後の文の値となる。
返し値として複数のオブジェクトを返したいときには、
それらのリスト
(同種の数値ベクトルなら配列でもよい。複雑な構造を返したいときには
"クラス"を定義してそれを返すことも行われる。)を
返すのが標準である。
また、途中で関数を終了して、値を返すために return(返し値)がある。
制御構造としては、for, while, repeatも揃っている。
> x <- 1:10 > rslt<-NULL; for(a in x) rslt<-c(rslt,ifelse(a %% 2 ==1,"o","e")) # x の要素が偶奇に応じた"o","e"のベクトルを作る。ここでは、
rsltを空
( NULL )として用意して、その後、
x の要素を一個ずつ a へ取り出し、その偶奇に応じた文字を rsltに連結している。
後で要素を反復的に追加してゆこうとする初期値(オブジェクト)を
用意するためにはこのように NULL を用いる。
セミコロン ; で結ぶことにより、複数の文を1行中に書ける。
また、制御構造のブロックを作るためには中括弧 { }を使う。
また、反復中に、反復を終了させる breakや、今回の反復を中止して、
次の反復へ進む nextもある
( breakや nextは、
使わないで済む算法を考えた方が良いが。いやその前に、
適切な配列などのデータ形式を利用していれば、
反復構文自体を使わないで書ける場合が多い。
( ifelse(x %% 2 ==1,"o","e") の結果は何になるか?)
そのためには
ベクトル/行列演算、
論理演算子と参照などの機能や、applyなどの
データ集約
にまとめた適当な関数などを用いれば良い。
(最初のうちは簡潔に書くことを気にしすぎなくても良いが。))。
for( el in x ){ #
if( el%%3 == 0 ) next # el が 3 の倍数だったら、反復の先頭へ戻れ。
if( el > 7 ) break # el が 7 より大きかったら、反復を終了せよ。
cat(" ",el,"\n") # el を表示せよ。
}
印刷制御文字は C言語と同じで、たとえば、"\n"で改行が行われる。
もう一つ関数の例をあげておく。
get.prime<-
function( n=20 ){ # n 番目までの素数を求める。
p <- 1
prime.list <- 2
while(length(prime.list) < n){ # 素数リストの長さが n 未満な限り以下を実行
p<-p+2 # 新規の素数候補を p とする。
if( all(p%%prime.list != 0 ) ) # 既存の素数での剰余が 全て 0 でないなら、
prime.list<-c(prime.list,p) # 素数リストに p を追加。
}
prime.list # 素数リストを返す。
}
以上詳しくは、マニュアル[1]を参照されたい。
5. グラフィック環境
グラフィックウィンドーを開くには以下をする。
> motif()そのあと、たとえば、
> hist( rnorm(1000) )とすれば、平均 0, 標準偏差 1 の正規乱数 1000個のヒストグラムをそのウィンドーに 表示してくれる。
breaksというオプションで階級の区分点を指定できる。
> hist( rbinom(200,10,0.4), breaks=-0.5:10.5 )表題や日付を図に付けたいときには、以下のようにすればよい。
> title( "200個の二項乱数 b(10,0.4)" ) > stamp()
関数 plot()は、同じ長さの引数を2つ取って、
座標を適当に設定してから、散布図を書いてくれる。
> x<-rnorm(500); y<-x+rnorm(500) > plot(x,y)表記する範囲を指定したければ、次のような指定を追加すれば良い。
> plot(x,y,xlim=c(-4,4),ylim=c(-6,6))
既存の座標上に線を追加したければ、
> lines(c(min(x),max(x)),2*c(min(x),max(x)))などの関数もある。 普通の線グラフも関数
plot()で書ける。
> x<-(-300:300)*0.01 > plot(x,dnorm(x),ylim=c(0,1),type="l") > lines(x,dnorm(x,sd=0.5))引数
typeはプロットの形式を定めるもので、
"p","l","b","n","h" などで、それそれ、
点プロット(省略時)、
線プロット、
点と線のプロット、
軸のみでプロットしない、
各点からの垂線プロット
を意味する。
その他のグラフの表示形式
(軸の引き方、軸名の付け方、表題の付け方などなど)
については、関数 par()をマニュアルで参照されたい。
motif()で描画した結果は、プリントアウトすることができる。
グラフィックウインドウで Graph のボタンを押し、
さらに開いたメニューで printのボタンを押すと、
描いた結果がpostscript形式に変換されプリンタに送られる。
( ただし、前提として、motif ウインドウ 上方にある Options の Printing での選択画面の
"Command"の指定場所にその計算機環境での出力用のコマンド["lpr -Pprt1"など]が
入力されている必要がある。)
また、このときのpostscript形式のファイルが
ps.out.0001.psなどの名前で、
Splusを開始したディレクトリーに作成されている。
単に、postscript形式のファイルだけが欲しい
(紙へのプリントアウトはいらない) 時には、
グラフィック画面 オプション(Options)のメニュー中の印刷方向(printing...)で、
プリントアウトを指定するコマンドを、
存在しないコマンド(たとえば"aaa")に変更しておいて、
グラフィック画面上のボタンでプリントアウトすればよい
(コマンドがないというエラーメッセージが出ますが、気にしなければ良い)。
> postscript(file="ファイル名") plot など作図のためのコマンドの実行 >dev.off()をすればよい。
dev.off()で作図機器(postscript)を閉じたときにファイル出力が実行される。
途中経過は motif の画面には出力されない(postscriptファイルに出力される)ので、
実施する plot など作図のためのコマンドで
どんな画面が作成されるのか、事前に motif の画面で描画して確認してから実行するとよい。
\epsfileなど用いれば良い。
(例)
6. 最後に(受講学生へ)
やって見てわからないこと、どうもシステムがマニュアルどおり動いていない?、
このメモの誤りを見付けたなど、なんでも関まで問い合わせて下さい。
わかる事はお答えしますし、不具合は善処します。
メイルでも結構です。
注意:
Splus を実行すると、~/MySwork/ なるディレクトリが作られる。
さらに、そこに .Data なるディレクトリが作成される。
.Data はピリオドで始まる名前なので、UNIX環境で、
ただ ls としても表示されない。以下で表示される。
% ls -al ~/MySwork
ディレクトリ ~/MySwork/.Data
で、 ls -al とすると、 .Audit なるファイル名がある。
このファイルは、Splus の実行履歴を記録しているので、
これを見れば、過去に自分が実行したコマンドを知ることができる。
ただし、上限までは実行しただけ大きくなる。上限の文字数は以下で知ることができる。
> options()$audit.size
[1] 5e+05 # デフォルトの大きさ( 5MB )
このまま放っておくと、他にも大きなファイルがある場合、
個人のディスク制限 20MB をオーバーしてしまい、
ログインできなることもある。以下を実行しておくことを勧める。
> options(audit.size=10000) # 10000文字分だけ取っておく
また、Splus 内で与えたオブジェクトも ~/MySwork/.Data 以下にファ
イルとして残る。特に、実験で一億個の乱数を生成して臨時のオブジェクトに付置した
xx <- rnorm(100000000)
というようなときに、xx が以後、必要なければ Splus環境から rm(xx)
などで消すこと。
消さないといつまでも残って、利用可能ディスク容量を減らすことになる。
ディスク制限(quota)は、ログイン時に表示されるので、注意すること。
また、以下でもわかるので注意すること。
% /usr/ucb/quota -v ログイン名
- Splus のプロセスについて
Splus を終了するときは、必ず、以下で、終了すること。
> q()
上で終了しなかったときには、必ず、
% ps -f
で Splus (/opt/splus6/cmd/Sqpe) がプロセスに残っていないことを
確認しなさい。
もし残っていたら、ps -f で見たプロセス番号(PID)を利用して、
% kill -9 プロセス番号
として、プロセスを残さないようにすること。