2018年4月9日(月)の統計解析Iという科目、今年はスパース推定の話をするということをシラバスに書いた。全体で学生90名(M1が75名程度)も聞きにきて、教室は立ち見になった。初回はガイダンスというのが普通のようだが、初回から演習課題24問を課した。厳しすぎるといって逃げ出す学生もいるかも知れないが、「手を動かさないと身につかない」という信念をわかってもらえると信じている。
- 説明変数と目的変数の$N$個の組$(x_{1},y_1),\cdots,(x_{N},y_N)$から、$S:=\sum_{i=1}^N(y_i-\beta_0-\beta_1x_{i})^2$を最小にする切片$\beta_0$と傾き$\beta_1$をもとめたい。
$\displaystyle \overline{x}=\frac{1}{N}\sum_{i=1}^Nx_{i}$, $\displaystyle \overline{x^2}=\frac{1}{N}\sum_{i=1}^Nx_{i}^2$, $\displaystyle \overline{y}=\frac{1}{N}\sum_{i=1}^Ny_i$, $\displaystyle \overline{xy}=\frac{1}{N}\sum_{i=1}^Nx_{i}y_i$,
を用いて、それらの解$\hat{\beta}_0$, $\hat{\beta}_1$を求めよ。 - 前問で、各$i=1,\cdots,n$に対して、最初に$x_i$から$\overline{x}$を、$y_i$から$\overline{y}$をひいて、新しい$x_i,y_i$に対して計算した$\overline{x}$, $\overline{y}$が0になるようにしてから$\hat{\beta}_1$をもとめるにはどうすればよいか。そして、求められた$\hat{\beta}_1$と、0にする前のもとの$\overline{x}$, $\overline{y}$とを用いて、$\hat{\beta}_0$はどのようにかけるか。以下では、$p$個の説明変数と目的変数の$N$個の組$(x_{1,1},\cdots,x_{1,p},y_1),\cdots,(x_{N,1},\cdots,x_{N,p},y_N)$から、切片$\beta_0$と各変数の傾き$\beta=(\beta_1,\cdots,\beta_p)$を推定したい。各$x_{i,j}$から$\displaystyle \bar{x}_j=\frac{1}{N}\sum_{i=1}^Nx_{i,j}$, $j=1,\cdots,p$をひき、各$y_i$から$\displaystyle \bar{y}=\frac{1}{N}\sum_{i=1}^Ny_i$をひき、それぞれの平均を0にしてから$\beta$を推定する(推定値を$\hat{\beta}$とする)。最後に$\hat{\beta}_0=\bar{y}-\sum_{j=1}^p \bar{x}_j\hat{\beta}_j$とする。また、$X=(x_{i,j})_{1\leq i\leq N, 1\leq j\leq p}\in {\mathbb R}^{N\times p}$, $y=(y_i)_{1\leq i\leq N}\in {\mathbb R}^N$とおくものとする。
- $X^TX$が正則のときに、
$$||y-X\beta||^2:=\sum_{i=1}^N(y_i-\beta_1x_{i,1}-\cdots-\beta_px_{ip})^2$$
を最小にする$\beta$が$\hat{\beta}=(X^TX)^{-1}X^Ty$であたえられることを示せ。 - 行列$X\in {\mathbb R}^{n\times p}$とベクトル$y\in {\mathbb R}$を入力し、切片の推定値$\hat{\beta}$と傾きの推定値$\hat{\beta}\in {\mathbb R}^{p}$の推定値を出力する関数linearをR言語で構成したい。以下の(1)(2)をうめよ。
inner.prod=function(x,y)sum(x*y)
linear = function(X,y){
n=nrow(X); p=ncol(X);
X=as.matrix(X); x.bar=array(dim=p); for(j in 1:p)x.bar[j]=mean(X[,j]);
for(i in 1:p)X[,j]=X[,j]-mean(X[,j]); ## Xの平均化
y=as.vector(y); y.bar=mean(y); y=y-y.bar ## yの平均化
beta=as.vector(## ここをうめる(1)
beta.0=## ここをうめる(2)
return(list(beta=beta,beta.0=beta.0))
} - $\lambda$を非負定数として、
$$\frac{1}{2N}||y-X\beta||^2+\frac{\lambda}{2}||\beta||_2^2:=\frac{1}{2N}\sum_{i=1}^N(y_i-\beta_1x_{i,1}-\cdots-\beta_px_{ip})^2+\frac{\lambda}{2}\sum_{j=1}^p\beta_j^2$$
を最小にする$\beta$をもとめよ(Ridge回帰)。 - 2.の関数にさらに、入力としてさらに非負定数$\lambda$を加えて、3.の処理を行うR言語の関数ridgeを求めよ。
- https://web.stanford.edu/~hastie/StatLearnSparsity/data.htmlの米国犯罪データをダウンロード(その画面でCtrl-a, Ctrl-cでコピーして、各自PCのエディタに新規ファイルオープンして、Ctrl-vでペースト)して、crime.txtというファイルに格納し、2.,4.で作成したプログラムを実行せよ。
列 説明/目的 変数の意味
1 目的 人口100万人あたりの犯罪率
2 (今回は用いない)
3 説明 警察官の年間給与
4 説明 25歳以上で高校を卒業した人の割合
5 説明 16-19歳で高校に通っていない人の割合
6 説明 18-24歳で大学生の割合
7 説明 25歳以上で4年制大学を卒業した人の割合
Rを実行する際には、そのファイルが格納されているディレクトリにアクセスできるようにする。
crime=read.table(“crime.txt”)
> X=crime[,3:7]
> y=crime[,1]
> linear(X,y)
$beta
[1] 10.8637397 -7.6530479 3.6495085 -0.5302667 8.0399178
$beta.0
[1] 606.2853
> ridge(X,y,100)
$beta
[1] 6.5983750 -3.7691561 1.0359781 -1.2917819 0.8428258
$beta.0
[1] 701.4541
8. $p=2$のとき、$X^TX$の各成分を、$v_1:=\sqrt{\frac{1}{N}\sum_{i=1}^Nx_{i,1}^2}$,
$v_2:=\sqrt{\frac{1}{N}\sum_{i=1}^Nx_{i,2}^2}$, $\displaystyle \rho:=\frac{\frac{1}{N}\sum_{i=1}^Nx_{i,1}x_{i,2}}{v_1v_2}$を用いて表示せよ。
9. $X^TX$の固有値が$\gamma_1,\cdots,\gamma_p$のとき、$X^TX+N\lambda I$の固有値が$\gamma_1+N\lambda,\cdots,\gamma_p+N\lambda$となることを示せ。また、$X^TX$に逆行列が存在しない条件を$\gamma_1,\cdots,\gamma_p$を用いてあらわせ。さらに、$\lambda>0$である限り、$X^TX+N\lambda I$には逆行列が必ず存在することを示せ。ただし、非負定値($X^TX$のように、ある行列$A$があって$A^TA$とかける行列として定義される)の固有値がすべて非負となることは証明しないで用いてよい。
10. 関数$f: {\mathbb R}^p\rightarrow {\mathbb R}$が、任意の$0<\alpha<1$と$u_1,u_2\in {\mathbb R}^p$について
$$f(\alpha u_1+(1-\alpha)u_2)\leq \alpha f(u_1)+(1-\alpha) f(u_2)$$
であるとき、関数$f$は凸であるという。この定義に基づいて、下記の各関数が凸であることを示せ。
(a) $f(u)=u^2$
(b) $f(u,v)=u^2+v^2$
(c) $f(u)=|u|$
11. $\displaystyle f(\beta_1,\cdots,\beta_p):=\sum_{i=1}^N(y_i-\sum_{j=1}^px_{i,j}\beta_j)^2$が凸であることを示せ。ただし、凸関数の和が凸関数になることは、証明無しで用いてよい。
12. 以下の各関数$f: {\mathbb R}\rightarrow {\mathbb R}$と$u_0\in {\mathbb R}$について、
$$f(u)\geq f(u_0)+z(u-u_0)\ ,\ u\in {\mathbb R}$$
が成立する$z\in {\mathbb R}$の集合、および$f(u)$の$u=u_0$での微分係数$f'(u_0)$を求めよ。今後、この集合を$\partial f(u)$とかくものとする(劣勾配)。また、要素が1個しかない場合、集合ではなく($\{$,$\}$をつけず)要素だけを書く場合がある。
(a) $f(u)=u$, $u_0=0$
(b) $f(u)=u^2+2u$, $u_0=1$
13. $f(u)=|u|$について
$$\partial f(u_0)=
\left\{
\begin{array}{ll}
-1,&u_0<0\\
{[-1,1]},&u_0=0\\
1,&u_0>0\\
\end{array}
\right.
$$
となることを示せ。
14. $\lambda$を非負定数として
$$f(u)=\frac{1}{2N}\sum_{i=1}^N(y_i-ux_{i,1})^2+{\lambda}|u|$$
について
$$\partial f(u_0)=
\left\{
\begin{array}{ll}
-\frac{1}{N}\sum_{i=1}^Nx_{i,1}(y_i-u_0x_{i,1})-\lambda,&u_0<0\\
-\frac{1}{N}\sum_{i=1}^Nx_{i,1}y_i+\lambda[-1,1],&u_0=0\\
-\frac{1}{N}\sum_{i=1}^Nx_{i,1}(y_i-u_0x_{i,1})+\lambda,&u_0>0\\
\end{array}
\right.
$$
となることを示せ。また、$0\in \partial f(0)$となる条件を求めよ。ただし、劣勾配の場合でも微分と同様、
関数の和の劣勾配がそれぞれの劣勾配の和になること、関数の定数倍の劣勾配がもとの関数の劣勾配の定数倍になることは、証明無しで用いて良い。
15. $p=1$, $\frac{1}{N}\sum_{i=1}^Nx_{i,1}^2=1$のとき、$\lambda$を非負定数として、
$$\frac{1}{2N}||y-X\beta||^2+{\lambda}||\beta||_1:=
\frac{1}{2N}\sum_{i=1}^N(y_i-\beta_1x_{i,1}-\cdots-\beta_px_{ip})^2+{\lambda}\sum_{j=1}^p|\beta_j|$$
を最小にする$\beta=(\beta_1)$をもとめよ(Lasso回帰)。
16.
$\lambda$を非負定数として、関数
$$S_{\lambda}(x):=
\left\{
\begin{array}{ll}
x-\lambda,&x>\lambda\\
0,&|x|\leq \lambda\\
x+\lambda,&x<-\lambda\\
\end{array}
\right.
$$
が、関数$(x)_+=\max\{x,0\}$,
$${\rm sign}(x)=\left\{
\begin{array}{ll}
-1,&x<0\\
0,&x=0\\
1,&x>0
\end{array}
\right.
$$
を用いて、$S_\lambda(x)={\rm sign}(x)(|x|-\lambda)_+$とかけることを示せ。
17. $\lambda>0,x\in {\mathbb R}$を入力として、$S_{\lambda}(x)$を出力する関数soft.thをR言語でかけ。
さらに、下記を実行し、関数の動作が正しいことを確認せよ(出力をpdfで提出する)。
curve(soft.th(5,x),-10,10)
場合分けで関数を定義したり、maxを用いると、R言語ではグラフを描いてくれない。sign, abs, pmaxを用いよ。
18. 12.の解を関数$S_\lambda(\cdot)$を用いてかけ。また、R言語で記述するとき、下記の空欄をうめよ。
linear.lasso.1=function(x,y,lambda){
x.bar=mean(x); y.bar=mean(y);
x.scale=scale(x) # 平均化するだけでなく、分散が1になるようにする
beta=## ここをうめる
beta.0=y.bar-beta*x.bar
return(list(beta=beta,beta.0=beta.0)
}
19. 下記は、$j$番目の説明変数を固定して、残差
$$y_i-\sum_{k\not=j}x_{i,k}\hat{\beta}_k\ ,\ i=1,\cdots,N$$
を計算して、係数$\hat{\beta}_j$を推定するということを、$j=1,\cdots,p$で繰り返し、
さらにそれら全体を繰り返して、$\hat{\beta}_1,\cdots,\hat{\beta}_p$の値の収束をまつ処理(Lasso推定の座標降下法)である。下記の空欄(残差を計算する箇所)をうめて、$\lambda=0$のときに、関数linearと近い値になることを確認せよ。
linear.lasso=function(X, y, lambda=0){
INF=1000
X=as.matrix(X); y=as.vector(y)
p=ncol(X)
X.bar=array(dim=p); for(j in 1:p){X.bar[j]=mean(X[,j]); X[,j]=X[,j]-X.bar[j];}
y.bar=mean(y); y=y-y.bar
beta=array(0, dim=p)
eps=1
beta.old=INF
while(eps>0.01){
for(j in 1:p){
r= ## ここをうめる
beta[j]=soft.th(lambda,cov(r, X[,j]))/cov(X[,j],X[,j])
}
eps=max((beta-beta.old)/beta.old)
beta.old=beta
}
beta.0=y.bar
for(j in 1:p)beta.0-x.bar[j]*beta[j]
return(list(beta=beta, beta.0=beta.0))
}
20. $\hat{y}=X\hat{\beta}$として、$X^T(y-\hat{y})=0$が成立すること、すなわち$\sum_ix_{i,j}(y_i-\hat{y_i})=0$, $j=1,\cdots,p$となることを示せ。次に、$p=2$, $c>\sum_{i=1}^n(y_i-\hat{y}_i)^2$として、$\beta=(\beta_1,\beta_2)$に関する方程式$||y-X\beta||^2=c$が$\beta=(\hat{\beta}_1,\hat{\beta}_2)$を中心とする楕円曲線になることを示せ。ただし、必要であれば、以下の等式を用いて良い。
$$||y-X\beta||^2=\sum_{i=1}^N
(y_i-\beta_1x_{i,1}-\beta_2x_{i,2})^2=\sum_{i=1}^N\{(y_i-\hat{y}_i)-(\beta_1-\hat{\beta}_1)x_{i,1}-(\beta_2-\hat{\beta}_2)x_{i,2}\}^2$$
21. 下記は、犯罪のデータについて、$\lambda$の値を変えながら、各変数の係数がどのように変化するかを表示したものである。plotとlinesの行の$\lambda$を$\log \lambda$にかえて、グラフを表示させよ。
df=read.table(“crime.txt”)
x=df[,3:7]
y=df[,1]
p=ncol(x)
lambda.seq=seq(0,2000,10)
plot(log(lambda.seq),coef.seq, ylim=c(-12,12), type=”n”, col=”red”) ##この行
for(j in 1:p){
coef.seq=NULL; for(lambda in lambda.seq)coef.seq=c(coef.seq,ridge(x,y,lambda)$beta[j])
par(new=TRUE)
lines(lambda.seq,coef.seq, col=j) ##この行
}
legend(“topright”,legend=
c(“警察官の年間給与”,”25歳以上で高校を卒業した人の割合”,
“16-19歳で高校に通っていない人の割合”,”18-24歳で大学生の割合”,
“25歳以上で4年制大学を卒業した人の割合”), col=1:p, lwd=2, cex =.8)
22. ライブラリglmnetをインストールして、linear.lassoをglmnetにおきかえて実行せよ。
23. 前問の関数ridgeをlinear.lassoにかえて実行し、グラフを表示させよ。特に、すべての変数の係数が0になるようにlambda.seqの範囲を調節せよ。
24. すべての変数の係数が0になるような$\lambda$の最小値が、$\displaystyle \lambda=\max_{1\leq j\leq p}\frac{\sum_{i=1}^Nx_{i,j}y_i}{\sum_{i=1}^Nx_{i,j}^2}$であたえられることを示せ。