diff --git a/R/idiomclassification_mk_pf.R b/R/idiomclassification_mk_pf.R
new file mode 100644
index 0000000..e8c0f87
--- /dev/null
+++ b/R/idiomclassification_mk_pf.R
@@ -0,0 +1,126 @@
+library(caret)
+library(tidyverse)
+library(DMwR)
+library(randomForest)
+library(FSelector)
+# library(randomForestExplainer)
+# may need to: options(expressions = 5e5) to avoid stackoverflow for installing package
+
+setwd(dirname(rstudioapi::getSourceEditorContext()$path))
+stopwords <- readLines(con = "../data/stopwords.txt",encoding="UTF-8")
+oringramme <- read.csv("../data/gold03_anno_ml_synfeat_nstopw.csv", header = TRUE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE)
+syfeaturenames <- read.csv("../data/syfeatures.tsv", header = TRUE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE)
+# syfeaturenames$navalue<-sapply(syfeaturenames$navalue,as.numeric)
+
+deleteStopwords =  function(wl, stopwords = NULL) {
+  wl[!(wl %in% stopwords)]
+}
+
+ngramme <- oringramme %>%
+  add_column(NSTOPW = sapply(oringramme$tokens,function(x) length(deleteStopwords(tolower(unlist(strsplit(x," "))),stopwords)))) %>%
+  # select(-matches("CO_TOKEN.*"), -tokens) %>%
+  select(-matches("CO_TOKEN.*")) %>% # keep tokens for interpretability
+  mutate(across(matches(".rank.*"), ~ replace_na(.x, 1000))) %>%
+  mutate(across(c("dice", "lfmd", "llr", "ld", "pmi"), ~ replace_na(.x, min(.x) - 1))) %>%
+  rename_at(syfeatures$innames, ~ syfeatures[syfeatures$innames==.x,]$synames ) %>%
+  mutate(across(everything(), ~ replace_na(.x, 0))) %>%
+  mutate(CO_IDIOM = as.factor(if_else(CO_IDIOM !=1, 0, 1))) # just two classes: 0 no idiom, 1 idiom
+
+covars <- c("CO_LL", "CO_Z", "CO_G", "CO_T", "CO_LOGDICE", "CO_PMI", "CO_MI3", "CO_DEREKO", "CO_SGT", "CO_WIN5_VEC","CO_WIN5_VEC_AUTOSEM")
+syvars <- c(syfeaturenames$synames,"NSTOPW")
+vars <- c(covars,syvars)
+
+fmla <- as.formula(paste("CO_IDIOM ~ ", paste(vars, collapse= "+")))
+fmlaco <- as.formula(paste("CO_IDIOM ~ ", paste(covars, collapse= "+")))
+fmlasy <- as.formula(paste("CO_IDIOM ~ ", paste(syvars, collapse= "+")))
+
+# Simple train/test split
+
+trainRows <- sample(nrow(ngramme), nrow(ngramme)*0.8, replace = FALSE)
+train <- ngramme[trainRows,]
+test <- ngramme[setdiff(1:nrow(ngramme),trainRows),]
+
+rf_classifier = randomForest(fmla, train, ntree=100, mtry=10, importance=TRUE)
+
+# only SY features
+# rf_classifier = randomForest(fmlasy, train, ntree=100, mtry=10, importance=TRUE)
+
+prediction_for_table <- predict(rf_classifier, test %>% select(-CO_IDIOM))
+
+# different cutoff for prediction
+# prediction_for_table <- predict(rf_classifier, test %>% select(-CO_IDIOM), cutoff = c(0.8, 0.2))
+
+confusion <- table(observed=test$CO_IDIOM,predicted=prediction_for_table)
+conf <- confusionMatrix(confusion)
+print(conf)
+varImpPlot(rf_classifier)
+
+# optional resampling with smote
+
+smoted.data <- SMOTE(fmla, subset(train, select = c("CO_IDIOM", vars)), perc.over = 1200, perc.under = 100)
+rf_classifier = randomForest(fmla, smoted.data, ntree=100, mtry=4, importance=TRUE)
+prediction_for_table <- predict(rf_classifier,test %>% select(-CO_IDIOM))
+confusion <- table(observed=test$CO_IDIOM,predicted=prediction_for_table)
+confusionMatrix(confusion)
+
+# Using estimates by random forest on entire dataset
+
+library(randomForest)
+rf_classifier_full = randomForest(fmla, data=ngramme, ntree=100, mtry=2, importance=TRUE, cutoff=c(0.8,0.2))
+rf_classifier_full
+varImpPlot(rf_classifier_full)
+
+# Feature ranking
+
+# rf features as table
+
+# correlated features seem to split their rankings
+
+rfranks<-importance(rf_classifier_full)[,3:4]
+
+# ttest
+
+idioms<-ngramme %>% filter(CO_IDIOM==1)
+nonidioms<-ngramme %>% filter(CO_IDIOM!=1)
+
+ttestPvalues<-sapply(vars,
+                     function(sel) t.test(idioms[sel],nonidioms[sel])$p.value)
+
+# information gain
+# multiply by 1000 to avoid undersized bins
+# features are ranked individually not matter their correlation
+igain<-information.gain(fmla, data=ngramme%>%mutate_at(vars, ~ . * 1000),unit="log2")
+
+featureRanks<-cbind(rfranks,igain,ttestPvalues)
+
+#randomForestExplainer::explain_forest(rf_classifier )
+
+# averate estimates and feature ranks over 10 runs
+errrate<-0
+conf<-matrix(0,2,3)
+featureRanks<-matrix(0,4,length(vars))
+for (i in 1:10) {
+  # rfc =randomForest(fmla, data=ngramme, ntree=100, importance=TRUE)
+  rfc =randomForest(fmla, data=ngramme, ntree=100, importance=TRUE, cutoff=c(0.8,0.2))
+  errrate<-errrate+rfc$err.rate[100,1]
+  conf<-conf+rfc$confusion
+  featureRanks<-featureRanks+
+    cbind(importance(rfc)[,3:4],
+          sapply(vars,
+                 function(sel) t.test(idioms[sel],nonidioms[sel])$p.value),
+          information.gain(fmla, data=ngramme%>%mutate_at(vars, ~ . * 1000),unit="log2"))
+  print(errrate/i)
+  conf1<-round(
+    rbind(
+      cbind(conf[,1:2]/i,(1-conf[,3]/i)*100),
+      c(100*diag(conf[,1:2])/colSums(conf[,1:2]),NA),
+      c(rowSums(conf[,1:2]/i),NA)),digits=2)
+  colnames(conf1)<-c("0","1","rec")
+  rownames(conf1)<-c("0","1","prec","sum")
+  print(conf1)
+}
+featureRanks<-featureRanks/10
+
+
+
+