■
test
Rでコンジョイント分析
Rにはコンジョイント分析用のパッケージが無いなー、とずいぶん前から思っていましたがいつの間にやら登録されていたようです。
コンジョイント分析とは何か？についてはアルベルト社の解説が分かりやすいと思います。
install.packages("conjoint") #パッケージのインストール。初回のみ。 library(conjoint) #パッケージの呼び出し。 data(tea) #サンプルデータの呼び出し # 各プロファイルの表示 # 13プロファイル（4属性（3水準、3水準、3水準、2水準）） print(tprof) # 各プロファイルの表示 # price variety kind aroma # 1 3 1 1 1 # 2 1 2 1 1 # 3 2 2 2 1 # 4 2 1 3 1 # 5 3 3 3 1 # 6 2 1 1 2 # 7 3 2 1 2 # 8 2 3 1 2 # 9 3 1 2 2 # 10 1 3 2 2 # 11 1 1 3 2 # 12 2 2 3 2 # 13 3 2 3 2 # 各水準のラベル print(tlevn) # levels # 1 low # 2 medium # 3 high # 4 black # 5 green # 6 red # 7 bags # 8 granulated # 9 leafy # 10 yes # 11 no # 選好マトリクス。数値が大きい方が高評価。 head(tprefm) #100名×13プロファイル # profil1 profil2 profil3 profil4 profil5 profil6 profil7 profil8 profil9 profil10 profil11 profil12 profil13 # 1 8 1 1 3 9 2 7 2 2 2 2 3 4 # 2 0 10 3 5 1 4 8 6 2 9 7 5 2 # 3 4 10 3 5 4 1 2 0 0 1 8 9 7 # 4 6 7 4 9 6 3 7 4 8 5 2 10 9 # 5 5 1 7 8 6 10 7 10 6 6 6 10 7 # 6 10 1 1 5 1 0 0 0 0 0 0 1 1 # 分析する。最小二乗法。 # 引数は順に選好マトリクス、プロファイル、水準のラベル Conjoint(tprefm, tprof, tlevn) # [1] "Part worths (utilities) of levels (model parameters for whole sample):" # levnms utls # 1 intercept 3,5534 # 2 low 0,2402 # 3 medium -0,1431 # 4 high -0,0971 # 5 black 0,6149 # 6 green 0,0349 # 7 red -0,6498 # 8 bags 0,1369 # 9 granulated -0,8898 # 10 leafy 0,7529 # 11 yes 0,4108 # 12 no -0,4108 # [1] "Average importance of factors (attributes):" # [1] 24,76 32,22 27,15 15,88 # [1] Sum of average importance: 100,01 # [1] "Chart of average factors importance"
各水準の部分効用値と各属性の寄与率のグラフも同時に出力される。
price
variety
kind
aroma
寄与率
また、結果の分析だけでなくプロファイルの作成もできる。
experiment <- expand.grid( price = c("low", "medium", "high"), variety = c("black", "green", "red"), kind = c("bags", "granulated", "leafy"), aroma = c("yes", "no")) design <- caFactorialDesign(data=experiment, type="orthogonal") # フルプロファイル法はtype="full" # 直交計画はtype="orthogonal" # を指定した計画type="fractional", cards=16 design # price variety kind aroma # 3 high black bags yes # 5 medium green bags yes # 10 low black granulated yes # 17 medium red granulated yes # 22 low green leafy yes # 27 high red leafy yes # 34 low red bags no # 42 high green granulated no # 47 medium black leafy no # 分析に使うときはこれをさらにコーディングする。 caEncodedDesign(design) # price variety kind aroma # 3 3 1 1 1 # 5 2 2 1 1 # 10 1 1 2 1 # 17 2 3 2 1 # 22 1 2 3 1 # 27 3 3 3 1 # 34 1 3 1 2 # 42 3 2 2 2 # 47 2 1 3 2
ちなみにできたプロファイルの直交性を確認したい時は、プロファイルをコーディングした後、相関係数を見ればよい。
cor(caEncodedDesign(design)) # price variety kind aroma # price 1 0 0 0 # variety 0 1 0 0 # kind 0 0 1 0 # aroma 0 0 0 1
他に
- 回答者ごとの各水準の部分効用値: caPartUtilities(y, x, z)
- 回答者ごとの各プロファイルの全体効用値: caTotalUtilities(y, x)
- 上記全体効用値に基づいたクラスタリング(k-means): caSegmentation(y, x, c)
- コンジョイントシミュレーション手法
- 最大効用: caMaxUtility(sym, y, x)
- Bradley-Terry-Luce(BTL): caBTL(sym, y, x)
- ロジット: caLogit(sym, y, x)
といった関数が実装されています。
評価の無いプロファイルがある場合に対応してないとか、最小二乗法以外の解法が無いとか物足りない点もありますが、けっこう使い勝手が良さそうじゃないですか？これ。
なお、選択型コンジョイント分析(CBC)についてはデータ解析環境Rによる選択型コンジョイント分析入門(PDF注意)を参照されたい。
数独をRで解く
http://apollon.issp.u-tokyo.ac.jp/~watanabe/sample/sudoku/index_j.htmlに載っている数独の問題をRで解きます。
といっても、例によって必要なパッケージを呼び出して、関数一発です。
そう、数独用のパッケージがあるのです、Rならね。
install.packages("sudoku") library(sudoku) HW1 <- matrix(c(0,6,1,0,0,7,0,0,3, 0,9,2,0,0,3,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,8,5,3,0,0,0,0, 0,0,0,0,0,0,5,0,4, 5,0,0,0,0,8,0,0,0, 0,4,0,0,0,0,0,0,1, 0,0,0,1,6,0,8,0,0, 6,0,0,0,0,0,0,0,0), 9, byrow=TRUE) solveSudoku(HW1, verbose=TRUE) Inkara1 <- matrix(c(0,0,5,3,0,0,0,0,0, 8,0,0,0,0,0,0,2,0, 0,7,0,0,1,0,5,0,0, 4,0,0,0,0,5,3,0,0, 0,1,0,0,7,0,0,0,6, 0,0,3,2,0,0,0,8,0, 0,6,0,5,0,0,0,0,9, 0,0,4,0,0,0,0,3,0, 0,0,0,0,0,9,7,0,0), 9, byrow=TRUE) solveSudoku(Inkara1, verbose=TRUE) Inkara2 <- matrix(c(8,0,0,0,0,0,0,0,0, 0,0,3,6,0,0,0,0,0, 0,7,0,0,9,0,2,0,0, 0,5,0,0,0,7,0,0,0, 0,0,0,0,4,5,7,0,0, 0,0,0,1,0,0,0,3,0, 0,0,1,0,0,0,0,6,8, 0,0,8,5,0,0,0,1,0, 0,9,0,0,0,0,4,0,0), 9, byrow=TRUE) solveSudoku(Inkara2, verbose=TRUE)
結果はネタバレになるので書きません。
処理時間はHW1が0.39秒、Inkara1が3.12秒、Inkara2が10.95秒でした。
第29回Ｒ勉強会＠東京( #TokyoR )の予習メモ
今週末の Tokyo.R #29 ( http://atnd.org/events/36417 ) はいつもにも増して興味深い発表が多いので、より理解できるようにちょっとだけ予習というかメモ書きをしておく。
予習するのは以下の３本。
- RでAHP ( [twitter:@Soultoru] )
- シリーズ前処理2013：外れ値（仮）( [twitter:@sfchaos] )
- Rによる独立成分分析( [twitter:@mikiya0417] )
★RでAHP
AHP( Analytic Hierarchy Process) は一対比較法の一種。
一対比較法は官能評価などで３個以上の対象物を２個の組み合わせの相対評価から、心理的、主観的な評価を数量化（尺度化）する方法。
いくつかの手法があるが、よく知られたものとしてサーストン法（プロビットモデル）、シェッフェ法、ブラッドレイ・テリー法（ロジットモデル）などがある。
AHPはこれを階層化したものでオペレーションズ・リサーチの世界で発達してきた。
RでAHPを行うには群馬大の青木先生による関数を使うのが手っ取り早い。
http://aoki2.si.gunma-u.ac.jp/R/AHP.html
他の一対比較法についても青木先生による関数がある。
サーストン法 http://aoki2.si.gunma-u.ac.jp/R/ThurstonePairedComparison.html
シェッフェ法 http://aoki2.si.gunma-u.ac.jp/R/ScheffePairedComparison.html
ブラッドレイ・テリー法 http://aoki2.si.gunma-u.ac.jp/R/ScheffePairedComparison.html
★外れ値
外れ値の検討 ([twitter:@isseing333] さん)
http://d.hatena.ne.jp/isseing333/20110901/1314883790
scratch-Rには mvoutlier パッケージの使用例。
http://eau.uijin.com/stats/anovaAssumptions.html
outliers パッケージのgrubbs.test関数の使用例。
http://sc1.cc.kochi-u.ac.jp/~murakami/cgi-bin/FSW/fswiki.cgi?page=R%A5%E1%A5%E2#p12
[twitter:@kingqwert] さんによる DMwR パッケージのlofactor関数の使用例。
http://d.hatena.ne.jp/kingqwert/20130103/p1
RjpWikiにもどなたかが外れ値の探索に便利な関数を載せてくれている。
http://www.okada.jp.org/RWiki/?R%A4%C7%A5%C7%A1%BC%A5%BF%A5%AF%A5%EA%A1%BC%A5%CB%A5%F3%A5%B0#h9b85c91
★Rによる独立成分分析
独立成分分析は次元縮約の方法の一つ。
Rでは fastICA パッケージで実行できる。
独立成分分析そのものについては下記が参考になった。
http://www.slideshare.net/zansa/121218-zansa13-for-web
主成分分析と独立成分分析の違いについては下記のページが参考になる。
http://meg.aalip.jp/ICA/index.html
Rでの実行については [twitter:@kzfm] さんの以下のエントリも参考になる。
http://blog.kzfmix.com/entry/1226406576
■
メモ
http://flowingdata.com/2010/08/31/how-to-visualize-data-with-cartoonish-faces/
http://blogs.iq.harvard.edu/sss/archives/2006/11/chernoff_faces_1.shtml
青木先生
source("http://aoki2.si.gunma-u.ac.jp/R/src/face_plot.R", encoding="euc-jp")
source("http://aoki2.si.gunma-u.ac.jp/R/src/face_data.R", encoding="euc-jp")
test <- cbind(1:9,50)
pos <- c(1,rep(2,17),1)
x <- face.data(test, pos)
old <- par(mfrow=c(3, 3))
for(i in 1:9) {
face.plot(x[i,])
}
par(old)
■
メモ
library(RFinanceYJ) library(aplpack) library(RColorBrewer) SINCE <- '2012-01-01' DATE.END <- '2012-12-31' Yahoo <- quoteStockTsData('4689.T',since=SINCE, date.end=DATE.END) DeNA <- quoteStockTsData('2432.T',since=SINCE, date.end=DATE.END) GREE <- quoteStockTsData('3632.T',since=SINCE, date.end=DATE.END) Mixi <- quoteStockTsData('2121.T',since=SINCE, date.end=DATE.END) CA <- quoteStockTsData('4751.T',since=SINCE, date.end=DATE.END) GungHo <- quoteStockTsData('3765.Q',since=SINCE, date.end=DATE.END) Intage <- quoteStockTsData('4326.T',since=SINCE, date.end=DATE.END) X <- format(DeNA[,1], "%m/%d") Y <- cbind( Yahoo[,7]/Yahoo[1,7], DeNA[,7]/DeNA[1,7], GREE[,7]/GREE[1,7], CA[,7]/CA[1,7], Mixi[,7]/Mixi[1,7], # GungHo[,7]/GungHo[1,7], Intage[,7]/Intage[1,7]) cols <- brewer.pal(7, "Dark2") par(xaxt="n") matplot( y=Y, lwd=2, col=cols, lty=1:6, type="l", xlab="date", ylab="per_close") par(xaxt="s") axis(1,at=1:length(X),labels=X) #abline(h=1:5,lty=2,col="gray") #legend(100, 4, c( "Yahoo","DeNA","GREE","CA","Mixi","GunHo","Intage"), # lwd=2,lty=1:6, col = cols) abline(h=seq(0.4,1.2,0.2),lty=2,col="gray") abline(h=1,col="gray") legend(1, 0.7, c( "Yahoo","DeNA","GREE","CA","Mixi","Intage"), lwd=2,lty=1:6, col = cols) library(animation) # http://www.imagemagick.org #Define path to convert.exe (http://www.imagemagick.org) ani.options( convert = shQuote('C:/Program Files/ImageMagick-6.8.1-Q16/convert.exe') ) test <- function(){ for (i in 1:nrow(DeNA)) { DAT <- rbind(DeNA[i,-1],GREE[i,-1]) faces(DAT,print.info=FALSE, face.type=1, main=Yahoo$date[i],labels=c("DeNA","GREE")) } } saveMovie(test(),interval=0.3,moviename="TEST",movietype="gif",outdir=getwd())
Chernoff face
WIKIPEDIA http://en.wikipedia.org/wiki/Chernoff_face
元論文(1973) http://www.apprendre-en-ligne.net/mathematica/3.3/chernoff.pdf
http://flowingdata.com/2010/08/31/how-to-visualize-data-with-cartoonish-faces/
library(aplpack)
http://www.wiwi.uni-bielefeld.de/com/wolf/software/aplpack.html
representation of the rows of a data matrix by faces. The features paramters of this implementation are:
1-height of face
2-width of face
3-shape of face
4-height of mouth
5-width of mouth
6-curve of smile
7-height of eyes
8-width of eyes
9-height of hair
10-width of hair
11-styling of hair
12-height of nose
13-width of nose
14-width of ears
15-height of ears.
http://aoki2.si.gunma-u.ac.jp/R/face.html
library(symbols)
symbol(iris,type="face",scheme=3,sortby=2,colin=5)
http://www.e-stat.go.jp/SG1/estat/List.do?bid=000001027245&cycode=0
アニメーションパッケージについては
http://d.hatena.ne.jp/EulerDijkstra/20120129/1327834316
http://www.slideshare.net/sleipnir002/animation10-11306609
http://d.hatena.ne.jp/ryamada/20110809/1312750005
http://www.slideshare.net/mickey24/r-de-animation
http://menugget.blogspot.jp/2013/01/producing-animated-gifs-and-videos.html#more
K-01 に KMZ INDUSTAR 50-2 50mm/F3.5 (M42) を付けて撮ったみた
先日購入したK-01が楽しくてたまらない。
何が楽しいって、マニュアルレンズで撮りやすいのが最高に楽しい。
で、ついにM42マウントのレンズにも手を出してしまった。
しかもロシア（ソ連）レンズ。
K-01にぴったりなパンケーキレンズ INDUSTAR 50-2。
もっとトイレンズっぽい写りかと思っていたら、いやいやスゲーよく写る！
今日はちょっと時間が取れたのでこの組み合わせで試し撮りしてきました。
風鈴
金魚
吽形
阿形
風車
注連縄
御神籤
手水
手水
手水
手水
窓
境内
空蝉
蛇口
ぐるぐるボケ
木漏れ日
っと、いい感じの写真だけ載せましたが、何しろ1959年に生産が始まった古いレンズです。
さらに元を辿れば1902年に設計されたTessarのデッドコピーなんだそうで100年以上前から使われているレンズ。
デジタル対応の最新レンズには敵わない部分ももちろんあるわけで。
逆光にはめっぽう弱く、こんな感じに盛大なフレアが出ます。
フードがあればだいぶましになりそうですが、そうすると小さくて薄いパンケーキレンズの利点が損なわれるような気もしますしねぇ。
