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-11-05

Rstan を変な環境にインストール

その1:processxで詰まる。

(1)ソースをもってくる。

(2)ソースを展開して、一つ一つ下のようにビルド

(3)ソースを圧縮、Rに読ませる。install.packages("processx.tar",repos=NULL,type="source")


他は、processx-vector.c はg++に変えるとかで対処できる。

Rstanインストール時の

In file included from processx.h(34),

from create-time.c(4):

unix/processx-unix.h(29): error: identifier "siginfo_t" is undefined

void processx__sigchld_callback(int sig, siginfo_t *info, void *ctx);

^

compilation aborted for create-time.c (code 2)

とか、

In file included from processx.h:34:0,

from create-time.c:4:

unix/processx-unix.h:29:42: error: unknown type name ‘siginfo_t’

void processx__sigchld_callback(int sig, siginfo_t *info, void *ctx);

のエラーは、

gcc -std=gnu99 -I"/ap/r/3.5.1/lib64/R/include" -DNDEBUG -I/usr/local/include -D_POSIX_C_SOURCE=199309L -fpic -g -O2 -std=c99 -c create-time.c -o create-time.o

を加える必要あり。

-D_POSIX_C_SOURCE=199309L

https://github.com/geommer/yabar/issues/18

その2:rstanで詰まる

Error in .shlib_internal(args) :

C++14 standard requested but CXX14 is not defined

https://qiita.com/kouki111x1/items/0ae9ef4fff617a8c0f55

ホームの下に、.Rに以下のRmarkvarsように作成して対処

dotR <- file.path(Sys.getenv("HOME"), ".R")

if (!file.exists(dotR))

dir.create(dotR)

M <- file.path(dotR, "Makevars")

if (!file.exists(M))

file.create(M)

cat("\nCXX14FLAGS=-O3 -Wno-unused-variable -Wno-unused-function -fPIC",

"CXX14 = g++ -std=c++1y",

file = M, sep = "\n", append = TRUE)

途中で -fPICをいれろとあったので、追加

CXX14がはいってないので、下のバージョンで

2018-10-29

Rで %in%  リストの要素が存在するかのチェック

リストの要素が存在するかのチェックできる。便利。

y<-c(1,2,3)
 x<-c(1,2,2,3,4)

 x %in% y

出力。

TRUE TRUE TRUE TRUE FALSE

応用、ある特定の列名をデータフレームから除く

dat <- dat[, !(colnames(dat) %in% c("X", "Y", "Z"))]

http://rmecab.jp/wiki/index.php?R_%CE%F3%CC%BE%A4%C7%CE%F3%A4%F2%BA%EF%BD%FC%A4%B9%A4%EB

2018-10-26

Windows8 のcygwinにpandas をいれる

まず、setup.exeで python3-pip や python2-pipをインストールしておく

pip で最初は失敗したが、rebaseallでできるようになった。

 pip3 install --upgrade pip
python3 -m pip install pandas