header

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS

データ解析 担当:間瀬茂 学部(情報科学科, (2013) 5学期) 前期(火曜日3,4限)

講義の狙い

データの統計的解析を,フリーの統計解析システムRを使って実習する。最終的にいくつかの実データを用いたデータ解析の実習を行う.

内容

  • 統計解析システムRの紹介とインストール
  • Rの使用法の基礎
  • Rを用いた最適化
  • Rを用いた回帰分析と変数選択
  • Rを用いた最尤推定
  • Rを用いたシミュレーション
  • 実データを用いた実習

テキスト

テキスト ,及び用意した資料を使用する。

参考情報

評価

出席及びRを用いたデータ解析実習のレポート.

注意

実習は各自のパソコンにインストールしたRを用いて行う.必要な資料は履修者へ メイルで指示するので,受講希望者は教員宛てに受講希望のメイルを予め送ること. またこのページの末尾に受講者への連絡や添付ファイルを順次掲載するので,注意してほしい.

データ解析の講義で使用したRコードやデータは データ解析の資料URL に順次アップロードしますので必要に応じて参照してください.

受講者への連絡

特別講演のお知らせ.新任の助教野村俊一さんは損保会社に勤務し,アクチュアリストの資格を持つという異色の経歴を持っています.野村さんにアクチュアリの試験や業務について個人的経験から話して頂く講演会が以下のように開かれます.関心のある方は自由に参加してください.

時間: 2013年7月18日(木) 13:20-14:50
場所: 大岡山西8号館W棟 W809号室

講演者: 野村 俊一 氏 (東京工業大学 数理・計算科学専攻・間瀬研究室 助教)
題目: アクチュアリーの業務、試験、課題について

概要: 「アクチュアリー」とは,保険・年金の業界で活躍する数理業務のスペシャリストである.
本発表では、自身の経験を踏まえてアクチュアリーに関する以下の内容を紹介する.
・アクチュアリーの役割と前職(損害保険会社)での経験
・アクチュアリー資格試験を通るまでの軌跡と勉強方法
・アクチュアリーが現在取り組んでいる主要な課題

(2012.07.16, データ解析第四回レポート課題)

提出期限はこれまでのを含め8月16日(金)厳守! 提出はこれまでのようにメイルを利用してください. 人に相談は構いませんが,丸写しは厳禁.非常識なほど似たレポートがあればすべてを採点しません. レポートには「データ解析第4回レポート」と明記,名前と学籍番号も必ず明記のこと. 複数のレポートを一緒にまとめて提出しないこと.

第四回レポートの課題:誕生日が一年365日,一様かつ独立なn=100人をシミュレーションし,その中の最大誕生日一致人数を返す関数BD2を定義せよ.BD2をN=10000回実行し,その結果を表にまとめよ.

ヒント: 100人の誕生日を無作為に選ぶには

x <- sample(1:365,100,replace=TRUE)  # 復元抽出

とする. 100人の誕生日の一致状況を判定するアルゴリズムは,素朴にコードすると複雑且つ時間がかかりがちである.次の例を参考にすること.これ以外の実行可能なアルゴリズムがあれば考えてください:

> x <- c(1,1,2,7,4,5,8,5,5)
> (y <- sort(x))
[1] 1 1 2 4 5 5 5 7 8
> z <- rle(y) # 関数 rle は数列 y の引き続く同じ値の数を集計する
> z
 Run Length Encoding
  lengths: int [1:6] 2 1 1 3 1 1  # つまりベクトルyには1が2回,2が1回,,,5が3回,,,8が1回と続く
  values : num [1:6] 1 2 4 5 7 8
> z$lengths
[1] 2 1 1 3 1 1

(2013.07.16 講義資料)

今週の講義のデモ用ファイルを ここ に置きました.L11.R です.

(2012.07.09, データ解析第三回レポート課題)

提出期限は二週間.(人に相談は構いませんが,丸写しは厳禁).提出は私宛のメイルの件名に「データ解析第三回レポート」と明記して,添付ファイルもしくはメイル本文にテキストとして直接書いてください.名前と学籍番号も必ず明記してください.

BMI指数は体重と身長から簡単な計算で得られることが,広く普及した理由の一つです.しかし簡単な計算で得られるものが有用であるという根拠はありません.今回の課題では,体脂肪率の予測という観点から,新しい BMI指数 nBMI=(体重Kg)^a/(身長m)^b を定義し,最適なパラメータ a,b, を人体計測値データから決定してみてください.

注意:人体計測データファイル bodysize.csv はここからダウンロードして使ってください.

X <- read.csv("bodysize.csv") # 人体の各種測定値データフレームを読み込む
Y <- X[,2:4]                      # 身長,体重,体脂肪率変数だけのデータフレーム
names(Y)[1:3] <- c("height","weight","fat") # 身長,体重,体脂肪率変数名を簡略化
Y[[1]] <- Y[[1]]/100              # 身長変数の単位をメートルに変換
attach(Y)                         # データフレーム Y 内の変数を直接アクセスできるように登録
BMI <- weight/(height)^2          # BMI指数のベクトルを計算
fm <- nls(fat ~ weight^a/height^b, data=Y, start=c(a=1,b=2)) # a,bの初期値は古典的BMI指数のそれを使用
###############################################################
# 以下を検討してください:
# (1) 非線形回帰結果 fm から,新しいBMI指数データ nBMI を求める (ヒント:fitted 関数)
# (2) パラメータ a,b の値を求める (ヒント:coef関数)
# (3) 相関係数 cor(fat,BMI), cor(fat,nBMI) を比較する
# (4) 新旧のBMI指数値とfat値の差の標本標準偏差 sd(fat-BMI),sd(fat-nBMI)を比較する

(2012.07.9 講義資料)

今週の講義のデモ用ファイルを ここ に置きました.L10.R です.

(2012.06.25 講義資料)

今週の講義のデモ用ファイルを ここ に置きました.L9.R です.

(2013.06.18 講義資料)

今週の講義のデモ用ファイルを ここ に置きました.L8.R です.

(2013.06.11 講義資料)

今週の講義のデモ用ファイルを ここ に置きました.L7.R です.データファイル sname.RData, Nelson.RData も必要です. sname.RData, Nelson.RData はRオブジェクトをそのまま記録したファイルで, load命令で読み込みます:

load("sname.RData") # オブジェクト sname が復元される

load("Nelson.RData") # オブジェクト Nelson が復元できる

(2012.06.11, データ解析第二回レポート課題)

提出期限は今回も二週間です(人に相談は構いませんが,丸写しは厳禁).提出は私宛のメイルの件名に「データ解析第二回レポート」と明記して,添付ファイルもしくはメイル本文にテキスト(こちらがベター)として直接書いてください.名前と学籍番号も必ず明記してください.長すぎる場合は適宜編集してください.

BMIは相対的な肥満度を表す簡便な指数として広く使われています.その根拠の一つとして体脂肪率と相関が高いことが挙げられています.しかし,本当でしょうか.今回のレポートの内容はBMIと体脂肪率の関係,及び体脂肪率を他の身体計測値を用いて表す(重回帰)です.以下の作業を実行するRのプログラムと最終的に得られた結果(summary関数の出力)です.まず次の作業を行なってください.各段階でのRの出力は含めても含めなくても結構です.

X <- read.csv("bodysize.csv") # 人体の各種測定値データフレームを読み込む(Windows OSの場合)

str(X)

Y <- X[,-1] # ID変数を取り除いたデータフレーム

str(Y)

names(Y)[1:3] <- c("height","weight","fat") # 身長,体重,体脂肪率変数名を簡略化

### BMI指数を計算

attach(Y) # データフレーム Y 内の変数を直接アクセスできるように登録

BMI <- weight/(height/100)^2 # BMI指数のベクトルを計算

cor(fat,BMI) # BMI と体脂肪率の相関係数

### 体脂肪率の箱型図とはずれ値

boxplot(fat, BMI) # 体脂肪率とBMIの並行箱型図を描く

plot(fat,BMI) # 体脂肪率とBMI指数の散布図を描く

names(BMI) <- as.character(1:813) # BMI指数ベクトルに通し番号ラベルを付ける

boxplot.stats(BMI)$out # BMI指数の外れ値を確認

### 体脂肪率をBMI指数に回帰

summary(y1 <- lm(fat ~ 1+BMI)) # 体脂肪率をBMI指数に回帰(定数項を含む)し要約

op <- par(mfrow=c(2,2)) # 描画領域を2x2に分割

plot(y1) # 4つの回帰診断図を一度に表示

matplot(fat,cbind(predict(y1),resid(y1))) # 体脂肪率と予測値,残差の同時プロット

abline(0,0) # x軸を重ね書き

### 体脂肪率を他のすべての変数に重回帰

summary(y2 <- lm(fat ~ ., data=XX)); plot(y2)

引き続き,step 関数で変数選択した結果,fat 変数を他のすべての変数に二次回帰した結果,そしてそれを step 関数で変数選択した結果を吟味してください.swiss データの解析プログラムが参考になるでしょう.考察も添えてください(統計的データ解析では,解析結果の吟味こそが本当の作業です).

(2012.06.04 講義資料)

今週の講義のデモ用ファイルを ここ に置きました.L5.Rです.

(2012.05.28 講義資料)

今週の講義のデモ用ファイルを ここ に置きました.L6.R です

(2012.05.21 講義資料)

今日の講義のデモ用ファイルを ここ に置きました.swiss.R です.

(2012.05.14, データ解析第一回レポート課題)

提出期限は二週間です.それ以降も受け付けますが遅延ペナルティーがつくと思ってください.(人に相談は構いませんが,丸写しは厳禁).提出は私宛のメイル maseshigeru+da@gmail.com の件名に「データ解析第一回レポート」と明記して,添付ファイルもしくはメイル本文にテキストとして直接書いてください.名前と学籍番号も必ず明記してください.

レポートの内容は以下の作業を実行するRのプログラムと最終的に得られた集計表です.各段階でのRの出力は含めても含めなくても結構です.長すぎる場合は適宜省略してください.最終結果の集計表は付け加えてください.

======================================

WHOの基準によれば,BMI(Body Mass Index)

BMI=体重(Kg)/(身長(m))^2  (身長はメートル単位に注意!)

は次の8段階に分類されます:

1  16未満         重度の痩身
2  16以上17未満   中度の痩身
3  17以上18.5未満 軽度の痩身
4  18.5以上25未満 普通
5  25以上30未満   肥満
6  30以上35未満   肥満クラス1
7  35以上40未満   肥満クラス2
8  40以上         肥満クラス3 

(1) Rで体重 w(Kg) と身長 h (cm) を引数とし,WHO基準クラス番号 1,2,...,8 を返す関数 BMI(w, h) を定義する.

(2) データ physical をデータ解析の資料URLからダウンロードし Rで読み込む(日本人の身体計測値データ)

physical <- read.csv("bodysize_s.csv")  # Windows, Mac なら
physical <- read.csv("bodysize_u.csv")  # Linux なら

(3) str(physical) でオブジェクトの構造を確認.特に総データ数と,身長と体重変数が何番目か確認.

(4) BMI関数(と必要なら for ループ,if 文)を使い全員の肥満度クラス番号を計算し,ベクトル bmi に入れる.(データフレーム physical から体重と身長変数を抜き出すには二重鈎括弧 physical[[.]] を使う)

(5) bmi を集計表にまとめる(つまり、クラス1,2,3,4,5,6,7,8がそれぞれ何人ずついるか。ヒント;関数 table を使用)

(2012.05.14 講義資料)

今週の講義のデモ用ファイルをアップロードしました.L4.R です

(2012.05.7 講義資料)

今日の講義のデモ用ファイルを データ解析の資料URL に置きました.L3.R です.

(2012.4.23)

今週の講義で使ったデモ用のRのソースファイル L2.R を データ解析の資料URL に置きました. 漢字コードは MS Windows 用の Shift JIS です.Linux や Mac OS を使用している人は適当なソフトで漢字コードを使用OS用に変換してください.

(2012.4.16)

今日の講義で使ったデモ用のプログラムのソース L1.R をデータ解析の資料URLに置きました.またデモ用の関数のソースファイルが Demo.R です.講義と同じ内容を実行するには,まず必要なパッケージ evaluate (とそれが依存する二つのパッケージ)を使っているパソコンにインストールする必要があります.これは実行中のR端末から命令

install.packages("evaluate", dependencies=TRUE)

を実行するとひとりでにインストールできます.途中で質問がありますが,CRANのミラーサイトは日本のどれかを選んでください.ファイル Demo.R を R の作業用フォルダー(ディレクトリ)においてから,

source("Demo.R")

で Demo 関数を読み込みます.ls()命令で確認してください.最後に

Demo2("L1.R") 

とするとデモが始まります.CRキーを押すたび,命令と対応実行結果が順に表示されます なお,Demo.R,L1.R を始めとして日本語を含むファイルは使っているOS毎に日本語のコーディングが異なります.利用者が最も多いと思われる Microsoft Windows 用のコーディングである Shift JIS 用のファイルだけを提供します.Linux 利用者は漢字コーディング変換用のソフト nkf を利用して例えば

nkf -e --overwrite Demo.R  

という命令で漢字コーディングを EUC に変更してください.

(2013.4.16)

CRAN (Comprehensive R Archive Network)の RjpWiki から,例えばWindows用のバイナリーが手に入ります.日本のRユーザーの情報交換サイト RjpWiki には各種OSへのRのインストール方法の紹介がありますが, 賞味期限切れの記事も結構ありますので注意.私の愛用している Ubuntu Linux ではRのインストールは

sudo apt-get install r-base

主要パッケージの一括インストールは

sudo apt-get install r-cran*

です.

Counter: 5268, today: 1, yesterday: 0