独立成分分析

```{r}

#p.71 

library(mvtnorm)

data=data.frame(anscombe)

X=as.matrix(data[,colnames(data) %in% c("x1")])

Y=data$y1

l=nrow(X);n=ncol(X)

rates=sample(1:10,l,replace=T)


log_lik=function(x){
  
return(log(dexp(x,rate=rates)))  
  
}


W=array(0.01,dim=c(l,l))

ite=10000

eta=0.000001;h=0.01

for(s in 1:ite){
  
pre_Y=W%*%X
  
dw=(diag(1,l)+t(t(c(log_lik(pre_Y+h)-log_lik(pre_Y))/h))%*%t(pre_Y))%*%W  
  
W=W+eta*dw  

print(mean(-log_lik(pre_Y)))

#print(mean(-log_lik(pre_Y))+log(dmvnorm(c(pre_Y))))

  
}

```

```{r}

#独立成分分析　村田昇

#不動点法　p.134

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x1","x2")]))/100

pis=data$y

w=rep(0.01,ncol(X))

pthi=function(W){
  
prob=1/(1+exp(-X%*%W))  
  
return(sum(abs(-pis*log(pis)+prob*log(prob))))  
  
}

ite=20

h=0.01

cost=c()

for(l in 1:ite){
  
vec=w;dw=c()
 
for(j in 1:length(w)){
   
vec_sub=vec;vec_sub[j]=vec_sub[j]+h

dw=c(dw,(pthi(vec_sub)-pthi(vec))/h)
   
}

 norms=sqrt(sum(dw^2))
 
 w=-dw/norms
 
 cost=c(cost,pthi(w))
 
print(pthi(w))
  
}


```

```{r}

#独立成分分析　村田昇

#不動点法　p.134 part2

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x1","x2")]))/100

pis=data$y

w=rep(0.01,ncol(X))

pthi=function(W){

prob=1/(1+exp(-X%*%W))  

return(sum(abs(-pis*log(pis)+prob*log(prob))))  

}

ite=2000

h=0.01

eta=0.001

cost=c()

for(l in 1:ite){

vec=w;dw=c()

for(j in 1:length(w)){

vec_sub=vec;vec_sub[j]=vec_sub[j]+h

dw=c(dw,(pthi(vec_sub)-pthi(vec))/h)

}

norms=sqrt(sum(dw^2))

w=w-eta*dw;w=w/sqrt(sum(w^2))

cost=c(cost,pthi(w))

print(pthi(w))

}





```

```{r}

#行列関数への拡張

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x1","x2","x3")]))/100

pis=data$y



pis=array(sample(pis,nrow(X)*ncol(X),replace=T),dim=c(nrow(X),ncol(X)))

w=array(0.01,dim=c(ncol(X),ncol(X)))

pthi=function(W){
  
prob=1/(1+exp(-X%*%W))  
  
return(sum(abs(-pis*log(pis)+prob*log(prob))))  
  
}

ite=30

h=0.01

cost=c()

for(l in 1:ite){
  
vec=w;dw=array(0,dim=dim(w))
 

for(i in 1:dim(w)[1]){
for(j in 1:dim(w)[2]){
   
vec_sub=vec;vec_sub[i,j]=vec_sub[i,j]+h

dw[i,j]=(pthi(vec_sub)-pthi(vec))/h
   
}}
 
 w=-dw;
 
 w.qr=qr(w)
 
 Q=qr.Q(w.qr);R=qr.R(w.qr)
 
 r_vec=ifelse(diag(R)>0,1,-1);S=diag(r_vec)
 
 w=Q%*%S
 
cost=c(cost,pthi(w))
 
print(pthi(w))
  
}





```

```{r}

#行列関数への拡張 part2

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x1","x2","x3")]))/100

pis=data$y



pis=array(sample(pis,nrow(X)*ncol(X),replace=T),dim=c(nrow(X),ncol(X)))

w=array(0.01,dim=c(ncol(X),ncol(X)))

pthi=function(W){

prob=1/(1+exp(-X%*%W))  

return(sum(abs(-pis*log(pis)+prob*log(prob))))  

}

ite=30000

h=0.01

eta=0.001

cost=c()

for(l in 1:ite){

vec=w;dw=array(0,dim=dim(w))


for(i in 1:dim(w)[1]){
for(j in 1:dim(w)[2]){

vec_sub=vec;vec_sub[i,j]=vec_sub[i,j]+h

dw[i,j]=(pthi(vec_sub)-pthi(vec))/h

}}

 w=w-eta*dw;

 w.qr=qr(w)

 Q=qr.Q(w.qr);R=qr.R(w.qr)

 r_vec=ifelse(diag(R)>0,1,-1);S=diag(r_vec)

 w=Q%*%S

cost=c(cost,pthi(w))

print(pthi(w))

}

```