論理の流刑地

流罪に遭い戸惑う世捨人の雑記

【小ネタ】SPSSファイルにおける数値⇔値の対応関係をR上でも参照する

世界で10人くらいにしか需要がないだろう備忘録シリーズ。

Motivation

一般的に言って、SPSSのデータ(.sav拡張子のやつ)をRで利用するとき、read.spss()関数を使ってデータフレームに変換する。
しかしこのやりかたでは、SPSS上では参照できてた数値⇔値ラベルの対応関係が参照できなくて困ることがごくまれにある。

そんなときは、sjlabelledパッケージを使えばよい、ということが分かったのでメモ。
いやぁニッチな需要に応えるpackageを作ってくれる人が世界にはおるんやなぁ...

How to

基本的な使い方は以下のURLがわかりやすいので要参照。

library(sjlabelled)
data1 <- sjlabelled::read_spss( file= "C:/Users/XXX/example.sav" , enc="UTF8")
 
#データ全体の変数ラベルを取得
var_lbl <- sjlabelled::get_label(data1)

#特定の変数の値一覧を取得
q1_values <- sjlabelled::get_values( data1$q1)

#特定の変数の値ラベル一覧の取得
q1_labels <- sjlabelled::get_labels( data$q1)

#[応用] 数値と値ラベルの対応がわかるかたちで値ラベル一覧の取得
q1_labels_v2 <- sjlabelled::get_labels( data$q1, value="p")

普段は、foreignパッケージのread.spssをつかうことが多い思う。
だがsjlabelledではread_spssなので注意すべし。

Enjoy!

【小ネタ】geom_point()やgeom_line()で二つ以上の変数の組み合わせを用いてgroup化する方法

This is a note for 備忘

ggplot2のgeom_line()で折れ線グラフをかくときに、系列を二つの変数の組み合わせ(性別×年齢とか、学歴×人種とか)でとらえたい時がある。
このやり方を知らなくて、ずっと、二つの変数を組み合わせた新たな水準変数を作っていたが、それは不要だと今更気づいたという話だ。

stackoverflow.com

このQ&Aにあるように、引数のmappingの指定のなかで、group=interaction(var1, var2)を使えばよい。
仮想のデータdf1で、年代yearをX軸として、性別sexと人種raceの組み合わせを系列としてY軸に所得incomeをplotしたいときは、
以下のように書くのだ。

library(ggplot2)
pl1 <- ggplot(data = df1 , mapping=aes( x = year , y = income , group=interaction( sex , race) , colour = race , linetype= sex , shape = sex) + geom_line() + geom_point()
plot(pl1)

このinteraction(sex, race)の部分により、系列を2水準の組み合わせで定義している。
線の色だったり点の形だったりは、この系列単位でなく元の変数(sex, race)により定義できるのがウレシイところである。

【小ネタ】Rでclipboard利用 in Windows

備忘用。わりかし実用的だが忘れがちな小ネタ。
Macだと若干やり方が違うので注意(という自分用備忘)

Motivation

Rでデータフレームをつくって、それをExcelにのせたいときがある。
また、Excelやブラウザ*1上でコピーしたtableをそのままdata frameとして使いたい時がある。

そんなとき、お手軽なのがクリップボードを使うことである。

参考記事

How to

Rのデータフレームをクリップボード

write.table(test_data,file="clipboard-16384",sep="\t",fileEncoding = "CP932")

やりかたは簡単。write.tableの引数fileに"clipboard"を指定するだけ。
ただし、自分は以前の記事で書いたようにデフォルトの文字コードUTF-8にしているので例えばExcelに移す時に文字化けが発生しがち。
したがって、fileEncoding="CP932"を利用することで文字化けを回避している(参考URL)。

また、データフレームの容量によっては"clipboard buffer is full and output lost"というおぞましいErrorが出る(とくに文字を含んだときに顕著)。
これを回避するためには、引数fileに"clipboard"ではなく"clipboard-16384"と指定する(参考URL)。
16384は上限ってだけでそれ以下ならばよい。

クリップボードからRのデータフレームへ

まぁ対象はなんでもいいんだけど、ブラウザからのほうが便利さ&目新しさはあると思うので、
とりあえず藤井七段の対局結果のページから、「最近の対局結果」のテーブルをコピーした前提で話をすすめる。

コードは以下の通りだ。このWebページはUTF-8に設定されているようだがなぜかfileEncoding="CP932"じゃないとエラーを吐く。理由はわからない。

rdf <- read.table(file="clipboard-15000", header=T, sep="\t",fileEncoding = "CP932")

すると、綺麗に元のテーブルがdata.frame化されている。

f:id:ronri_rukeichi:20190625122953p:plain

Enjoy!!

*1:時間があるときは再現性の問題を含めてもExcelのほうはread.csvなどでやるのが王道だと思うが、このブラウザコピーのほうは割と有用だと思う。スクレイピングのコードは結構書くのに時間がかかるし。

ブール代数分析をふわっと

ふわとろ高級オムライスがたべたい

Motivation

「どうせ簡単だろ」とタカを括ってあんまり勉強してなかった手法シリーズなので。
なんとなくの感触で過信してちゃんと時間かけないのは恥ずべきことなんだよな...

てことで下らへん見て勉強する...

◆文献list

  1. 石田淳(2017)『集合論による社会的カテゴリーの展開』(第3-4章)
  2. Caren & Panofsky , 2005 , "A Technique for Adding Temporality to Qualitative Comparative Analysis"(URL)

基本定理

全部網羅的に書いてもあれなので、自分が躓いたところだけメモ

まず基本的な定理(前提となるもの)

交換律: x+y = y + x, x\cdot y = y \cdot x
分配率: x + ( y \cdot z) = (x + y ) \cdot ( x + z)
同一律x\cdot I = x , \ \ x + O = x
補元律:x+\bar{x} = 1, \ \ x\cdot \bar{x} = O

基本的定理のコロラリー

べき等律: x+x = x, \  \  , x\cdot x = x
有界律: x + I = I , \ x \cdot O =O
吸収律: x\cdot(x+y) = x, \  x + (x \cdot y ) = x
結合律(省略)、対合律(省略)

覚えとくと便利なこと

相対原理
ブール代数の任意の性質は、+と・,x\bar{x}を入れ替えても成立する

ブール代数における順序関係
 xy=x \Leftrightarrow x+y =y
が成り立っているとき、順序関係 x \leq yが成り立っているとする。

ブール関数の簡単化(いわゆるcsQCAのやってることはだいたいコレ)

「ブール関数の簡単化」とは、標準積和形(最小項*1の和によるブール関数の表現式)からスタートして、できるだけ項数と各々の項を構成する変数の数が少ない式を求めること、である。

その過程は、

  1. 主項の導出
  2. 最小積和形の導出

の2ステップを踏む*2

主項の導出

Definition:
まず、最初の段階では標準積和形から、最小化定理とベキ等律を繰り返し適用して、これ以上簡単化できないところまで式を変形する。
こうして得られる式の各項を主項(prime implicant)という(石田 2017: 65-66)

具体的には

  1. ハミング重み(最小項に含まれる1の数)の算出
  2. ハミング重みが1つだけ異なる最小項の総当たり比較→最小化定理による変数の削減
  3. 最小化定理が適用できなくなるまで、第二段階の変数削減を繰り返し行う

という形で行う。
 D = ABC + ABc + Abc + abc = AB + Ac +bcみたいな感じで。

最小積和形の導出

Definition
最小積和形とは、主項数が最も少なく、かつ主項数が同じであれば変数の数が最も少なくなるような表現式のことである
(石田 2017: 66-67)

この導出は、主項選択表を用いることでおこなう。

主項選択表とは、表側に前段階で導出した主項、表頭に最小項をおいた表である。

[STEP1]
行(主項)が列(最小項)を包含している場合、該当するセルに〇をつける。
i行j列のセルの場合「最小項j\leq主項i、が成り立ってるとき〇をつける」ということである。

[STEP2]
必須項(最小項をもれなく包含するために不可欠な主項)に◎をつける
「主項選択表を縦にみて、1つしか〇のつけられない最小項が必須項である」(同:67)
必須項を選ぶことで、求める最小積和項となる(必須項がない場合や複数の組み合わせがあるときもある)。

発展的な話題

上で挙げたCaren&Panofsky(2005)はQCAに時間軸(イベントの順番)を導入したTQCAなるものを提唱している。
一通り論文を読んだが、①順番を定義すること、②「イベントが起きなかった」ことの取り扱いがやや特殊であること、を除けば難しくはなさそうだ。

おわりに

この記事を書き終わったら、羽生さんが若手の伸び盛り佐々木大地五段相手に深夜の千日手指し直し局を制しておられた。
わたしのような老兵もまだ頑張らねばなー。


Relaxing Breath of the Wild music with rain

Enjoy!!

*1:ブール関数を構成するブール変数or補元で構成される項。3変数だったらXyZとかxyzとか。

*2:クワイン・マクラスキー法と呼ぶらしい

R in Windowsの環境設定(.Rprofileまわり)

力を入れることの何倍も力を抜くことのほうが難しい。だが重要。

Introduction

備忘 of 備忘。

久しぶりにRに触る仕事をしているが、
なんか無理やりMacWindowsにファイルを移して環境を構築した*1ので、色々齟齬がでている(主に文字コードまわり)

やったこと

ファイル読み込みの際の文字コードの設定

基本的にMacで作ったコードファイルはUTF-8で保存されているため、options()$encodingを変更する。
具体的には、以下のコードを.Rprofileに追記する

options(encoding="UTF-8")

日本語でのプロット用の設定

ちなみに上の修正をしようと、.Rprofileを触ったら、「setHookとquartzFontががなんちゃら...」というエラーが出た。
前回Macで環境構築したのが昔すぎて忘却の彼方なので、なんだっけこれ?と調べた。
下の記事等を見るに図の出力の際に、日本語を使うための設定だったようだ。

qiita.com
marui.hatenablog.com


ちなみに、これはMacでの設定なので、Windowsでの日本語プロットについては、以下参照して設定するのがよい(参考URL1, 参考URL2
自分は基本的にメイリオがすきなので、下のコードではじめからメイリオを使えるようにしてる。

setHook(packageEvent("grDevices","onLoad"), function(...){
  grDevices::windowsFonts(Meiryo=grDevices::windowsFont("Meiryo"))  
} )

意味としてはパッケージgrDevicesがロードされたときに呼び出される処理として、windowsFontsのMeiryoパラメータとしてメイリオを設定している。
なんかイメージとしてはjavascriptのイベントトリガとかが近いのかなと思った(適当)。

これにより、par(family=Meiryo)とか、いきなり書ける。

おまけ

.Rprofileに書く作図上の文字コード設定上の処理と、その意味については以下のgithubで公開されている.Rprofileのコメントが参考になる

.Rprofile · GitHub


Enjoy!!


*1:今更初めてRstudio使ってるけどあんまり自分の感覚にfitしないな...

どのようにして分位点回帰の推定は線形計画法でなしうるのか

食費を節約したいし痩せたいけどついついコンビニでドリアとか買っちゃう社会の闇

Introduction

分位点回帰モデルに関する解説を読むと、「推定のところでは線形計画法を使う」と書いてある(Koenkerの有名な解説書等)。
分位点回帰の推定において、何が最小化されているかというのはそれなりにとっつきやすいんだが、
それがなぜ線形計画法で解けるのかということは何となくわかった気になってそのままにしていた。

なぜか急に気になって電車の中でスマホで調べたら、
以下のStackExchangeでの質問への回答がなかなかにわかりやすかった*1ので、自分用にメモしておく。

stats.stackexchange.com

線形制約問題自体については、最近たまたま見つけた以下の岩波講座「最適化法」の第2章(著者:今野浩=ヒラノ教授)がこれまで読んだ中で一番過不足なく説明されていると感じた。さすが第一人者。

前提

線形計画法の標準形

線形計画問題の標準形は、以下のようにあらわされる。

 min \ \  c^T z  , Az = b

ここで、 z ,\ c , \ bはいずれもk × 1、Aはn × kの行列であるとする。
また、 z >=0 となっている。

解き方(単体法、内点法など)はともかく、この形で示すことができれば、
線形計画法として解にたどりつく(無限解や解が存在しないという判定も含め)ことができる。

分位点回帰における「最小にすべきもの」

一般に、回帰モデルとは、個人iに関する従属変数y_iを独立変数と回帰係数の線形結合と誤差の和として表すものである。
y_i = \beta^T x_i + \epsilon_i

さてここで、教科書的な説明を振り返っておくと、通常の線形回帰モデル(OLS)においては、
「誤差の二乗和を最小化する」という条件を満たすように回帰係数が求められる。すなわち、
   \beta \in R^k ,  ( s.t.  \ \ min \  \Sigma_i^N \epsilon_i^2 )

となる。一方、分位点回帰モデル(以下QRMと略記)においては、「誤差の重み付き和を最小化する」という条件を満たすように回帰係数が求められる。すなわち、

 \beta \in R^k ,  ( s.t.  \ \ min \  \Sigma_i^N  \rho_i |\epsilon_i |)

ここで、weightである\rhoは以下のように定義される
 \rho_i = \tau \ (  \rm{if} \ \epsilon_i >=0 ) , 1 - \tau \ ( \rm{if} \  \epsilon_i < 0 )

ここで \tauは、独立変数が与える影響を推定したい分位点であり、0 < \tau < 1である。
よりわかりやすく日本語で説明すると、各ケースiについて「観測値 - 予測値」の符号によって重みを変えているということである。
たとえば20%点に対してのQRの場合、観測値が予測値を超える場合の重みは0.2、そうでないときは0.8となる。

何がわからなかったのか

これまで何がわからなかったのか、すなわち「分位点回帰は線形計画法で解ける」と言われてもどこが引っ掛かっていたのかというと、

  1. そもそも回帰係数\betaが決まらないと\epsilon_iの符号は決まらないのに、「\epsilon_iを最小化するように係数を推定する」というのは困難なのでは?(相互依存性の問題1)
  2. \epsilon_iの符号が決まらないとウェイト\rho_iも決められないが、\betaが未決定だと前述のように\epsilon_iもわからないのでウェイトがどうなるのかは事前にわからないのでは?(相互依存の問題2)

ということで、つまるところ
回帰係数\betaの推定値 ⇔ \epsilon_iの値 ⇔ \epsilon_iの符号 ⇔ ウェイト\rho_iの値
という相互依存関係をどのように扱うのかということがわからなかった
ということである。

相互依存性への解決法としての同時的決定

今回勉強してみての結論からいうと、上記の問題は
\epsilon_i\beta\rho_iを同時決定されるパラメータとして解く」
によって解消される。

線形制約側の話

まず線形制約側について考える。回帰モデル
 Y_i = \beta^T X_i + \epsilon_i
をさきほどみた線形制約Az=b, z >=0の形に変えることを考える。
とりあえず縦にn個つなげて、
 Y = \beta^T X + \epsilon
とする、左辺右辺ともに要素n個のベクトルとなる。

さてここで行列Aにあたるのは、観測された独立変数Xとなる。
「変数」という言葉で惑わされそうになるが、観測値は固定なので、スカラの制約行列となるのである。
そして線形制約式の右辺bにあたるのは、従属変数列Yである。

したがって、線形制約式における変数zにあたるのは、残った二つ。
係数\betaと、誤差\epsilonである。

しかし、これらは定義上負の値をとりうる(これを線形計画法の用語で「自由変数」という)ので、
以下のように二つの非負変数ベクトルの差として表現しなおすこととする。
 \epsilon = u - v , \ (u >=0 ,\  v >= 0 )
 \beta =  \beta' - \beta'' ( \beta' >=0 , \ \beta'' >=0 )

すると、
 Y = \beta^T X + \epsilon
 \iff \ Y = ( \beta' - \beta'' )^T X  + u- v
であり、ここで次元nの単位行列Iを導入すると
\iff  \ Y = ( \beta' - \beta'' )^T X  + (u- v)^T I
となるから結局、
 A = \begin{pmatrix} X & -X & I &  -I \end{pmatrix}
 z = \begin{pmatrix} \beta' & \beta''  & u & v \end{pmatrix}
とおくと、
 Az = Y という綺麗な線形制約の形になっているのである。

最小化対象式側の話

残りは \textrm{min} \  c^Tz のほうだが、
変数zは先ほど定義したので、c(スカラのベクトル)を適切に定義すればいい
 c = \begin{pmatrix}  \textbf{0}  & \textbf{0} & \tau 1_N & (1- \tau) 1_N \end{pmatrix}
ただしここで1_NはN×1ですべての要素が1の行列である。

cの定義から明らかであるが、最小化対象(=cとzの内積)には、回帰係数が出てこないようになっている。
実測値>予測値となるときは \epsilon_i = u_i 、実測値<予測値となるときは\epsilon_i = - v_iとなる*2ので、
 c^T z =  \begin{pmatrix}   \tau 1_N & (1- \tau) 1_N \end{pmatrix}^{T} \begin{pmatrix} u & v \end{pmatrix}
は最小すべき重み付き絶対値和 \Sigma_i^N  \rho_i |\epsilon_i |になっている。

あとは、これを普通の線形計画問題として解けばQRの推定値が得られるのである。
めでたしめでたし。

最後に

なんか「一見相互依存に見えていてデッドロックを起こしてしまい)解けない問題を同時決定問題として解く」というのは、
統計学や工学など理系の問題だけでなく、社会科学や、政策決定などのより一般的な問題解決においても適宜必要な知恵ではないかと思った。
囚人のジレンマがジレンマなのは同時決定できる超越的な主体が存在しないからだよね、って話である。

線形計画もわかったような気になってまた忘れるを繰り返している気がするが
時間があったら前述の岩波講座の忘れちゃいけない部分だけまとめときたい。

Enjoy!

*1:なのに質問者は「長くてわからん」って言っててああ無情。

*2:厳密にはこれは別途証明する必要があるが省略する

最近印象に残った言葉

社会は人に一つの人格を望む。
人格が一貫していればしているほど誠実な人間だと捉えられる。
だから皆「表向きの自分」を演じてそこからはみ出る部分を隠したがる。
隠してきたものを消してやれば一貫した人生になる。
反対に暴き始めれば色んなものが壊れる。

ーードラマ「dele」第7話

プラトンの考えでは真理は永久のもので、静的抽象概念の中にのみ存在する。真理はけっして実際的現実的な現象の中には存在しない。これらは流転の過程にあるのだからね」
「でも、私たちの真理の概念は全く逆よね」とケートが言った。
「私たちにとって真理は現実の多数のかけらから成り立っている。流転と自己変化はその本質よ。変化は重大な真実だから、私たちの考えでは過程こそ事物の本質なのよ」
「科学者らしい言い方だね」とアームブラスターが言った。

ーーJ・ジェイコブス「市場の倫理 統治の倫理」pp.367-368

市場の倫理 統治の倫理 (ちくま学芸文庫)

市場の倫理 統治の倫理 (ちくま学芸文庫)

「僕は時々、ラ・マシアで話をしてくれと言われるけど、よくこういう例を引き合いに出すんだ。
『毎晩寝るときに、サッカーが好きかどうかを自分自身に問いかけてほしい。今この瞬間にボールを抱えてベッドから起き、プレーができるだろうか』とね。
もし答えが『ノー』なら、サッカー以外に目を向けたほうが良いと思う」

ーーグイレム・パラゲ「知られざるペップ・グアルディオラ」p.72

知られざるペップ・グアルディオラ サッカーを進化させた若き名将の肖像

知られざるペップ・グアルディオラ サッカーを進化させた若き名将の肖像