バンディット問題の理論とアルゴリズム (機械学習プロフェッショナルシリーズ)



```{r}
#multi arms bandit

#理論限界 p.26

t=100;n=20

X=array(0,dim=c(t,n))

for(j in 1:n){
  
X[,j]=rpois(t,sample(10:100,1))  
  
}

ave=apply(X,2,mean)

divergence=function(k,l){
  
return(log(l/k)+(k/l)-1)
  
}

E_regret=c()

for(j in 1:t){
  
ave_vec=ave[ave!=min(ave)]  
  
E_regret=c(E_regret,sum((ave_vec-min(ave))*log(j)/(log(min(ave)/ave_vec)+ave_vec/min(ave)-1)))
  
}


```

```{r}


#UCB方策 p.31

t=100;n=20

X=array(0,dim=c(t,n));N=array(0,dim=c(t,n))

N_vec=array(0,dim=c(1,n))

for(j in 1:n){
  
X[,j]=rpois(t,sample(10:100,1))

}

for(j in 1:t){
  
N_vec=N_vec+rpois(n,sample(10:20,1))   

N[j,]=N_vec
    
}

K=sample(1:t,1)

no_vec=c()

for(j in (K+1):t){
  
mu_hat=apply(X[1:(j),],2,sum)/N[(j),]+sqrt(log(j)/(2*N[j,]))

mu_hat_data=data.frame(no=1:n,mu_hat=mu_hat)

no_vec=c(no_vec,min(mu_hat_data$no[mu_hat_data$mu_hat==max(mu_hat)]))
  
}

ave=apply(X,2,mean)

E_regret=log(c(1:t))/sum(min(ave)-ave[ave!=min(ave)])


```

```{r}

t=100;n=20

X=array(0,dim=c(t,n));N=array(0,dim=c(t,n))

N_vec=array(0,dim=c(1,n))

for(j in 1:n){
  
vec=sample(10:100,2,replace=T)  
  
X[,j]=sample(0:1,n,prob=c(vec/sum(vec)),replace=T)

}

N=rep(0,n);M=rep(0,n)

a=0.1;b=0.2

no_vec=c()

for(j in 1:20){
  
mu=c()  
  
for(i in 1:n){
  
mu=c(mu,rbeta(1,a+M[i],b+N[i]-M[i]))  
  
}

mu_data=data.frame(no=1:n,mu=mu)

no=min(mu_data$no[mu_data$mu==max(mu_data$mu)])

no_vec=c(no_vec,no)

N[no]=N[no]+1;M[no]=M[no]+X[j,no]

}


```

```{r}

#Hedge algorithm p.74

t=100;n=20

X=array(0,dim=c(t,n))

for(j in 1:n){
  
X[,j]=rpois(t,sample(10:100,1))  
  
}

w=rep(1/n,n);eta=0.001

no_vec=c()

p_mat=array(0,dim=c(t,n));w_mat=array(0,dim=c(t,n))

for(j in 1:t){
  
p=w/sum(w)  

p_data=data.frame(no=1:n,p=p)

no_vec=c(no_vec,min(p_data$no[p_data$p==max(p_data$p)]))

X_vec=X[j,]

w=w*exp(eta*c(X_vec))

w_mat[j,]=w;p_mat[j,]=p

}



```

```{r}

#Exp3 p.74

t=1000;n=20

X=array(0,dim=c(t,n))

for(j in 1:n){
  
val=rpois(2,sample(5:10,1))  
  
X[,j]=rpois(t,sample(10:100,1))*sample(c(-1,1),t,replace=T,prob=c(val/sum(val)))  
  
}

w=rep(1/n,n);eta=0.001;r=0.90

no_vec=c()

p_mat=array(0,dim=c(t,n));w_mat=array(0,dim=c(t,n))

for(j in 1:t){
  
p=(1-r)*w/sum(w)+r/n  

p_data=data.frame(no=1:n,p=p)

no=min(p_data$no[p_data$p==max(p_data$p)])

no_vec=c(no_vec,no)

X_vec=rep(0,n);X_vec[no]=X_vec[no]+X[j,no]

X_vec=X_vec/p

w=w*exp(eta*c(X_vec))

w_mat[j,]=w;p_mat[j,]=p

}



```

```{r}

#Exp3.p p.79

t=1000;n=20

X=array(0,dim=c(t,n))

for(j in 1:n){
  
val=rpois(2,sample(5:10,1))  
  
X[,j]=rpois(t,sample(10:100,1))*sample(c(-1,1),t,replace=T,prob=c(val/sum(val)))  
  
}

w=rep(1/n,n);eta=0.001;r=0.90;beta=0.1

no_vec=c()

p_mat=array(0,dim=c(t,n));w_mat=array(0,dim=c(t,n))

for(j in 1:t){
  
p=(1-r)*w/sum(w)+r/n  

p_data=data.frame(no=1:n,p=p)

no=min(p_data$no[p_data$p==max(p_data$p)])

no_vec=c(no_vec,no)

X_vec=rep(0,n);X_vec[no]=X_vec[no]+X[j,no]

X_vec=X_vec/p;X_vec[no]=X_vec[no]+beta/p[no]

w=w*exp(eta*c(X_vec))

w_mat[j,]=w;p_mat[j,]=p

}



```

```{r}

#PolyInf p.93

t=100;n=5

X=array(0,dim=c(t,n))

for(j in 1:n){
  
val=rpois(2,sample(5:10,1))  
  
X[,j]=rpois(t,sample(10:100,1))*sample(c(-1,1),t,replace=T,prob=c(val/sum(val)))  
  
}

G=rep(0,n);eta=0.1;r=0.1;

no_vec=c()

p_mat=array(0,dim=c(t,n));G_mat=array(0,dim=c(t,n))

C_fun=function(x){
  
return(log(sum(exp(eta*x))/(1-r))/eta)  
  
}

for(j in 1:t){
  
p=exp(eta*(G-C_fun(G)))+r/n

p_data=data.frame(no=1:n,p=p)

no=min(p_data$no[p_data$p==max(p_data$p)])

no_vec=c(no_vec,no)

X_vec=rep(0,n);X_vec[no]=X_vec[no]+X[j,no]

X_vec=X_vec/p;

G=G+X_vec

G_mat[j,]=G;p_mat[j,]=p

}

regret=min(sqrt(n*t),t)/20


```

```{r}

#逐次削除方策

t=30;n=10

X=array(0,dim=c(t,n))

for(j in 1:n){
  
val=rpois(2,sample(5:10,1))  
  
X[,j]=rpois(t,sample(1:2,1))*sample(c(-1,1),t,replace=T,prob=c(val/sum(val)))  
  
}

colnames(X)=as.character(c(1:n));XX=X

ave=sort(apply(X,2,mean),decreasing = T)

delta=ave[1]-ave

epsilon=1

H_e=1/(2*(ave[1]-ave[2]+epsilon))+sum(1/(2*sum((delta+epsilon*length(delta))^2)))


delta=0.0001;stop=0

samples=c(floor(log(1/delta)/(2*(ave[1]-ave[2]+epsilon)^2)),floor(128*log(4*n/delta)/(ave[1]-ave[-1])))

no_vec=c()

for(j in 1:t){
  
if(stop==0 & ncol(X)>1){  
  
beta=log((4*n*j^2)/delta)

if(j==1){
  
mu_hat=c(X[1,])  
  
}else{

mu_hat=apply(X[1:j,],2,mean)

}

und=mu_hat-sqrt(beta/(2*j))

upp=mu_hat+sqrt(beta/(2*j))

mu_data=data.frame(no=1:n,mu=mu_hat)

no=min(mu_data$no[mu_data$mu==max(mu_data$mu)])

no_vec=c(no_vec,no)

if(und[no]+epsilon>max(upp[-no])){
  
stop=1  

print(paste0("stop_times:",j))
  
}else{
  
 nos=data.frame(no=1:n,sign=ifelse(und[no]>upp,1,0))
 
 #nos=nos[-no,]
 
 nos=nos %>% filter(sign==0);numbers=sort(nos$no)

 ave=ave[numbers];
 
 X=X[,numbers];n=ncol(X)
 
 print(n)
   
}
  
}

}


```

```{r}

#p.118

library(mvtnorm)

t=1000;n=20;m=5

a_vec=array(0,dim=c(1,m+1))

for(i in 1:t){

a_mat=array(0,dim=c(n,m))

for(j in 1:m){
  
val=rpois(2,sample(6:10,1))  
  
a_mat[,j]=sample(0:1,n,replace=T,prob=c(val/sum(val))) 
  
}

a_mat=cbind(rep(i,n),a_mat)

a_vec=rbind(a_vec,a_mat)

}

a_vec=tail(a_vec,nrow(a_vec)-1)

a_mat=data.frame(a_vec)



sigma0=10;sigma=5;b=rep(0,n)

inv_A=diag(sigma0/sigma,n)

sita=rmvnorm(1,inv_A%*%b,sigma*inv_A)

alpha=0.1

for(j in 1:t){
  
sita_pre=sita  
  
sita=rmvnorm(1,inv_A%*%b,sigma*inv_A)  

a_mat_sub=a_mat %>% filter(X1==j);a_mat_sub=as.matrix(a_mat_sub[,2:ncol(a_mat_sub)]) 

mu=c()

for(i in 1:m){
  
vec=a_mat_sub[,i] 

mu=c(mu,as.numeric((vec)%*%t(sita)+alpha*sqrt(log(j))*sigma*sqrt(t(vec)%*%inv_A%*%vec)))
  
}

m_data=data.frame(no=1:m,values=mu)

no=min(m_data$no[m_data$values==max(m_data$values)])

inv_A=inv_A-(inv_A%*%a_mat_sub[,no]%*%t(a_mat_sub[,no])%*%inv_A)/as.numeric(1+t(a_mat_sub[,no])%*%inv_A%*%a_mat_sub[,no])
  
b=b+a_mat_sub[,no]*max(max(m_data$values))

print(sum(abs(sita-sita_pre)))

}

```

```{r}

#p.121

library(mvtnorm)

t=1000;n=20;m=5

a_vec=array(0,dim=c(1,m+1))

for(i in 1:t){

a_mat=array(0,dim=c(n,m))

for(j in 1:m){
  
val=rpois(2,sample(6:10,1))  
  
a_mat[,j]=sample(0:1,n,replace=T,prob=c(val/sum(val))) 
  
}

a_mat=cbind(rep(i,n),a_mat)

a_vec=rbind(a_vec,a_mat)

}

a_vec=tail(a_vec,nrow(a_vec)-1)

a_mat=data.frame(a_vec)



sigma0=10;sigma=5;b=rep(0,n)

inv_A=diag(sigma0/sigma,n)

sita=rmvnorm(1,inv_A%*%b,sigma*inv_A)  

no_vec=c()

for(j in 1:t){
  
sita_pre=sita  
  
sita=rmvnorm(1,inv_A%*%b,sigma*inv_A)  

a_mat_sub=a_mat %>% filter(X1==j);a_mat_sub=as.matrix(a_mat_sub[,2:ncol(a_mat_sub)]) 

m_data=data.frame(no=1:m,values=apply(a_mat_sub*c(sita),2,sum))

no=min(m_data$no[m_data$values==max(m_data$values)])

no_vec=c(no_vec,no)

inv_A=inv_A-(inv_A%*%a_mat_sub[,no]%*%t(a_mat_sub[,no])%*%inv_A)/as.numeric(1+t(a_mat_sub[,no])%*%inv_A%*%a_mat_sub[,no])
  
b=b+a_mat_sub[,no]*max(max(m_data$values))

print(sum(abs(sita-sita_pre)))

}

```

```{r}


library(mvtnorm);library(dplyr)

t=1000;n=30;m=3

a_vec=array(0,dim=c(1,m+1))

vals=array(0,dim=c(m,2))

for(j in 1:m){
  
val=rpois(2,sample(6:10,1))

vals[j,]=val
  
}

for(i in 1:t){

a_mat=array(0,dim=c(n,m))

for(j in 1:m){
  
val=vals[j,]  
  
a_mat[,j]=sample(0:1,n,replace=T,prob=c(val/sum(val))) 
  
}

a_mat=cbind(rep(i,n),a_mat)

a_vec=rbind(a_vec,a_mat)

}

a_vec=tail(a_vec,nrow(a_vec)-1)

a_mat=data.frame(a_vec)

X_s=c()

for(j in 1:t){
  
val=rpois(2,sample(5:10,1))  
  
X_s=c(X_s,sample(c(0:1),1,prob=c(val/sum(val)),replace=T))  
  
}



sigma0=0.1;sigma=5;b=rep(0,n)

sita=t(t(rep(0,n)));no=1;no_vec=c(no)  

prob_train=array(0,dim=c(t,m))

for(j in 1:t){
  
sita_pre=sita  
  
H=diag(1/sigma0,n);G=sita/sigma0

for(i in 1:j){
  
no=no_vec[i]  
  
a_mat_sub=a_mat %>% filter(X1==i);a_mat_sub=as.matrix(a_mat_sub[,2:ncol(a_mat_sub)]) 
 
H=H+exp(as.numeric(t(sita)%*%a_mat_sub[,no]))*(a_mat_sub[,no]%*%t(a_mat_sub[,no]))/((1+exp(as.numeric(t(sita)%*%a_mat_sub[,no])))^2) 

G=G+exp(as.numeric(t(sita)%*%a_mat_sub[,no]))*a_mat_sub[,no]/((1+exp(as.numeric(t(sita)%*%a_mat_sub[,no]))))

if(X_s[i]>0){
  
G=G-a_mat_sub[,no]  
  
}

}

sita=sita-solve(H)%*%G

sita=t(rmvnorm(1,sita,solve(H)))

no_data=data.frame(no=1:m,value=t(a_mat_sub)%*%c(sita))

no=min(no_data$no[no_data$value==max(no_data$value)])

no_vec=c(no_vec,no)

val=exp(t(a_mat_sub)%*%c(sita))/(1+exp(t(a_mat_sub)%*%c(sita)))

prob_train[j,]=c(val)

#print(sum(abs(c(sita_pre)-c(sita))))

#print(paste0("p:",val))
#print(no)

}

no_vec=no_vec[-length(no_vec)]

#result

no_test=c()

prob_test=array(0,dim=c(t,m))

for(i in 1:t){
  
a_mat_sub=a_mat %>% filter(X1==i);a_mat_sub=as.matrix(a_mat_sub[,2:ncol(a_mat_sub)]) 

val=exp(t(a_mat_sub)%*%c(sita))/(1+exp(t(a_mat_sub)%*%c(sita)))

prob_test[i,]=c(val)

no_data=data.frame(no=1:m,value=val)

no=min(no_data$no[no_data$value==max(no_data$value)])

no_test=c(no_test,no) 

print(paste0("p:",val))
  
}

sum(no_vec==no_test)

```