草薙の研究ログ

英語の先生をやってます。

モデルの中で何が捨象できるかを語らない科学

数理モデルというものは,その記述の仕方の形式性の割には,数理モデルということばに親しみを感じないほとんどのひとが思うより,本来結果主義的で効用主義的なものだ。

モデルは,もちろん現象それ自体ではないし,その現象を大幅に捨象していて,しかしそこから得られる予測や知見が有益だと見込まれているものだ。この「有益だ」という考え方は一部の学術分野にはないこともあるが。

人文社会系のほとんどの数理モデルは,世界が数字に支配されていて,その世界の斉一的な決定則を表すものだなんてことを意味しない。せいぜいが「観測がうまく当てはまる」,「うまい数理的な近似になっている」という程度の含意である。しかし,そのモデル(世界の決定則それ自体ではない)について考えることで,人が適切に意思決定をできたり,そして個人間の合意が得られ,判断の公共性が発生する(たとえばある種のエビデンスになる)というに考えられている。

 

たとえば,TOEICのスコアは大学が自前でやっている単語テストの成績から予測できるとする。

y = ax + b

という簡単なモデルを考えて,

TOEIC = 8.5×単語テストの得点 - 120

とか,そんなふうに。

このモデルは,それがうまくTOEICのスコアを予測できるとか,つまりこの現象のいい数理的な近似になっていれば,まあいいモデルだといえる。

 

ただし,このとき,TOEICは単語の力だけではないのだから,このモデルは世界の斉一的な決定則の記述として不完全だ!というのはちょっと勇み足だと思う。

TOEICは単語のテストだけではない!とか。

もっというなら,TOEICのときの部屋の温度をモデルに取り入れるべきかだとか,世界の斉一的な決定則を考えれば無限にそんな要因はありえる。

数理モデルを扱う人がもつ,結果主義というか,もっと大きくいってプラグマティズム的な考え方の下では,単にこのモデルはそれらを捨象しているのである。線形モデルなら,それら全部をひっくるめて誤差(残差)として考える。

 

ここで問題は,何が捨象できるかそして何を積極的に捨象するべきかということになる。

 

…外国語教育研究やその関連分野である第二言語習得研究が扱ってきた変数や現象は莫大であり,確かに複雑極まりない。もはや還元主義の傾向が強いテーマに関しては指数的に用語が増えていっていて,収拾付かない例もある。

だからこそ,なにが捨象できるか,そしてなにを捨象するべきか,という考えが今後重要になってくると思う。

そして個人的には,大局的に人間の意思決定や判断の合意形成や公共性の創出という点から見れば,捨象できるものはずっとずっと多いと思っている。

…世界の斉一的な決定則それ自体を考えるならば,むしろなにも捨象できないとも思っている。

Rで日付データの処理

自分用のメモ。

#日付クラスへ変換
d<-"2016-1-1"
d2<-as.Date(d)
class(d2)

#日付データの足し引き(日付クラスだとこれができるようになるのが最高)
d2-1
d2+1
d2-1000
d2+1000

#基準日から1日毎にログイン回数を累積計算
#datは時間とログイン回数のデータフレーム

#基準日を設定
start<-as.Date("2016-04-08")
r<-numeric(130)

#計算
for(i in 1:130){
 r[i]<-sum(dat[dat[,1]<start+i,2])
 }

一般化パレート分布をデータに当てはめる

一般化パレート分布は所得の分布などに使われるそうだ。
外国語教育研究でもこういった分布になる変数を私はひとつだけ知っている(いわないwww)。

Rにいろいろあると思うけど,ここではactuarパッケージとfitdistrplusパッケージを使う。
actuarパッケージに関数があるから,それをfitdistrplusパッケージのfitdist関数当てはめるというわけ。

#パッケージの準備
library(actuar)
library(fitdistrplus)

#乱数作っちゃう
#第一形状母数が3,第二形状母数が3,尺度母数が1000
set.seed(1)
dat<-rgenpareto(1000,3,3,scale=1000)

#経験分布を可視化
par(mfrow=c(1,2))
plot(ecdf(dat),main="")
hist(dat,col="lightblue",main="")


f:id:kusanagik:20170131170612p:plain

#最尤推定(初期値は適当)
fit<-fitdist(dat,"genpareto",start=list(shape1=1,shape2=1,scale=1000))
summary(fit)

Fitting of the distribution ' genpareto ' by maximum likelihood 
Parameters : 
          estimate  Std. Error
shape1    2.973575   0.3122382
shape2    2.811662   0.2800328
scale  1070.970526 236.9332116
Loglikelihood:  -8230.236   AIC:  16466.47   BIC:  16481.2 
Correlation matrix:
           shape1     shape2      scale
shape1  1.0000000 -0.6723317  0.9114599
shape2 -0.6723317  1.0000000 -0.9020353
scale   0.9114599 -0.9020353  1.0000000

plot(fit)


f:id:kusanagik:20170131170939p:plain


うむ。

MCMCを使って指数正規合成分布(ex-Gaussian)の母数を推定

RのMCMCpackにはMCMCmetrop1Rっていう関数があって,これは任意(自作)の対数尤度の関数をいれてMCMCでサンプリングすることができる。なので,結構手軽にMCMCを使ってデータに好きな分布を当てはめることが可能。

ここでは(まったくそんなことはしなくてもいいんだけど),指数正規合成分布の母数(μ,σ,τ)をMCMCを使って推定してみる。

#準備
library(retimes)
library(MCMCpack)

#データの例
dat<-rexgauss(1000,3000,100,800)

#関数の準備
llf<-function(beta,x){
sum(log(dexgauss(x,beta[1],beta[2],beta[3])))
}

#MCMCしてみる,初期値は全部適当,burninとかmcmcの数とかデフォルトのまま
m<-MCMCmetrop1R(llf,theta.init=c(2000,100,100),x=dat)

#結果を見てみる
summary(m)

Iterations = 501:20500
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 20000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

       Mean    SD Naive SE Time-series SE
[1,] 3005.7 13.00  0.09189         0.3102
[2,]  110.8 10.92  0.07724         0.2642
[3,]  779.5 27.79  0.19647         0.6677

2. Quantiles for each variable:

        2.5%    25%    50%    75%  97.5%
var1 2979.96 2997.0 3005.7 3014.3 3031.0
var2   90.49  103.5  110.5  117.8  133.7
var3  726.38  760.6  778.9  798.1  836.6

plot(m)


f:id:kusanagik:20170128185732p:plain

非線形最小二乗法で学習コンテンツ消化曲線をモデル化

まあ結構いろんなことに汎用的に当てはまることだと思うのだけど,ソフトウェアの品質管理とかの分野では,ソフトウェア信頼度成長曲線という手法があるそうだ(古い友人に教えてもらった)。

これは,ソフトウェア開発において,バグの発見数や残ったバグの数の予測などに使われるものだそうで,基本的には横軸に時間や工数,縦軸にはバグの数とかをプロットする。すると,S字カーブのような曲線になるそうだ。で,これを非線形最小二乗法を使って何かの関数に当てはめてモデル化するというはなし。ふむふむ。とっても面白そう。品質管理というだけあって,実用的だし,外国語教育の業務改善にもすごく通じるものがありそうだ。

外国語教育では,オンライン教材におけるコンテンツ消化率とかがパッと思いつくような例だ。オンライン教材で単語テストを繰り返すとか,そういう場面では,横軸に時間を,縦軸に消化率をプロットすると,おそらく同じように成長曲線を描くことができる。こんな感じだということしよう。ある学生があるオンライン教材を15週にかけて勉強した様子で,横軸には週,縦軸にはそのコンテンツの消化率をプロットした,と。(こういうデータは外国語教育にはかなりお蔵入りしてる)

#例の作成
t<-seq(1,15,1)
d<-c(0.04,0.02,0.03,0.06,0.12,0.22,0.36,0.50,0.60,0.72,0.85,0.91,0.97,0.98,1.00)

#作図
plot(t,d,xlab="時間(単位:週)",ylab="コンテンツ消化率")

f:id:kusanagik:20170128150313p:plain

このデータになんらかの関数を当てはめたい,と(最近欲求不満なのかなんでも当てはめたくなる)。
えっと,いろんな関数が当てはまるだろうけど,ここでは,ゴンペルツと4母数ロジスティック関数にしよう。

Rでやるには,デフォルトで使えるstatsパッケージのnlsという関数を使えばいい。非線形最小二乗法でフィッティングをする関数だ。いろいろ便利。

#ゴンペルツに当てはめ
fit.gom<-nls(d~SSgompertz(t,Asym,b2,b3))

#結果
summary(fit.gom)

Nonlinear regression model
  model: d ~ SSgompertz(t, Asym, b2, b3)
   data: parent.frame()
   Asym      b2      b3 
 1.0985 12.1508  0.7114 
 residual sum-of-squares: 0.004886

Number of iterations to convergence: 0 
Achieved convergence tolerance: 5.022e-06

AIC(fit.gom)
BIC(fit.gom)

#当てはめた曲線を描き足し
plot(t,d,xlab="時間(単位:週)",ylab="コンテンツ消化率")
lines(t,predict(fit.gom),lty=2,col=2)

#4母数ロジスティック関数に当てはめ
fit.fpl<-nls(d~SSfpl(t,A,B,xmid,scal))

#結果
summary(fit.fpl)

Formula: d ~ SSfpl(t, A, B, xmid, scal)

Parameters:
      Estimate Std. Error t value Pr(>|t|)    
A    -0.009647   0.015871  -0.608    0.556    
B     1.023469   0.017314  59.114 4.01e-15 ***
xmid  8.223108   0.108640  75.692 2.66e-16 ***
scal  1.791406   0.112481  15.926 6.06e-09 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01943 on 11 degrees of freedom

Number of iterations to convergence: 0 
Achieved convergence tolerance: 3.818e-06

AIC(fit.fpl)
BIC(fit.fpl)

#当てはめた曲線をさらに描き足し
plot(t,d,xlab="時間(単位:週)",ylab="コンテンツ消化率")
lines(t,predict(fit.gom),lty=2,col=2)
lines(t,predict(fit.fpl),lty=2,col=4)

#2つのモデルをANOVA
anova(fit.gom,fit.fpl)

Analysis of Variance Table

Model 1: d ~ SSgompertz(t, Asym, b2, b3)
Model 2: d ~ SSfpl(t, A, B, xmid, scal)
  Res.Df Res.Sum Sq Df     Sum Sq F value Pr(>F)
1     12  0.0048859                             
2     11  0.0041547  1 0.00073119  1.9359 0.1916

#予測のために
#ゴンペルツ曲線の関数を定義
gomf<-function(x,Asym,b2,b3){
  y<-Asym*exp(-b2*b3^x)
  y
}

#4母数ロジスティック曲線の関数を定義
fplf<-function(x,A,B,xmid,scal){
  y<-A+(B-A)/(1+exp((xmid-x)/scal))
  y
}


この例だとピッタシだ。


f:id:kusanagik:20170128151403p:plain



外国語教育研究では,これと似たような技術の応用に,キーログを使った川口ほか(2016, Language Education & Technology)がある。

CiNii 論文 -  エッセイライティングにおける増加語数の時系列推移傾向とエッセイ評価の関係 : モデルフィッティングを用いた検討

混合正規指数合成分布モデル(?)を最尤法で…

聞いたこともないけど,要素数2の混合正規指数合成分布(ex-Gaussian)モデルというのを考えてみる。ま,2つの異なる認知プロセスが混合したときの反応時間の分布だとか,そんなそれっぽいことを考えてみる。そんなものは多分ない。

ま,でもこの確率密度関数は,

dmixexgauss<-function(x,lambda,mu1,sigma1,tau1,mu2,sigma2,tau2){
	y<-lambda*(1/tau1)*exp(mu1/tau1+(sigma1^2)/(2*tau1^2)-x/tau1)* pnorm(x,mu1+(1/tau1)*sigma1^2, sigma1)
+(1-lambda)*(1/tau2)*exp(mu2/tau2+(sigma2^2)/(2*tau2^2)-x/tau2)* pnorm(x,mu2+(1/tau2)*sigma2^2,sigma2)
return(y)
	}

こうなはず。これで尤度関数を作って,bbmleパッケージのmle2とかで最尤推定してみれば…と思ったんだけど,やっぱ上手く推定できなかった。methodもBFGSとかいろいろ試したけど,ううむ。なんと7母数の確率密度関数だもんな。俺の間違いってのもあるし。

ま,でもこの自作の確率密度関数自体はあってるぽい。暇ができたらちょっとやってみよう。

library(retimes)
set.seed(1)
dat<-c(rexgauss(400,1000,300,800),rexgauss(100,4000,800,500))
hist(dat,main="",breaks=20,col="lightblue",freq=F,ylim=c(0,.0006),xlim=c(0,8000))
x<-seq(0,8000,1)
lines(x,dmixexgauss(x,.80,1000,300,800,4000,800,500),lty=2,col=2)


f:id:kusanagik:20170127211001p:plain


*追記(1/28)
標本サイズ,母数,初期値次第では結構うまくいく場合もあるようだ。たとえば,

#あからさまな例
set.seed(1)
dat<-c(rexgauss(500,1000,500,100),rexgauss(500,10000,500,100))

#mle2で最尤推定
#対数尤度を返す関数
lf<-function(l,m1,s1,t1,m2,s2,t2){
-sum(log(dmixexgauss(dat,l,m1,s1,t1,m2,s2,t2)))
}

#適当な初期値で最尤推定
fit<-mle2(lf,start=list(l=.01,m1=1000,s1=100,t1=100,m2=10000,s2=100,t2=100))
summary(fit)

Maximum likelihood estimation

Call:
mle2(minuslogl = lf, start = list(l = 0.01, m1 = 1000, s1 = 100, 
    t1 = 100, m2 = 10000, s2 = 100, t2 = 100))

Coefficients:
     Estimate Std. Error  z value     Pr(z)    
l  5.0001e-01 1.5811e-02  31.6234 < 2.2e-16 ***
m1 8.8352e+02 5.3678e+01  16.4598 < 2.2e-16 ***
s1 4.2277e+02 2.8256e+01  14.9618 < 2.2e-16 ***
t1 2.5258e+02 5.1428e+01   4.9114 9.043e-07 ***
m2 9.8383e+03 5.6083e+01 175.4242 < 2.2e-16 ***
s2 4.4467e+02 2.8829e+01  15.4245 < 2.2e-16 ***
t2 2.4854e+02 5.3408e+01   4.6536 3.263e-06 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

-2 log L: 16646.24 

coef(fit)
vcov(fit)
AIC(fit)

当てはめた分布と推定した母数から平均,分散,歪度,尖度をもとめる

なんかある論文で,分布が強く歪んでいることが理論的に明白だった変数に,ガンマ分布か対数正規分布かなにかを当てはめて,その推定母数と適合度指標のみを(きつい紙幅の関係もあって)報告したときに,「標本の記述統計を報告しないとはけしからん」というようなことを査読者さまにご教示されたことがある。ま,ごもっともであるけど,紙幅の関係もあって,ううむ,標本自体がそんな重要なのかな…みたいに思いつつも,どうしたんだったかな。結構前の話だ。

でも,もちろん標本の記述統計量はわからんけど,この当てはめた分布の母数の推定値から,理論的には,そしてその分布モデルがしっかりと当てはまっているなら,標本のおおよその記述統計であれば,あとからでもとめることもできる。たとえばガンマ分布はかなり簡単で,こんな自作関数でも作ればいい。

gammadescriptive<-function(shape,scale){
	m<-shape*scale
	v<-shape*(scale^2)
	s<-2/sqrt(shape)
	k<-6/shape
	result<-list("mean"=m,"variance"=v,"skew"=s,"kurtosis"=k)
	result
}

たとえば,この関数に,形状母数3,尺度母数(scale)4を入れてやる。

gammadescriptive(3,4)

$mean
[1] 12

$variance
[1] 48

$skew
[1] 1.154701

$kurtosis
[1] 2

するとこんな感じになる。もちろん,「あくまでも理論値」なんだけど。

こんな感じで実験してみよう。

#パッケージの準備
library(psych)
library(fitdistrplus)

#例の作成
set.seed(1)
dat<-rgamma(1000,shape=3,scale=4)

#理論値
gammadescriptive(3,4)

$mean
[1] 12

$variance
[1] 48

$skew
[1] 1.154701

$kurtosis
[1] 2


#この例の記述統計を見てみる
describe(dat)
   vars    n  mean   sd median trimmed  mad  min  max range skew kurtosis   se
X1    1 1000 11.95 7.09   10.4   11.19 6.59 0.74 50.2 49.46 1.12     1.86 0.22

#大体近い
#また,この乱数にガンマ分布をフィットさせる
model<-fitdist(dat,"gamma")
model
Fitting of the distribution ' gamma ' by maximum likelihood 
Parameters:
       estimate Std. Error
shape 2.8585896 0.12111673
rate  0.2392977 0.01108254

#ちっ,これrateで出してきやがる,でもscale=1/rate
#でも,大体形状母数3,尺度母数4くらい

sh<-model$estimate[[1]]
sc<-1/model$estimate[[2]]

#またこの母数の理論値を見てみる
gammadescriptive(sh,sc)
$mean
[1] 11.94575

$variance
[1] 49.92004

$skew
[1] 1.182917

$kurtosis
[1] 2.098937


もちろんだけど,おおよそ標本の記述統計にも合っている。
だから,ある変数に狙った分布がよく当てはまっているといえるのであれば,平均,分散,歪度,尖度などは報告しなくても,「後からおおよその値は計算できる」こともあることは事実。ただ,この記述統計4つを報告したところで,この経験分布はガンマ分布に近いとか適合度はどれくらいだとか,そういうのはわからないか,または技術的により困難だってこと。

もちろん,標本の記述統計は要りません,って話じゃないけど。