データ科学(2008年度後期 開講)で作ったスクリプト
()内は使ったデータファイル名
スクリプト内で, データの読み込みを行っている部分(基本的に一行目)は保存場所を指定してください
データ読み込みからボックスプロット (HOUSE.txt)
HIVと車の台数の回帰～見かけの相関がある回帰分析～ (HIV_Car.txt)
住宅価格の値段と面積の単回帰分析 (HOUSE.txt)
住宅価格の値段に対する各データでの重回帰分析 (HOUSE.txt)
新生児の体重に対する重回帰分析～stepAICを用いた変数選択法～(作りかけ) (NeoWeight.txt)
菖蒲のデータに対する主成分分析とbiplot～多次元データを２次元にプロット～ (IrisData.txt)
DBHのデータに対する主成分分析とbiplot～多次元データを２次元にプロット～ (DBHData.txt)
住宅価格のデータに対する主成分分析とbiplot～多次元データを２次元にプロット～ (HOUSE.txt)
菖蒲のデータに対するクラスター分析～kmeans法を用いた手法～ (IrisData.txt)
DBHのデータに対するクラスター分析～kmeans法を用いた手法～ (DBHData.txt)
USACrimeDataに対するクラスター分析～階層的クラスタリング～ (USACrimeData.txt)
三国志の武将のデータに対するクラスター分析～階層的クラスタリング～ (Sangokushi.txt)
ICUデータに対する一般化線形回帰～ロジスティック回帰～ (ICUData.txt)
カブトムシのデータに対する一般化線形回帰～ポアソン回帰～ (BeetleData.txt)
金のデータに対する回帰(作りかけ) (Gold.txt)
Rコマンダーについて
Rの命令は, すべてコマンド入力…
→コマンドをある程度覚える必要.
→GUI(マウスによる操作；Graphical User Interface)で何とかできないものか…
→Rコマンダーというパッケージ！！
参考になりそうなPDF
注意!バージョンが古いので, 若干設定が違う…(と思う)
→16ページ目からはRコマンダーの使用法なので参考になる！！
Rコマンダーのインストール
Rコマンダーのパッケージをインストールした後, Rを起動
library(Rcmdr) と入力
「パッケージが足りない」のような警告が出る
“はい(Y)”を選択
足りないパッケージのインストール開始
時間がかなりかかります
終了後、Rコマンダーが起動します。 毎回やるの？ 一度、インストールすると、あとはRを起動し, library(Rcmdr) と入力すればよい
PDFの13ページのウィンドウが出てくる
→使い方は, そのPDFの16ページ目からを参照するとよい
使ったデータ
HOUSE.txt
HIV_Car.txt
NeoWeight.txt
IrisData.txt
DBHData.txt
USACrimeData.txt
Sangokushi.txt
ICUData,txt
BeetleData.txt
Gold.txt
メール
(Matlab上からのメールの送り方→コチラ)
Matlabで出来るんだから、Rでもなんとか出来ないものか・・・と調べた結果、出来そうです。
"sendmailR"というパッケージを取ってくる (たとえばCＴAN)
使うメールサーバのSMTPを調べる
sendmail(from="送り元メールアドレス",to="送り先メールアドレス",subject="件名",msg="本文",control=list(smtpServer="SMTPサーバー"))
が基本構文
・送り元メールアドレス；送り元として設定するメールアドレス
・送り先メールアドレス；送り先のメールアドレス
・件名；送るメールの件名
・本文；メールの本文
・SMTPサーバ；SMTPサーバのアドレス
補足；そんなに細かい設定が不要ならば、"mail"パッケージ (CRAN) を使うことも可能
資料
sendmailRの説明 (英語)
mailの説明 (英語)
ベクトルの最小値・最大値のある場所を探したい (最小値や最大値を返す関数はある)
(Matlabだと, min・max関数の戻り値として得られる)
あるベクトルAの最小値・最大値を取っている要素の番号を探したい.
式で書くなら, argmin(A), argmax(A)とかで書かれる話.
(二か所で最小値・最大値を取る場合は, 最初の要素番号を見つける)
最小値について
if(sum(min(A)==A)==1){
argmin<- sum(1:length(A)*(min(A)==A))
}else{
argmin<-min(1:length(A)*(min(A)==A)+(length(A)+1)*((1:length(A)*(min(A)==A))==0))
}
・sum(min(A)==A)==1は, 最小値が一つかどうかを調べている.
・ sum(1:length(A)*(min(A)==A))は, (min(A)==A) で最小になる要素の部分をTRUE(1), それ以外をFALSE(0)のベクトルにする.
それと 1:length(A) を掛けて, 最小になる要素の部分だけその要素番号で他が0のベクトルを作る.
そして, その和を取れば, 最小になる要素の要素番号が得られる.
・else以下は, 最小値を取る要素が二つ以上ある場合を調べている.
・1:length(A)*(min(A)==A)は前と同じで, 最小になる要素のところに要素番号を入れる.
・(1:length(A)*(min(A)==A))==0)は, 最小値じゃないところをTRUE(1), 最小値をFALSE(0)にしている.
そこに, ベクトルの長さ (length(A))+1を書けることで, 最小じゃないところは要素数+1の番号となる.
その和を取ることで, 最小値の部分は要素番号, 最小値以外の部分は要素数+1のベクトルになるので, その最小値をとればよい.
（もしかしたら, if文での分岐は不要で, else以下のみでも行ける気がする・・・）
最大値について
最小値と同様に考えていく
if(sum(A==max(A))==1){
argmax<-sum(1:length(A)*(max(A)==A))
}else{
argmax<-min(1:length(A)*(max(A)==A)+(max(A,length(A))+1)*((1:length(A)*(max(A)==A))==0))
}
・sum(A==max(A))==1；最大値となる要素が一つかどうか
・sum(1:length(A)*(max(A)==A))は, (max(A)==A) で最大になる要素の部分だけTRUE(1), それ以外がFALSE(0)のベクトルにする.
それと 1:length(A) を掛けて, 最大になる要素の部分だけその要素番号で他が0のベクトルを作る.
そして, その和を取れば, 最大になる要素の要素番号を得る.
・else以下は, 最大値を取る要素が二つ以上ある場合を調べている.
・1:length(A)*(max(A)==A)は, 最大値の要素に対応する部分はその要素番号が残って, 他が0のベクトルを作る
・((1:length(A)*(max(A)==A))==0)で最大値じゃないところがTRUE(1), 最大値のところがFALSE(0)になるベクトルを作る
そこに, max(最大値,要素数)+1の値を掛ける (最大値じゃない部分が最大となる要素番号を必ず超えるようにするため)
この和を取ることで, 最大値に対応する部分は要素番号となり, 最大値以外の部分は max(最大値,要素数) の値のベクトルになる
・その最小値を取ることで, 最初の最大値の要素番号を返す
（もしかしたら, if文での分岐は不要で, else以下のみでも行ける気がする・・・）
プロキシを使った接続をしている場合
使っているネットワークによっては, インターネットに接続する場合にプロキシを通した接続をしていることがある (特に学校関係). その場合, Rにもプロキシの設定を"教えなければならない". これが分からない場合は, 次のようなエラーが出る
警告： リポジトリー http://cran.md.tsukuba.ac.jp/bin/windows/contrib/2.15 に対する索引にアクセスできません
警告： リポジトリー http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/2.15 に対する索引にアクセスできません
警告メッセージ：
In open.connection(con, "r") :
'cran.r-project.org' をポート 80 でコネクトできません
(日本語以外の言語を用いている場合は対応している言語での警告が出る)
これは, "アクセスしようと思ったけど出来なかった"という警告メッセージで, これはプロキシの設定をRに教えていないことが原因である.
そこで, Windowsユーザでショートカットを使っている場合は, ショートカットアイコンを右クリックし, プロパティを選択する. そこのショートカットタブの中の「リンク先(T):」という項目にショートカット先が書かれているので, そこのボックスの最後に「 --internet2」を追加する. すると, 「"～～～.exe" --internet2」という形になるので, OKを押して, (管理者権限じゃなければ管理者権限でOKですか?が出てくるので, そのままOK押して) 終了.
もう一回立ち上げなおすと, プロキシの設定を教えた状態での接続になるっぽい. 詳しくは, 参考URLを見てください.
参考URL; 脳のなかのこびと軍団 Rのプロキシ設定の(2)
アニメーションの作成 2015/6/17-修正 (一部変数名のミスで分散の設定が上手くできなかったため)
animation libraryを使うとアニメーションが作れる
相関係数を色々変えて見たアニメーションを描画するスクリプト
散布図を描画・回帰直線を描画 を止めつつアニメーションにして、軸を固定したアニメーションのためのスクリプト２、３
スクリプトのコメントで色々書いたので、詳細はそちらを
使うものは
ani.options関数、ani.record関数、ani.replay関数
ani.options関数では、アニメーション描画間隔の変更で使用（不要）
ani.record関数により、「前にplotしたものを保存する」
注意；ani.record関数の後にplotしたものを保存するのではなく、ani.record関数の前にplotしたものが保存される
ani.record()→plot だと「ここのplotで描いた図」がアニメーションに保存されない
plot→ani.record() だと「ここのplotで描いた図」がアニメーションに保存される
これが分からず、アニメーションが変な順で保存されてた
ani.replay関数でアニメーションを描画
このスクリプトを作るときに関して参考にしたＨＰ
plot関係；とりあえずplot()、グラフのいろいろな設定と重ね描き、低水準作図関数、
アニメーション関係；これからの可視化は動画の時代～Rでanimationパッケージで動画を作成する方法、その内部のスライド
回帰直線だけではなく楕円を描写させる
とりあえずで作ったものはこちら
ざっくり作ったので、よりよくするためには修正必須
授業で見せる分にはこれで十分・・・？
とりあえずのものなので、説明変数の分散や反応変数の分散を大きくすると上手く描写されない。
恐らく、10程度が限界？
楕円を描く際に、この辺も注意すればもっとよいものになりそう。
やってることは、上のものに加え、楕円の描写を入れている
楕円を描く際に、データの主成分分析を行い、第一主成分と第二主成分の軸に対する標準偏差を計算し、それを長軸・短軸とした。
また、楕円の中心は（説明変数の平均, 反応変数の平均）にした
もうちょっとやりようがあると思う
参考にしたＨＰ
主成分分析について；「Ｒと主成分分析」http://www1.doshisha.ac.jp/~mjin/R/24/24.html
楕円について；「■楕円」http://www.geisya.or.jp/~mwm48961/kou3/quadratic_1.htm
Rでの楕円のプロット；「図形描写関数群」http://aoki2.si.gunma-u.ac.jp/R/plot.html
ネットに繋がっていれば、ここのソースから自動で引っ張ってくれる
ネットに繋がっていない場合、ソースの部分だけコピペして保存→絶対パスで読み出す
ネットに繋がっていない場合に、ソースの細かいところを弄るのが面倒なので、そこをスクリプトの頭のほうでやるために、文字列の操作を使っている
参考にしたＨＰ；「17.文字列を操作する」http://cse.naro.affrc.go.jp/takezawa/r-tips/r/17.html
楕円のプロットの関数についての解説；「2 R で図を書いてみよう」http://www.htc.nagoya-u.ac.jp/~yamamoto/lecture/11/mc_cu11_R110930.pdf
図形の回転をさせるための角度を求めるときに使用；「03.簡単な計算」http://cse.naro.affrc.go.jp/takezawa/r-tips/r/03.html （この中のatanを使って回転角度を作っている）
気付いたこと
for文で回しているときに, sprintfを使った場合は表示されない(?)
→cat関数を使う；cat("---",a,"\n")とやると, --- (aの値) と出力され改行("\n"を入れているので).
現在時刻は, Sys.time()関数で取れる.
処理時間の測定は, 最初にproc.time()を何かに入れて一旦現在の(PC内の)時刻を保存しておいて, 処理が終わってからproc.time()と最初に入れたものの差を取ればよい. (単位は秒だったはず)
“paste”を使えば, 簡単に色んな文字列の結合が可能. 途中に変数を文字に変換して結合することも可能.
たとえば, s<-paste(FN,"p=",p,", k=",k,".csv") とすると, sは「FNの値 p=pの値 , k=kの値 .csv」となる.
色々組み合わせることで, 設定を反映した文字列を作って, それをファイル名として使って結果を書き出すことが可能.
なぜか, "--"の前に半角スペースが入るので, 注意が必要.
参考ページ