Hatena::ブログ(Diary)

arupaka-_-arupakaの日記

2018-09-13

多変量回帰と相関

多変量解析では、

(0)単独変数の調整

 (0-1)カテゴリ化すると楽 cut関数

(0-2)少なすぎる変数の除去

 

(1)説明変数間の相関を調べる

  →とりあえず、主要っぽい変数を10個くらいとりだすとわかりやすい

http://kusanagi.hatenablog.jp/entry/2014/07/23/220951

(1-1)ライブラリ

     library(psych)
     pairs.panels(dat)

cor.plot(tmp2)

→各変数がどういう情報をもっているか把握

(1-2)因子分析や主成分分析等で整理

 (1-3) 交互作用の調査

(1-4) 相関係数、ファイ、相関比を使う

 (1−5)3変数の依存

i tmp2<-d2
     d<-data.frame(y=tmp2$y,x1=as.numeric(paste(tmp2$MENSEKI)),x2=as.numeric(paste(tmp2$TIKUNENSUU)),x3=as.numeric(paste(tmp2$TINRYOU)))
     
     d2<-d%>% dplyr::filter( x1<=300 & x2<=100) %>% dplyr::group_by(x1=20*round(x1/20),x2=2*round(x2/2)) %>% dplyr:: summarise(z=length(y[!is.na(y)]),y=median(y,na.rm=T)) %>% dplyr::arrange(x2,x1)
     
     d3<-data.frame(d2[-3]) %>% tidyr::spread(key=x1,value=y)
     x_name<-d3[,1]
     d3[d3>=200]<-200
     #d3[d3<=10]<-10
     #image(x=1:length(d3[[1]]), y=as.numeric(names(d3[-1])),z=as.matrix(d3[,-1]),col=terrain.colors(100))
     #image.plot(x=1:length(d3[[1]]), y=as.numeric(names(d3[-1])),z=as.matrix(d3[,-1]),col=terrain.colors(100),zlim=c(0,200))
     image.plot(x=1:length(d3[[1]]), y=as.numeric(names(d3[-1])),z=as.matrix(d3[,-1]),col=terrain.colors(100),zlim=c(0,200),xlab="Ages",ylab="Area (m^2)")
      mtext(at=1:length(d3[[1]]),side=1,text=x_name)
     text(x=as.numeric(factor(d2[[2]])),y=d2[[1]],labels=d2[[3]])
     #contour(x=as.numeric(rownames(d3[,1])),y=as.numeric(colnames(d3[1,])[-1]),as.matrix(d3),nlevel=10,add=T)
     contour(x=1:length(d3[[1]]),y=as.numeric(names(d3[-1])),z=as.matrix(d3[,-1]),nlevel=10,add=T)
     ###

pairs.panels(tmp[,1:10])

fa.parallel(tmp[,1:10])

fa(r=tmp2[-1,-1],nfactors=7)

cv.test(y, v1)

interaction.plot

(2)説明変数と被説明変数の関係

(2-1) とりあえず、すべてカテゴリ変数化すると非線形までいける

(2-2) glm を使う anova 等で説明力チェック

(2-3) stan 等でベイズ的に調整

(例1):隣同士は同じあたりb(t+1)=b(t)+e(t))

(例2): 事前分布→データ数少ない場合は平均でおきかえ

(2-4)カテゴリ変数は一つ0にしなきゃいけない。

(2-5)ラッソやグループラッソ。

(2-6)ダミー変数化にはdumiesパッケージが便利

2018-09-10

Rでgroup lasso の実験 一般化線形モデルの正則化

x1<-sample(1:100)
x2<-sample(1:100)
y<-2*x1+0.01*x2+rnorm(100)
p<-data.frame(y,x1,x2)
lm(y~x1+x2,data=p)

f<-function(r){
  sum((y-r[1]*x1-r[2]*x2)^2)
}

optim(c(1,3),f)

##Lasso
f2<-function(r){
  sum((y-r[1]*x1-r[2]*x2)^2)+20000*sum(abs(r))
}

optim(c(1,3),f2)$par

#Ridge
library(dummies)
x<-rbinom(100,size=1,prob=0.5)
y<-100*x+rnorm(100)
xx<-dummy(x)
p<-data.frame(y,x1)
lm(y~.,p)
x1<-xx[,1]
x2<-xx[,2]
f<-function(r){
  sum((y-r[1]*x1-r[2]*x-r[3])^2)+0000*sum(abs(r))
}

optim(c(1,3,1),f)$par

f2<-function(r){
  sum((y-r[1]*x1-r[2]*x2-r[3])^2)+100*sum(abs(r))
}
optim(c(1,3,1),f2)$par

f3<-function(r){
  sum((y-r[1]*x1-r[2]*x2-r[3])^2)+100*sum((r)^2)
}
optim(c(1,3,1),f3)$par

#group lasso
f4<-function(r){
  sum((y-r[1]*x1-r[2]*x2-r[3])^2)+100*((sum((r[1:2])^2))^0.5+(r[3]^2)^0.5)
}
optim(c(1,3,1),f4)$par

一般化線形モデル

正則化

―逸脱度/データ数+λ×正則化項?

によると、

http://bayes.sigmath.es.osaka-u.ac.jp/ftanaka/T/modeling/2016chap12_2.pdf

対数尤度+正則化項でいいっぽい。


逸脱度(二乗誤差みたいなやつ)=-2×対数尤度比?

https://mumu.jpn.ph/forest/computer/2017/08/20/11019/


###GLM check

x<-rbinom(1000,size=1,prob=0.5)
mu<-exp(log(100)*x+log(3))

y<-rgamma(1000,shape=mu,scale =0.5)
xx<-dummy(x)
#p<-data.frame(y,x1)
x1<-xx[,1]
x2<-xx[,2]
#group lasso
f5<-function(r){
    ll<-exp(r[1]*x1+r[2]*x2+r[3])
    l<- -log(dgamma(y,shape=ll,scale=(r[4])))
    l[l==Inf]<-10^8
    mean(l)+0.01*((sum((r[1:2])^2))^0.5+(r[3]^2)^0.5+(r[4]^2)^0.5)
}

ans<-optim(c(1,3,1,1),f5)$par
ans
#optim(c(1,3,1),f5)$par
y_th<-rgamma(1000,shape=exp(ans[1]*x1+ans[2]*x2+ans[3]),scale=ans[4])
qqplot(y_th,y)
plot(density(y))
points(density(y_th),col=2)


###Genraization###

2018-09-02

EXELEファイルのRの読み込み。着工件数の年次データをえる

特定の都道府県の賃貸物件の着工数をえる。

国交省の【住宅】利用関係別 構造別 建て方別 都道府県別 戸数 (Excel ファイル 3,225KB) 

ダウンロード

http://www.mlit.go.jp/sogoseisaku/jouhouka/sosei_jouhouka_tk4_000002.html

readxlパッケージの利用

http://anpontan382.hatenablog.com/entry/2017/06/22/070514

install.packages("readxl")
library(readxl)
tyakou_y<-NULL
tyakou_x<-NULL
for(i in 1:30){
  xls <- read_excel("../001218891.xls",sheet=i,range="C161:C161",col_names = FALSE,col_types = "text")
  xls2 <- read_excel("../001218891.xls",sheet=i,range="A3:A3",col_names = FALSE,col_types = "text")
  
  tyakou_y[i]<-as.vector(xls[[1]])
  tyakou_x[i]<-as.vector(xls2[[1]])
  cat(i,"\n")
}