good turing estimation（メディア情報処理の基礎数理　尾関和彦先生

```{r}

#good turing estimation

#example 1

library(dplyr);library(stringr)

word1=c("Example input and output of SAS statements to do the Seasonal Mann-Whitney-Wilcoxon rank sum test and slope estimator pro- cedureonawater-qualitydatatimeseries")

word2=c("Example input and output of SAS statements to do the Seasonal Kendall test and slope estimator procedure on a water-quality datatimeseries")

train_word=unique(unlist(strsplit(word1," ")))

word=train_word

train_word=data.frame(word=train_word,freq_train=0,freq_test=0)

for(j in 1:nrow(train_word)){
  
train_word$freq_test[j]=sum(str_count(as.character(unlist(strsplit(word2," "))),as.character(word[j])) )
 
train_word$freq_train[j]=sum(str_count(as.character(unlist(strsplit(word1," "))),as.character(word[j])) )
 
}

N=sum(train_word$freq_train)

#P( a=a_i | C(a)=r )

min=min(train_word$freq_train);max=max(train_word$freq_train)

r=1

prob=r/sum(train_word$freq_train[train_word$freq_train==r])

#E_N[Nr]

sum(train_word$freq_train[train_word$freq_train==r]/N)

#E[P(a)|C(a)=r]

((r+1)/(N+1))*(sum(train_word$freq_train[train_word$freq_train==(r+1)]/(N+1))/sum(train_word$freq_train[train_word$freq_train==r]/N))


#r*

N*((r+1)/(N+1))*(sum(train_word$freq_train[train_word$freq_train==(r+1)]/(N+1))/sum(train_word$freq_train[train_word$freq_train==r]/N))


```

```{r}

#example 2

library(dplyr);library(stringr)

word1=c("This report describes procedures in the Statistical Analysis System (SAS) thatcanappropriatelybeusedtodetecttrendsinwater-qualitydata The report describes the use of both the standard procedures provided with SAS and twoadditionalprocedures,SEASKENandSEASRS AworkingknowledgeofSASby thereaderisassumed Thereportwaswrittenprimarilyforuser'softhe
U.S. Geological Survey's WATSTORE data base and Amdahl computer system; however, the source code for the SAS macros and procedures are included in the appendices. The user's guide for the SEASKEN and SEASRS procedures are also included in the appendices")

word2=c("Forpurposesofthistest,themeanstreamflow maybeusediftheinstantaneousstreamflowismissing. PROCQWSAScreatesa
SAS data set named DATA1 containing the station identification number and the eleven parameters requested in the QWRETR procedure. The parameter values are stored in variables named Pnnnnn where nnnnn is the WATSTORE parameter codes given in the retrieval list. In addition, the variables YEAR, MONTH, DAY, DATE, DECTIME, and SNAME were requested in the QWSAS statement.")

train_word=unique(unlist(strsplit(word1," ")))

word=train_word

train_word=data.frame(word=train_word,freq_train=0,freq_test=0)

for(j in 1:nrow(train_word)){
  
train_word$freq_test[j]=sum(str_count(as.character(unlist(strsplit(word2," "))),as.character(word[j])) )
 
train_word$freq_train[j]=sum(str_count(as.character(unlist(strsplit(word1," "))),as.character(word[j])) )
 
}

N=sum(train_word$freq_train)

#P( a=a_i | C(a)=r )

min=min(train_word$freq_train);max=max(train_word$freq_train)

r=2

prob=r/sum(train_word$freq_train[train_word$freq_train==r])

#E_N[Nr]

sum(train_word$freq_train[train_word$freq_train==r]/N)

#E[P(a)|C(a)=r]

((r+1)/(N+1))*(sum(train_word$freq_train[train_word$freq_train==(r+1)]/(N+1))/sum(train_word$freq_train[train_word$freq_train==r]/N))


#r*

N*((r+1)/(N+1))*(sum(train_word$freq_train[train_word$freq_train==(r+1)]/(N+1))/sum(train_word$freq_train[train_word$freq_train==r]/N))

#Nr and r are linear 

linear_data=train_word %>% group_by(freq_train) %>% summarise(n=n()) %>% mutate(log_n=log(n),log_r=log(freq_train))

cor(linear_data$log_n,linear_data$log_r)

a=cov(linear_data$log_n,linear_data$log_r)/cor(linear_data$log_r,linear_data$log_r)

b=mean(linear_data$log_n)-a*mean(linear_data$log_r)

linear_data=linear_data %>% mutate(predict=a*log_r+b)

linear_data=linear_data %>% mutate(S=0,N=0)

for(j in 2:(nrow(linear_data)-1)){
  
linear_data$S[j]=2*linear_data$n[j]/(linear_data$freq_train[j+1]-linear_data$freq_train[j-1])  
  
}

S_linear_data=linear_data %>% filter(S>0)

#linear corresponce(S)

cor(S_linear_data$log_r,log(S_linear_data$S))

#linear corresponce(log_n)

cor(linear_data$log_n,linear_data$predict)

log_n=linear_data$log_n;log_r=linear_data$log_r;predict=linear_data$predict;log_S=log(S_linear_data$S);s_log_r=S_linear_data$log_r

a_s=cov(s_log_r,log_S)/cov(s_log_r,s_log_r)

b_s=mean(log_S)-a_s*mean(s_log_r)

for(j in 2:(nrow(linear_data)-1)){
  
#N:fixed  
  
linear_data$N[j]=(linear_data$S[j]/2)*(linear_data$freq_train[j+1]-linear_data$freq_train[j-1])  
  
}

```

```{r}

#log_r plot

plot(log_r,log_n,type="o",col=2,xlim=c(min(log_r),max(log_r)),ylim=c(min(log_n),max(log_n)),xlab="log_r",ylab="log_n")

par(new=T)

plot(log_r,predict,type="o",col=3,xlim=c(min(log_r),max(log_r)),ylim=c(min(log_n),max(log_n)),xlab="log_r",ylab="log_n")


```

```{r}

#log_S plot

plot(s_log_r,log_S,type="o",col=2,xlim=c(min(s_log_r),max(s_log_r)),ylim=c(min(log_S),max(log_S)),xlab="s_log_r",ylab="log_S")

par(new=T)

plot(s_log_r,a_s*s_log_r+b_s,type="o",col=3,xlim=c(min(s_log_r),max(s_log_r)),ylim=c(min(log_S),max(log_S)),xlab="s_log_r",ylab="log_S")


```