Menzerath-Altmann法則と一般化ガンマ分布(Eroglu, 2013)
とても興味深いこと。Eroglu(2013)は,Menzerath-Altmann法則と一般化ガンマ分布の関係について説明している。
Eroglu, S. (2013). Menzerath-Altmann law for distinct word distribution analysis in a large text. Physica A, 392, 2775 – 2780
Menzerath-Altmann法則は,
(a)
という風に書ける。
ここで,b = 0,c ≒ 0のとき,
(b)
となる。これは普通の指数関数。
で,b ≒ 0,c = 0のとき,
(c)
これはZipfの法則。
ここで,一般化ガンマ分布の確率密度関数を(d)のように書くとする(parameterizationもいろいろあるが…)。
(d)
β = 1のとき,(d)式はガンマ分布の確率密度関数になるし,β=θのとき,ワイブル分布の確率密度関数になる。その上,ガンマ分布は,母数の値によっては,指数分布にもなるし,アーラン分布にもなる。
で,いちばん大事なんだけど,のとき,(d)式は(a)式になるらしい(Eroglu, 2013, p. 2776)。なるほど。
…なんだか,とっても安心する。
Menzerath-Altmann法則をRのnls関数で
Menzerath-Altmann法則(またはMenzerath法則; e.g., Altmann, 1980; 基本的アイデアはMenzerathが発見,Altmannが定式化)とは,ある言語学的単位の増加が,その単位の構成要素の減少に帰着すること。またはその逆。
もっと一般的にいうと,Altmann(1980)は,「ある言語学的単位の長さは,その構成要素の長さの関数である」としている。たとえば,ある形態素長(形態素内の音節数などとして測定する)の値が高ければ,その形態素内の平均的な音節長(音節内の音素の数などとして測定される)の値は低くなるということ。ある言語において,1音節からなる形態素の平均的な音素数は3,2音節からなる形態素(ひとつつあたりの)の平均的な音素数は2.5,3音節からなる形態素の平均的な音素数は2…というように。
それで,基本的にはさまざまな単位,たとえば時間などでもこの法則が当てはまるとされ,実証研究も少なくないし,ゲノムとかタンパク質とか,そういう分野にも応用されているんだそうだ。で,その関係は,簡単な方程式で表されると。
もっとも一般的な式は,
(a)
ここでのyは構成要素(小さい方)の長さ,xは対象とする言語学的単位の長さ,a,b,cは自由母数。eはネイピア数。この式を使い,観測から,a,b,cを推定する,というわけだ。
一番単純な方法だと,最小二乗法を使用することができる。
Rには,非線形最小二乗法用の関数nlsがあるので,これを使えばいい。Altmann(1980)にインドネシア語の形態素と音節の関係についてのデータが示されている。これは,こんな感じ。
Gabriel Altmann (1980). "Prolegomena to Menzerath's law". Glottometrika. 2: 1–10.
Morpheme Length(Syllables) | Syllable Length(Phonemes) |
1 | 3.10 |
2 | 2.53 |
3 | 2.29 |
4 | 2.12 |
5 | 2.09 |
やってみる。
#データ入力 s<-1:5 m<-c(3.1,2.53,2.29,2.12,2.09) #可視化 plot(s,m,xlab="Morpheme Length",ylab="Syllable Length",pch=20,cex=2) #モデルをフィット fit<-nls(m~a*s^(-b)*exp(1)^(-c*s),start=list(a=3,b=0,c=0)) #サマリーなど summary(fit) AIC(fit) BIC(fit) #決定係数 cor(m,predict(fit))^2 #予測値をプロット plot(s,m,xlab="Morpheme Length",ylab="Syllable Length",pch=20,cex=2) points(s,predict(fit),col=2,pch=20)
赤がモデル理論値。もちろん原著の推定結果ともほぼ一致している。
このデータでは,以下のようになった。
(b)
この法則がどれだけ法則と呼ぶに相応しい一般化可能性をもつかといった面と,具体的に(a)式がどれだけ観測にフィットするかということは置いといて,ある単位と,その構成素の関係に一般的な法則をもとめる姿勢というのは,とても正当な考え方だと思う。
ジニ係数とワイブル分布・対数正規分布・ガンマ分布の母数の関係
どうでもいいことなんだけど,ジニ係数はいくつかの分布であれば,分布の母数から直接的にもとまるんだそうだ。
ワイブル分布,対数正規分布,ガンマ分布を例に取ると,ワイブル分布とガンマ分布なら形状母数k(またはα)のみから,対数正規分布ならσのみから直接ジニ係数がもとまる。これらを計算するRの関数を作ってみる。
#ガンマ分布 gamma.gini<-function(a){ gini<-gamma((2*a+1)/2)/k/gamma(a)/sqrt(pi) gini } #ワイブル分布 weibull.gini<-function(k){ gini<-1-2^(-1/k) gini } #対数正規分布 lognorm.gini<-function(s){ gini<-2*pnorm(s/2*sqrt(2))-1 gini }
あるデータの分布がワイブル分布に従っているとして,その形状母数kが2だったとしたら,
weibull.gini(2) 0.2928932
この関数にはベクトルを入れられるので,こんなグラフも作れる。
x<-seq(0.1,5,.01) plot(x,weibull.gini(x),xlab="Parameter",ylab="Gini Coefficient",type="l")
ワイブル分布では,形状母数が大きくなればなるほどジニ係数は小さくなる。
…どうでもいいことなんだけど,ワイブル分布の形状母数について,MCMCのサンプルを得て,ジニ係数の信用区間を得ることもできるんじゃないだろうか。たとえば,
set.seed(1) d<-rweibull(1000,3,1) library(MCMCpack) llf<-function(beta,x){ sum(log(dweibull(x,beta[1],beta[2]))) } m<-MCMCmetrop1R(llf,theta.init=c(3,1),x=d) gini<-weibull.gini(m[,1]) plot(density(gini),type="n",xlab="Gini",ylab="Density",main="") polygon(density(gini),col=rgb(0,0,1,.2)) #信用区間 quantile(gini,c(.025,.975)) 2.5% 97.5% 0.1963506 0.2140616
こんな感じで。
みたいな。
GAMLSSとBAMLSSで分布をデータにフィット
GAMLSSとBAMLSS
新幹線に乗ってコンビニにいくようなものだけど,GAMLSSとそのベイズ的モデルであるBAMLSSを使って,観測のデータに対して分布をフィットさせる。
GAMLSSというのは,「位置,尺度,形状,それぞれの母数のための一般化加法モデル(Generalized additive model for location, scale and shape)」のこと。詳しくはwikiや,開発者らのページを。
Generalized additive model for location, scale and shape - Wikipedia
gamlss | for statistical modelling
大雑把にいえば,分布ベースの回帰モデルのひとつ。
基本的には,(ひとつの)応答変数になんらかの分布を仮定する。正規分布とかガンマとか。
たとえばガンマ分布だったら,α(形状)とβ(比率)母数のふたつがあると。
で,応答変数のそれぞれの値は,これらの母数によってランダムに発生していると。
GAMLSSでは,これら,それぞれの母数の値自体を予測する。
予測変数の値という条件における,応答変数が従うであろう分布の母数を予測して,そこから応答変数ができてるみたいな。
この母数に対する予測は,線形な関係であってもいいし,加法モデルというくらいだから,スプラインとかでもいい。
基本,4母数までの分布ならできるらしく(μ,σ,ν,τ),Rのgamlssパッケージには莫大な数の分布がサポートされている。
変量効果も入れられる。とっても自由。
BAMLSSというのは,これとほぼ同じようなことをベイズ的にやるんだそうだ(細かい事前分布の設定とかアルゴリズムはよくわからない)。
これが面白いのは,平均値を予測するだけでなくて,散布度や形状を予測することができるようになるところ。
たとえば,実験条件Aは,応答変数の平均値には影響を及ぼさないけど,分散を大きくする効果をもつ,とか。
歪んだ分布を仮定して,実験条件Bは,応答変数の形状に影響を及ぼす,とか。そういうモデリングができるようになる。
「分位点回帰のパラメトリック版を非線形でもできるようにしてる」,みたいな感じ。
…それは置いておいて,このパッケージは特に便利なので,まずは単変量の分布の母数をこれらのパッケージで推定する。
ある分布に従う単変量の母数の推定
MASSパッケージのfitdistrやfitdistrplusパッケージでもoptim関数でも,なんなら場合によってはglm関数でも,同じことをやることができる場合は多いのだけど,gamlssパッケージやbamlssパッケージは,普通に単変量のデータの分布を推定できる。
まずは,データにワイブル分布を当てはめて母数を推定する例をやってみる。
#読み込み library(gamlss) #データの作成 set.seed(92) d1<-rweibull(1000,2,1) #gamlssの基本関数を使って推定 #最初の回帰式がμ(最初の母数),次の回帰式をσと当てている使用 WEIは,ワイブル分布 m1<-gamlss(d1~1,sigma.formula=~1,family="WEI") #これでダン #サマリーしてみる summary(m1)
切片の推定値はそれぞれ,μ(ここでは尺度母数) = -0.01,σ(ここでは形状母数) = .677くらい。
これにはリンク関数としてlogが噛んでいるので,係数をexpする。
すると,まあ,1と2くらい。よく推定できている。もちろん,これはただの最尤推定と同じなので,他のパッケージとかとも一致する。
で,もっと楽な関数があって,gamlssMLという最尤推定の関数。
m2<-gamlssML(d1,"WEI") summary(m2)
これもおなじことになる。わざわざ切片からの回帰式とか書かなくてもよい。
もっと便利なのがあって,histDistという関数。
これは,ヒストグラムとフィットした確率密度曲線を描いてくれるだけでなくて,上のサマリーも返してくれる。
histDist(d1,"WEI")
図はこんな感じ。オプションもあるので,色々いじれるけど,自分はやらないだろうな。
BAMLSSでも基本は同じ。少なくともこの状態で,それぞれの事前分布は基本的に無情報っぽい。
MCMCサンプラーはJAGSとか色々指定できるみたいだけど,デフォルトではGMCMCというのを使っているそう(これはよくわからない)。
library(bamlss) #いろいろなオプションを省略して… m3<-bamlss(list(mu=d1~1,sigma=~1),family="weibull") summary(m3)
それぞれの母数の95%信用区間(パーセンタイル法)も出してくれる。
以下のようにしたら,MCMCサンプルも可視化できる。
par(mfrow=c(1,2)) plot(density(exp(unlist(m3$s[,1]))),main="Scale",xlab="Value") plot(density(exp(unlist(m3$s[,5]))),main="Shape",xlab="Value")
簡単なモデル
えっと,せっかくだから,ただ生データをフィットさせるだけでなくて,GAMLSSやBAMLSSの中でも,もっとも簡単な回帰モデル的なものをやってみよう。
アイリスデータの,がく片の長さ(sepal length)を応答変数,あやめの種類(3水準)を予測変数としよう。がく片の長さは正規分布に従うとして,種類は,がく片の長さの平均と標準偏差にそれぞれどのような影響をあたえるだろうか。ここであやめの種類はfactor型。
データとしてはこんな感じ。
見た感じ,平均はsetosa < versicolor < virginicaという大きさのようになっている。標準偏差も同じ。
普通の線形モデルなら,
summary(lm(iris[,1]~iris[,5])) anova(lm(iris[,1]~iris[,5]))
みたいにする。
一方,gamlssでは,こんな感じにする。
m4<-gamlss(iris[,1]~iris[,5],sigma.formula=~iris[,5],family="NO") summary(m4) ****************************************************************** Family: c("NO", "Normal") Call: gamlss(formula = iris[, 1] ~ iris[, 5], sigma.formula = ~iris[, 5], family = "NO") Fitting method: RS() ------------------------------------------------------------------ Mu link function: identity Mu Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.00600 0.04935 101.44 <2e-16 *** iris[, 5]versicolor 0.93000 0.08751 10.63 <2e-16 *** iris[, 5]virginica 1.58200 0.10179 15.54 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ------------------------------------------------------------------ Sigma link function: log Sigma Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.0528 0.1000 -10.528 < 2e-16 *** iris[, 5]versicolor 0.3814 0.1414 2.697 0.00783 ** iris[, 5]virginica 0.5900 0.1414 4.172 5.2e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ------------------------------------------------------------------ No. of observations in the fit: 150 Degrees of Freedom for the fit: 6 Residual Deg. of Freedom: 144 at cycle: 2 Global Deviance: 206.9715 AIC: 218.9715 SBC: 237.0353 ******************************************************************
普通の線形モデルと違って,標準偏差の大きさの違いも同時にモデリングしている。
まあ,見た感じの結果になっている。これは,こんな感じに可視化するとわかりやすい。
term.plotという関数もある。
par(mfrow=c(1,2)) term.plot(m4,"mu",main="Mu") term.plot(m4,"sigma",main="Sigma")
もちろん,診断プロットもある。
plot(m4)
それぞれの予測された母数から,確率密度曲線を描くとこんな感じ。
もちろん,標準偏差には影響を及ぼさないとか,標準偏差には影響を及ぼすけど,平均値には及ぼさない,といったモデルを作り,適合度を比較したり尤度比検定をしたりすることもできる。
まあでも,これ予測変数が因子型だったら,教師ありの混合分布モデルみたいなもん。
blmssパッケージでは,
m5<-bamlss(list(mu=iris[,1]~iris[,5],sigma=~iris[,5]),family="gaussian")
で,(当然だけど)同じような結果を得ることができる。MCMCサンプルを可視化したり,信用区間を構築してもよい。
とっても便利。
平均は変わらないけど,標準偏差だけ違う2群のモデル
回帰モデルは自由に作れるので,ある因子が標準偏差だけに影響を及ぼしているモデルとかも作れる。
たとえば,こんなデータだとしよう。
group<-as.factor(c(rep("A",100),rep("B",100))) score<-c(rnorm(100,50,10),rnorm(100,50,30)) d2<-data.frame(group,class) m6<-gamlss(score~1,sigma.formula=~group,data=d2) summary(m6) term.plot(m6,"sigma",main="Sigma")
結果はこんな感じ。
ちょっと複雑なモデル
えっと,実質的に意味があるかはわからないけど(ないだろう),がく片の長さの分布を正規分布だと仮定して,がく片の幅と花びらの長さから,その予測値という条件下における分布の母数(平均,標準偏差)を推定するとしよう。で,その予測値と予測される母数の関係は非線形だと。
#ここでは罰則付きスプラインの一種 m7<-gamlss(Sepal.Length~pb(Sepal.Width)+pb(Petal.Length),sigma.formula=~pb(Sepal.Width)+pb(Petal.Length),data=iris,family="NO") summary(m6) par(mfrow=c(2,2)) term.plot(m7,"mu",main="Mu") term.plot(m7,"sigma",main="Sigma")
可視化するとこんな感じ。
がく片の幅と花びらの長さが長くなれば,がく片の長さの分布における平均はそれぞれ単調に増加しそう(花びらの長さはちょっと指数的に見える)。
がく片の幅は,がく片の長さの分布における標準偏差とはあまり関係なさそうだ。しかし,花びらの長さは,がく片の長さの分布における標準偏差に非線形の影響を及ぼしているかも(意味ないと思うけど)。…みたいな。
で,これに種類を変量効果と考えて…というようにやっていくと,大分自由なモデリングができそう。
(最尤推定による分布のフィットから書き始めて,やっとここまで来たって感じ)
オンライン外国語学習における実測学習時間の格差:ジニ係数とローレンツ曲線
同じ学習教材を与えても,それを使用して勉強するひともいれば,勉強しないひともいる。格差,ということばが正しいかはわからないけど,たとえば大学生が,ある一定期間内において,大学が提供するオンライン外国語学習教材を使用して学習する時間は,ちょうど所得の分布に似た分布を示す。
これは歪んだ分布になる。(ワイブル分布の例;草薙(2017)は個人のWBTに対する個人のログイン時間は,ワイブル分布によって近似を与えられることを報告している)
これについて,ローレンツ曲線を描くとこんな感じ。
ちなみにこのデータのジニ係数は,だいたい0.30くらい。おおざっぱにいって,日本人の所得のジニ係数と同じくらい(?)。
全体のアクセス時間に対して,上半分の人が7割,下半分の人が3割の時間くらいといった比になる。下二割の人は,全体の7%,上8割の人の全体の93%。
この係数を減らすというのも,ひとつの業務的観点になりうる。
JansonとOllssonのイオタ
JansonとOllssonのイオタは,評定変数が多変量で複数評定者間の信頼性ないし一致率を評価する係数。Rでは,irrパッケージのiota関数を使う。
たとえばルーブリックを使用した評定を二人が行い,観点が5つあるとか,そんなときに使用する。
#こんなデータに [[1]] Rater1 Rater2 1 Old Old 2 New New 3 Old New 4 New New 5 Old Old 6 New New 7 New New 8 Old Old 9 New New [[2]] Rater1 Rater2 1 Slow Slow 2 Fast Fast 3 Slow Slow 4 Fast Fast 5 Slow Slow 6 Fast Fast 7 Slow Fast 8 Slow Slow 9 Fast Fast iota(dat)[] $method [1] "iota for nominal data (2 variables)" $subjects [1] 9 $raters [1] 2 $irr.name [1] "iota" $value [1] 0.775 $detail NULL
キーボード用リストレストを作ってしまえ
あまりタイピングが上手じゃないのがいけないのだけど,よく腱鞘炎になってしまう。私みたいにもの書きをまったくしないやつに限ってすぐ腱鞘炎になる。そんな軟弱な自分に嫌悪感を覚えながらも,やれキーボードが悪い,そしてリストレストが悪い,この企業の製品はダメだのと因果関係がそれほど自明でないことをいってしまうのが人間ってもの。本当は自分が軟弱なだけなのにね。自分の軟弱さを認めるくらいなら,顔の見えない製造企業を責めるほうが心が休まるものね。アマゾンのレビューもそういうのばっかり。
キーボードについても,世の人はやれメンブレンでも構わぬとか,パンタグラフで十分だの,メカニカルでないと死ぬ,それも赤軸がいいの,青軸がいいの,俺は黒軸でないと触らないだの,入門は茶軸がいいよとか,いいや無接点が最強だの,配列についてもUS配列でないとコードが書けないだとか日本人なんだから日本語配列でいこうぜとか,HHKBの配列でないとダメだとか,さらに細かいとControlはaの左だとかFnは右左に置けだとか,BSを一段下げろだとか矢印キーは甘えだとか,そもそもフルサイズキーボードがいいに決まっているとか,テンキーレスがいいだとか,いやいや60%コンパクトサイズがいいとか,もう開き直って40%キーボードでいいとか。
もはやリストレストでさえも,やれ幅は30cmがいいの奥行きは8cmがいいの,高さはどうのこうの,材質はどうの…
でもこれって難しい,自分は,HHKB Lite2 US配列と全く同じ配列で,フレームレスの赤軸(Minila ではいけない)が1万円きったら6台買おうと思う。しかしそういうのは存在しない。癖の強い自分みたいなものが欲しいものが企業にとって商機だとは限らない。
そうだ。なんか思い出した。そういえば,私が育った昭和時代の百姓の世界は,とても単純で複雑だった。なにか道具がほしければ,自分で作ったもんだった。なにか問題があれば,自分の工夫で解決したもんだった。さすがにキーボードは無理にしても,リストレストくらいだったら,誰でも簡単に自分で好きなように好きなだけ作れる。そうだ,仕様が気に入らないなら,仕様を責める前に自分で作ればいい。いつから俺は仕様を責めて自分ではなにも作らない人間になってしまったのか。
…というわけで,息抜きにキーボードのリストレストを作ってみることにした。
まずは家の近くのホームセンターにいって道具を集める。基本的にそれっぽい寸法の木材を買えばいい。自分の手には横30cm,奥行き8cm,高さは適当でいい,といって決めた。大概のホームセンターでは寸法を頼めば安く切ってくれるサービスがある。で,集めた材料とかはこんな感じ。
- 木材 30cm×8cm×5cm ×3 300円位
- 黒色スプレー
- やすり(持っている)
- サンドペーパー(100円)
- 木工用ボンド(100円)使わなかった
- すべり止めシート×3 300円
で,適当やすりで角を取って,適当にスプレーで黒く塗ると,どん。
おお,いいじゃないの。
で,裏にすべり止めシートを張ると,どん。
完成。20分くらい。技術も知識も手間もなにもいらない。
やば,高さよく考えてなかったから高すぎて到底使えないwwwww
高すぎるwwwww
あっれー!やっちまたー!
…次は高さもよく考えて作ろう。そうやって,トライアンドエラーで自分の好きな仕様を探っていこう。自分の手で作りながら。それって多分,とても単純で,そして心を捉える複雑なこと。
(続く)