Hatena::ブログ(Diary)

arupaka-_-arupakaの日記

2018-12-16

rstan で回帰 連続変数のカテゴリ化バージョン

rstan で回帰 連続変数のカテゴリ化バージョン

data {
  int N;
  int K;
  int C;
  // int G;
  matrix[N,K] X;
  int y[N];
  int start_flg[K];
  int order_flg[K];
  
  //int group[N];
  int Level_no;
  //number of room
  int m[N];
  
}

parameters{
  
  real a;
  vector<lower=-5,upper=5>[K-C] b;
  //real<lower=0> sigma;
  real<lower=0,upper=0.2> sigma_f;
  //real<lower=0,upper=0.5> lambda_f;
  real<lower=0> sigma3;
  vector<lower=0,upper=5>[C] sigma_c;
  vector<lower=0,upper=30>[N] rrr;
  real<lower=0.0,upper=1.0> q;
  //real<lower=0,upper=1>[N] p0;
  vector[N] f;
  //vector[G] f_c;
}

transformed parameters{
  
  vector[N] mu;
  vector[K] c;
  vector<lower=0>[K-C] sigma_b;
  real p0;
  real<lower=0,upper=1> pp0b;
  vector<lower=0,upper=1>[N] pp0;
  //matrix<lower=0,upper=1>[Level_no,N] theta;
  simplex[Level_no] theta[N];
  
  for(i in 1:K){
    if(order_flg[i]==0){
      
      if(start_flg[i]!=0){
        c[i]=0;
      }else{
        c[i]=b[i-sum(start_flg[1:i])];
        sigma_b[i-sum(start_flg[1:i])]=sigma_c[sum(start_flg[1:i])];
      }
    }else{
      
      if(start_flg[i]!=0){
        c[i]=0;
      }else{
        c[i]=c[i-1]+b[i-sum(start_flg[1:i])];
        sigma_b[i-sum(start_flg[1:i])]=sigma_c[sum(start_flg[1:i])];
      }
      
    }
  }
  //c[K]=-sum(b[1:(K-1)]);
  //for(i in 1:N){
    //  f[i]=f_c[group[i]];
    //
      //}
  
  
  mu=X*c+a+1*f;
  
  
  for(i in 1:N){
    //print(i);
    //print(f[i]);
    //print(mu[i]);
    //pp0[i]=inv_logit(mu[i]);
    pp0[i]=1.0/(1+exp(-mu[i]));
    pp0b=1-(1-pp0[i])^rrr[i];
    //pp0b=1.0/(1+exp(-mu[i]));
    
    //print("pp0b");
    //print(pp0b);
    p0=binomial_lpmf(1|m[i],pp0b);
    // print("exponential");
    //print(exp(p0));
    //print("rr");
    //print(rrr[i]);
    //print("q");
    //print(q);
    //theta[i][1]=q+(1-q)*(1-(1-exp(p0))^(rrr[i]));
    theta[i][1]=q+(1-q)*(exp(p0)/(1-exp(binomial_lpmf(0|m[i],pp0b))));
    //print("theta");
    // print(theta[i][j]);
    //pp0[i]=inv_logit(mu[i]);
    
  }
  
  
  for(i in 1:N){
    for(j in 2:Level_no){
      if(j<=m[i]){
        //pp0[i]=inv_logit(mu[i]);
        pp0[i]=1.0/(1+exp(-mu[i]));
        //pp0b=pp0[i];
        pp0b=1-(1-pp0[i])^rrr[i];
        //  print("pp0b");
        //print(pp0b);
        p0=binomial_lpmf(j|m[i],pp0b);
        //print("ij");
        //print(i);
        //print(j);
        //print("mm");
        //print(m[i]);
        //print("rr")
        //print(rrr[i]);
        //print("p0");
        //print(exp(p0));
        //if(p0>=-9){
          theta[i][j]=(1-q)*(exp(p0)/(1-exp(binomial_lpmf(0|m[i],pp0b))));
          //}else{
            //   print("PPP000");
            //   theta[i][j]=(1-q)*0.0001;
            
            //}
      }else{
        theta[i][j]=0;
        
      }
      //print("theta");
      // print(theta[j,i]);
      
    }
  }
  
}

model{
  sigma3 ~ student_t(4,0,3);
  sigma_c ~ normal(0,sigma3);
  rrr ~ gamma(5.6,1.0);
  //rrr ~ normal(4,0.00001);
  //sigma_f ~ exponential(lambda_f);
  b ~ double_exponential(0,sigma_b);
  //b ~ double_exponential(0,1);
  //b ~ double_exponential(0,0.3);
  //b ~ student_t(4,0,0.5);
  //sigma_f ~ student_t(4,0,0.3);
  sigma_f ~ student_t(4,0,3);
  
  //for(i in 1:(K-1)){
    //b ~ normal(0,sigma2);
    //}
  //b ~ rnorm(0,sigma2);
  //f_c ~ normal(0,sigma_f);
  f ~ normal(0,sigma_f);
  //y ~ neg_binomial_2_log(mu, 1.0);
  for(i in 1:N){
    //print("smaple");
    //print(theta[i]);
    y[i] ~ categorical(theta[i]);
    //print(y[i]);
  }
}

2018-11-11

多変量回帰のチェック

多変量回帰で共通した情報を持つ場合は、

共通した情報程度の係数になる。 ただし、AIC Step wise などを使うとノイズによって

不安定になり、片方に集中することも。。 Ridge回帰なら両方の変数を選んでくる。

例えば、

a1=情報1+0.2情報2

a2=情報1+0.3情報3

で、y=情報1+情報2

なら、a1の回帰係数のが少し大きくなる。


もし、片方に回帰係数が寄った場合、ランダム誤差か、それとも、情報1の割合が

大きく、情報3の割合がちいさかったからか? yへの

情報もってるからといって、回帰係数がおおいとは限らない。

回帰係数が真に大きい潜在変数に近いほうが回帰係数が大きくなる。

b1が増えて、b2=0が増えたとき、a1とa2がどれくらい増えるか連立方程式で計算。それをもとに係数が決められる。

library(glmnet)
a1<-rexp(1000)
a2<-rexp(1000)
a3<-rexp(1000)
b1<-0.9*a1+0.05*a2
b2<-0.9*a1+0.06*a2

y<-3*a1+0.1*a2+0.0*a3+1*rnorm(1000,sd=1)

ans<-lm(y~b1+b2)
coef(ans)
step(ans)

ans0<-lm(y~1)

step(ans0,direction=c("forward"),scope=list(upper="y~b1+b2",lower="y~1"))


#Lasso
fit.glmnet.lasso.cv <- cv.glmnet(as.matrix(cbind(b1,b2)),y,alpha = 1)
coefficients(fit.glmnet.lasso.cv)

#Ridge
fit.glmnet.lasso.cv <- cv.glmnet(as.matrix(cbind(b1,b2)),y,alpha = 0)
coefficients(fit.glmnet.lasso.cv)


a1<-rexp(1000)
a2<-rexp(1000)
a3<-rexp(1000)
b1<-0.9*a1+0.03*a2
b2<-0.9*a1+0.03*a3

y<-3*a1+0.1*a2+0.1*a3+1*rnorm(1000,sd=1)

ans<-lm(y~b1+b2)
coef(ans)
step(ans)

ans0<-lm(y~1)

step(ans0,direction=c("forward"),scope=list(upper="y~b1+b2",lower="y~1"))


#Lasso
fit.glmnet.lasso.cv <- cv.glmnet(as.matrix(cbind(b1,b2)),y,alpha = 1)
coefficients(fit.glmnet.lasso.cv)

#Ridge
fit.glmnet.lasso.cv <- cv.glmnet(as.matrix(cbind(b1,b2)),y,alpha = 0)
coefficients(fit.glmnet.lasso.cv)




a1<-rexp(1000)
a2<-rexp(1000)
a3<-rexp(1000)
b1<-0.9*a1+0.5*a2
b2<-0.9*a1+0.03*a3

y<-3*a1+0.1*a2+0.1*a3+1*rnorm(1000,sd=1)

ans<-lm(y~b1+b2)
coef(ans)
step(ans)

ans0<-lm(y~1)

step(ans0,direction=c("forward"),scope=list(upper="y~b1+b2",lower="y~1"))


#Lasso
fit.glmnet.lasso.cv <- cv.glmnet(as.matrix(cbind(b1,b2)),y,alpha = 1)
coefficients(fit.glmnet.lasso.cv)

#Ridge
fit.glmnet.lasso.cv <- cv.glmnet(as.matrix(cbind(b1,b2)),y,alpha = 0)
coefficients(fit.glmnet.lasso.cv)

3 x 1 sparse Matrix of class "dgCMatrix"
                     1
(Intercept) -0.1048879
b1           0.7910052
b2           2.2927133

2018-10-07

RStan でカテゴリカル変数の単回帰2

data {
  int N;
  int K;
  matrix[N,K] X;
  vector[N] y;
}

parameters{
  
  real a;
  vector[K] b;
  real<lower=0> sigma;
  //vector<lower=0>[K] sigma2;
  real<lower=0>[K] sigma2;
}

transformed parameters{
  vector[N] mu;
  mu=X*b+a;
 // mu=X*b;
  
  
  
}

model{
  //b ~ normal(0,sigma2);
  b ~ double_exponential(0,sigma2);
  //b ~ double_exponential(0,sigma2);
  y ~ normal(mu,sigma); 
}

 x<-tmp3[,5]
X<-dummy(x)
y<-as.numeric(tmp3[,1])
tmp<-list(x=X,y=y,N=length(x),K=length(table(x)))
 fit<-stan(file="test0_lm2.stan",data=tmp,seed=1235,chains=2)

ms<-extract(fit)

 plot(factor(names(X[1,])),apply(ms$b,2,mean),ylim=c(-30,30))
 points(factor(names(X[1,])),apply(ms$b,2,function(x){quantile(x,0.9)}),col=2)
 points(factor(names(X[1,])),apply(ms$b,2,function(x){quantile(x,0.1)}),col=2)
 names(X[1,])
  abline(h=0,lty=2)
data.frame(factor(names(X[1,])),apply(ms$b,2,mean),apply(ms$b,2,function(x){quantile(x,0.1)}),apply(ms$b,2,function(x){quantile(x,0.9)}))

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-07-24

RStan でカテゴリカル変数の単回帰

係数の和=0の制約をいれる場合。

##Fit with categorical variable 2
x0<-sample(1:4,100,replace=T)
library(dummies)
x<-dummy(x0)
x1<-x[,1]
x2<-x[,2]
x3<-x[,3]
x4<-x[,4]

a1<-2
a2<--2
a3<-10
a4<--10

y<-a1*x1+a2*x2+a3*x3+a4*x4+rnorm(100,sd=0.03)

#library("rstan")

data<-list(y=y,x2=x2,x1=x1,x3=x3,x4=x4,N=100,K=4)
fit<-stan(file="lm4.stan",data=data,seed=1234)
summary(fit)$summary


Rstan lm4.stan

data{
  
  int N;
  int K;
  real y[N];
  real x1[N];
  real x2[N];
  real x3[N];
  real x4[N];

}


parameters {
  
  real a[K-1];
  real b;
  real <lower=0> sigma;
  real <lower=0,upper=0.1> sigma_a;
}

transformed parameters{
  real y_base[N];
  real aa[K];
  real sum1=0;


  for(i in 1:(K-1)){
    aa[i]=a[i];
    sum1=sum1+aa[i];
  }
  aa[K]=0-sum1;

  for(i in 1:N){
   
   
    y_base[i]=aa[1]*x1[i]+aa[2]*x2[i]+aa[3]*x3[i]+aa[4]*x4[i]+b;
    
    
  }
  
  
}

model {
  
  for(i in 1:N){
    y[i] ~ normal(y_base[i],sigma);
  }  
}


カテゴリが0の制約を入れる場合

a1x1+a2x2+a3x4+a4x4=0;

##Fit with categorical variable 2
x0<-sample(1:4,100,replace=T)
library(dummies)
x<-dummy(x0)
x1<-x[,1]
x2<-x[,2]
x3<-x[,3]
x4<-x[,4]

a1<-2
a2<--2
a3<-10
a4<--10

y<-a1*x1+a2*x2+a3*x3+a4*x4+rnorm(100,sd=0.03)


data<-list(y=y,x2=x2,x1=x1,x3=x3,x4=x4,N=100,K=4)
fit<-stan(file="lm5.stan",data=data,seed=1234)
summary(fit)$summary

ms<- rstan::extract(fit)

apply(ms$aa2,2,mean)

lm5.stan

data{
  
  int N;
  int K;
  real y[N];
  real x1[N];
  real x2[N];
  real x3[N];
  real x4[N];

}


parameters {
  
  real a[K-1];
  real b;
  real <lower=0> sigma;
  real <lower=0,upper=0.1> sigma_a;
}

transformed parameters{
  real y_base[N];
  real aa[K];
  real sum1=0;


  for(i in 1:(K-1)){
    aa[i]=a[i];
  //  sum1=sum1+aa[i];
  }
  aa[K]= -(aa[1]*sum(x1)+aa[2]*sum(x2)+aa[3]*sum(x3))/sum(x4);

  for(i in 1:N){
   
   
    y_base[i]=aa[1]*x1[i]+aa[2]*x2[i]+aa[3]*x3[i]+aa[4]*x4[i]+b;
    
    
  }
  
  
}

model {
  
  for(i in 1:N){
    y[i] ~ normal(y_base[i],sigma);
  }  
}