線形回帰と一般逆行列を比較する

• NxM行列Xに係数ベクトルRをかけてPが出るような関係があるとする
• 
• Xの一般逆行列(こちらやこちら)を出すことで、XとPからRを逆算することを考える
• 今、のように誤差項も考えて、そのときにもRを逆算することを考える
• 同じようなことを一般線形回帰でもできる
• 両者を比較する
  • 誤差がないときは、両方ともきれいに推定。誤差が入ると、両者に差が見えてくる

N<-50 # サンプル数
M<-40 # 項目数
# 適当にデータを作ろう
X<-matrix(sample(0:2,N*M,replace=TRUE),N,M)

OriP<-X%*%R
OriV<-var(OriP)
# 一般化逆行列
	G <- ginv(X)
	X %*% G %*% X # 条件(i)の確認
	G %*% X %*% G # 条件(ii)の確認
	t(G %*% X) # 条件(iii)の確認
	t(X %*% G) # 条件(iv)の確認
sum((G%*%OriP-R)^2)
H<-seq(from=0,to=1,length=11)
Diffs<-rep(0,length(H))
Diffs2<-Diffs
for(i in 1:length(H)){
	if(H[i]==0){
		P<-OriP
	}else{
		P<-OriP+rnorm(N,sd=sqrt(H[i]*OriV))
	}
	Diffs[i]<-sum((G%*%P-R)^2)
	Diffs2[i]<-sum((lm(P~X)[[1]][-1]-R)^2)
}
plot(H,Diffs)
plot(Diffs,Diffs2)