Posted at
RDay 23

最尤法 ~LogLogisticDistributionを例に~

More than 1 year has passed since last update.


前置き:Log Logistic Distributionについて

log Logistic Distributionという確率分布があります。(以下loglogistic分布と略)



赤色:$\alpha$=1 , $\beta$=1

黄色:$\alpha$=2 , $\beta$=1

空色:$\alpha$=1 , $\beta$=2

青色:$\alpha$=3 , $\beta$=2

緑色:$\alpha$=3 , $\beta$=5

「非負の二つのパラメータ」$\alpha$ , $\beta$からなる上図のような裾の重い確率分布で、その確率密度関数を$f(x;\alpha,\beta)$とすると

$$

f(x;\alpha,\beta)=

\frac{\beta(\frac{x}{\alpha})^{\beta-1}}

{\alpha[1+(\frac{x}{\alpha})^{\beta}]^{2}}

$$

またn個の観測値がこのloglogistic分布から得られたものとして、その対数尤度をlとすると、

$$

l=n(log(\beta)-log(\alpha))+\sum_{i=1}^n [(\beta-1)log(\frac{x_i}{\alpha})-2\log(1+(\frac{x_i}{\alpha})^\beta)]

$$

が得られます。


本題:riversデータとLoglogistic分布の最尤推定

Rに標準実装されているデータセットにriversというものがあります。

北米にある141の主要な河川の長さ(単位:mile)のデータで、density()で密度推定してみると下図のようになっています。

plot( density(rivers) )

このriversデータに対してloglogistic分布を当てはめ最尤推定をしてみます。

但し重要な点が一つあります。

通常であればoptim関数を使って対数尤度を最大化すればそれでよいのですが、logligistic分布の場合はそうはいかない。

というのも、二つのパラメータα,βはともに非負であるという制約があるので単純にoptim関数で最適化するには不都合なわけです。

そこで非負制約下の最適化問題(不等式制約下の最適化問題)に対応できるconstrOptimの出番。


#対数尤度#
loglikelihood<- function(x){
return(
function(parameter){
alpha <- parameter[1]
beta <- parameter[2]
n<-length(x)
return( n*(log(beta)-log(alpha))+sum( (beta-1)*log(x/alpha)-2*log(1+(x/alpha)^beta) ) )
}
)
}

#制約条件:α>0,β>0に対応#
ui <- matrix(c(1, 0, 0, 1),ncol=2,nrow=2)
ci <- c(0, 0)

#制約付き最適化#
constrOptim(c(1, 1), loglikelihood(rivers), grad=NULL, ui=ui, ci=ci, control=list(fnscale=-1))



さて得られた推定結果を基に曲線を当てはめて見ましょう。

parameter<-constrOptim(c(1, 1), loglikelihood(rivers), grad=NULL, ui=ui, ci=ci, control=list(fnscale=-1))$par

x<-seq(0,4000,1)
plot( density(rivers) )
lines( dllogis(x,parameter[2],parameter[1]),col="red" )



※赤が最尤推定に基づくLogLogistic分布

まずまず良い結果が得られたのではないでしょうか。


オマケ:loglogistic分布をstanに実装してみる

loglogistic分布をstanに実装してみます。

※rstan:Version 2.16.2

library(rstan)

y<-rivers
data<-list(y=y,N=length(y))

log_logistic_code <- "
functions{
real log_logistic_lpdf(real y, real alpha, real beta){
return( log(beta)+(beta-1)*log(y/alpha)-log(alpha)-2*log(1+(y/alpha)^beta) );
}
}

data {
int<lower=0> N;
real y[N];
}
parameters {
real<lower=0> alpha;
real<lower=0> beta;
}
model {
for(t in 1:N){
y[t] ~ log_logistic(alpha,beta);
}
}
"
log_logistic<-stan(model_code=log_logistic_code,data=data,seed=3141592,chains=4,warmup=1000,iter=3500)
print(log_logistic)



最尤推定値とおおむね一致しますね。


参照

wikipedia:Log-logistic distribution

Package 'eha' - CRAN.R-project.org