packageを使わないでやってみた（情報理論の応用）

library(dplyr)

#data1:応募をよくするひと,data2:応募をあまりしない人

data1=data2=data.frame(num=1:100,monday=0,tuesday=0,wednesday=0,thursday=0,friday=0,saturday=0,sunday=0)

#だいたい2か月分
N=9;

sample_data=data.frame(num=1:N,monday=0,tuesday=0,wednesday=0,thursday=0,friday=0,saturday=0,sunday=0)

for(j in 2:ncol(data1)){

value=sample(1:9,1);p=value/10;q=1-p

data1[,(j)]=sample(0:1,nrow(data1),replace=TRUE,prob=c(p,q)) 

value=sample(1:9,1);p=value/10;q=1-p

data2[,(j)]=sample(0:1,nrow(data2),replace=TRUE,prob=c(p,q)) 

value=sample(1:9,1);p=value/10;q=1-p

sample_data[,(j)]=sample(0:1,N,replace=TRUE,prob=c(p,q))  


}

p2=apply(data1[,2:ncol(data1)],2,sum)/nrow(data1);q2=1-p2

p1=apply(data2[,2:ncol(data2)],2,sum)/nrow(data2);q1=1-p1


#epsilon=q;1-epsilon=p;x1=0;x2=1;S=0;F=1;

P_y_x1=data.frame(monday=0,tuesday=0,wednesday=0,thursday=0,friday=0,saturday=0,sunday=0)

P_y_x2=data.frame(monday=0,tuesday=0,wednesday=0,thursday=0,friday=0,saturday=0,sunday=0)

P_y_x1[1,]=((1-q1)^apply(sample_data[,2:ncol(sample_data)],2,sum))*(q1^(N-apply(sample_data[,2:ncol(sample_data)],2,sum)))  

P_y_x2[1,]=(p2^apply(sample_data[,2:ncol(sample_data)],2,sum))*((1-p2)^(N-apply(sample_data[,2:ncol(sample_data)],2,sum)))

P_F_y=P_y_x1/(P_y_x1+P_y_x2);P_S_y=P_y_x2/(P_y_x1+P_y_x2)

P_y_x1>P_y_x2

func=function(k){

z<-((1-q1)^k)*q1^(N-k)-(p2^k)*(1-p2)^(N-k)

return(z)

}

#p_y_x1=p_y_x2となる頻度

k_vec=c()

func_data=data.frame(terms=colnames(sample_data[,2:ncol(sample_data)]),upper=0,under=0)

for(l in 1:length(p1)){

times=100
newton_raphson=data.frame(times=1:times,value=0,f=0,df=0)
value=0
for(j in 1:nrow(newton_raphson)){
f<-func(value)[l];h=0.0001
df<-(func(value+h)[l]-func(value)[l])/h

  newton_raphson$f[j]=f
  newton_raphson$df[j]=df
  value=-f/df+value

  newton_raphson$value[j]=value
}

solution=min(newton_raphson$value[!(is.na(newton_raphson$value)>0)])

k_vec=c(k_vec,solution)

h=0.1

plot_data=data.frame(k=seq(0,N,0.1),p_y_x1=0,p_y_x2=0)

for(j in 1:nrow(plot_data)){

plot_data$p_y_x1[j]=((1-q1[l])^plot_data$k[j])*q1[l]^(N-plot_data$k[j]) 

plot_data$p_y_x2[j]=(p2[l]^plot_data$k[j])*(1-p2[l])^(N-plot_data$k[j])

}

func_data$upper[l]=func(solution+h)[l]

func_data$under[l]=func(solution-h)[l]

plot(plot_data$k,plot_data$p_y_x1,col=2,type="o",xlim=c(0,N),ylim=c(0,max(plot_data[,2:ncol(plot_data)])),xlab="k",ylab="prob",main=paste0(colnames(sample_data[,2:ncol(sample_data)][l])))

par(new=T)

plot(plot_data$k,plot_data$p_y_x2,col=3,type="o",xlim=c(0,N),ylim=c(0,max(plot_data[,2:ncol(plot_data)])),xlab="k",ylab="prob",main=paste0(colnames(sample_data[,2:ncol(sample_data)][l])))

}

k_vec

P_e_F=P_e_S=data.frame(terms=colnames(sample_data[,2:ncol(sample_data)]),prob=0)

pbinom=function(p,q,r){

z<-0  

for(n in p:q){

z<-z+(factorial(N)/(factorial(N-n)*factorial(n)))*(r^n)*((1-r)^(N-n))

} 

return(z)  

}

for(l in 1:nrow(P_e_F)){

if(k_vec[l]>=0 & k_vec[l]<=N){

 if(func_data$upper[l]>func_data$under[l]){

  P_e_S$prob[l]=pbinom(0,floor(k_vec[l]),1-q1[l])

  P_e_F$prob[l]=pbinom(floor(k_vec[l])+1,N,p2[l])

 }else{

  P_e_S$prob[l]=pbinom(floor(k_vec[l])+1,N,1-q1[l])

  P_e_F$prob[l]=pbinom(0,floor(k_vec[l]),p2[l])

 } 

}else{

if(k_vec[l]<0){

if(func_data$upper[l]>func_data$under[l]){

P_e_F$prob[l]=pbinom(0,N,1-q1[l])  

}else{

P_e_S$prob[l]=pbinom(0,N,p2[l])  

}

}else{

if(func_data$upper[l]>func_data$under[l]){

P_e_S$prob[l]=pbinom(0,N,p2[l])  

}else{

P_e_F$prob[l]=pbinom(0,N,1-q1[l])  

}  


}  



}


}

colnames(P_e_S)=c("terms","prob_S");colnames(P_e_F)=c("terms","prob_F");

P_e_F_S=left_join(P_e_S,P_e_F)

P_e_F_S=P_e_F_S %>% dplyr::mutate(P_e=(1/2)*(prob_S+prob_F))

#Rate distortion

p=1/3;

#図1.4より、(0,0):0,(0,1):a,(1,0):a,(1,1):0

a=5;

D=max(p*a,(1-p)*a)

q0=function(s){

z<-(p-(1-p)*exp(a*s))/(1-exp(a*s))

return(z)

}

q1=function(s){

z<-((1-p)-p*exp(a*s))/(1-exp(a*s))  

return(z)  

}

p_row=function(x,y){

if(x!=y){

z<-a  

}else{

z<-0  

}  

return(z)  

}

lambda0=function(s){

z<-(q0(s)*exp(s*p_row(0,0))+q1(s)*exp(s*p_row(0,1)))^(-1)  


return(z)  

}

lambda1=function(s){

z<-(q0(s)*exp(s*p_row(1,0))+q1(s)*exp(s*p_row(1,1)))^(-1)

return(z)

}


D=function(s){

  z<-lambda0(s)*p*q1(s)*exp(s*p_row(0,1))*p_row(0,1)+lambda1(s)*(1-p)*q0(s)*exp(s*p_row(1,0))*p_row(1,0)

return(z)

}

Q=function(s){

z<-array(0,dim=c(2,2))

z[1,1]=q0(s)*exp(s*p_row(0,0));
z[1,2]=q0(s)*exp(s*p_row(0,1));
z[2,1]=q1(s)*exp(s*p_row(1,0));
z[2,2]=q1(s)*exp(s*p_row(1,1));

z[1,]=z[1,]/sum(apply(z,1,sum)[1])
z[2,]=z[2,]/sum(apply(z,1,sum)[2])

return(z)


}

R=function(s){

z<-s*D(s)+p*log(lambda0(s))+(1-p)*log(lambda1(s))

return(z)  

}

data=data.frame(s=seq(0,7,0.1),R=0,D=0,lambda1=0,lambda0=0,q0=0,q1=0)

for(j in 1:nrow(data)){

data$R[j]=R(data$s[j])  

data$lambda0[j]=lambda0(data$s[j])

data$lambda1[j]=lambda1(data$s[j])

data$q0[j]=q0(data$s[j])

data$q1[j]=q1(data$s[j])

data$D[j]=D(data$s[j])

}

D_max=p*a

D1=0.1;D2=0.5

convex_data=data.frame(sita=seq(0.1,1,0.1),under_R=0,upper_R=0)

for(j in 1:nrow(convex_data)){

convex_data$under_R[j]=R(convex_data$sita[j]*D1+(1-convex_data$sita[j])*D2)

convex_data$upper_R[j]=convex_data$sita[j]*R(D1)+(1-convex_data$sita[j])*R(D2)


}

plot(data$s,data$R,type="o",col=2,main="R")

plot(data$s,data$D,type="o",col=2,main="D")

plot(data$s,data$lambda0,type="o",col=2,main="lambda1")

plot(data$s,data$lambda1,type="o",col=2,main="lambda2")

plot(data$s,data$q0,type="o",col=2,main="q0")

plot(data$s,data$q1,type="o",col=2,main="q1")

s_max= (1/a)*log(p/(1-p));

print(paste0("s_max=",s_max));print(paste0("伝送率：",log(2)/a))