diff --git a/README.md b/README.md
index 71aabd3..c44a4e2 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,3 @@
 # CTRF
-Causal Transfer Random Forests (CTRF)
+# This is an anonymous repository for the paper: "Causal Transfer Random Forest: Combining Logged Data and Randomized Experiments for Robust Prediction". It contains the Supplementary Material "WSDM_Supplementary_Material.jpg" (for anonymous reason, we only keep it as image yet will release the pdf version upon acceptance) and code for reproducing the results. The details are included in the Supplementary Material. We do not make any further changes after submitting the paper.
 
-Simple starting point for CTRF method.
\ No newline at end of file
diff --git a/WSDM_Supplmentary_Material.jpg b/WSDM_Supplmentary_Material.jpg
new file mode 100644
index 0000000..9d022a7
Binary files /dev/null and b/WSDM_Supplmentary_Material.jpg differ
diff --git a/ctrf/auction.py b/ctrf/auction.py
index 13bd409..edbfd29 100644
--- a/ctrf/auction.py
+++ b/ctrf/auction.py
@@ -14,7 +14,14 @@ def run_selection(seed, n_samples, auction_size, n_auctions):
     ind = np.random.randint(0, n_samples, size=auction_size*n_auctions)
     return ind, seed+1
 
-def run_auction(dataset, seed, model, epsilon, auction_size, n_auctions, max_slate):
+def sample_true(x):
+    if any(x):
+        return np.random.choice(np.where(x)[0])
+    else:
+        return np.random.choice(np.where(x==False)[0])
+    
+
+def run_auction(dataset, seed, model, epsilon, auction_size, n_auctions, max_slate, position_effect=0):
     seed += 1
     np.random.seed(seed)
 
@@ -67,19 +74,33 @@ def run_auction(dataset, seed, model, epsilon, auction_size, n_auctions, max_sla
         df['Layout'] = df.groupby('AuctionId')['AuctionId'].transform('count')
 
     # Rank by PClick, then cascade to generate clicks
+    # Effect for Position is Zero
     df['Uniform'] = np.random.uniform(size=len(df))
     df['WouldClick'] = np.where(df['Uniform'] <= df['TruePClick'], 1, 0)
     df['Click'] = 0
-    df.loc[df["WouldClick"].ne(0).groupby(df['AuctionId']).idxmax(),'Click']=1
-    df['Click'] = df['Click'] * df['WouldClick']
+    if position_effect==1:
+        df.loc[df["WouldClick"].ne(0).groupby(df['AuctionId']).idxmax(),'Click']=1
+    else:
+        sample_id = df['WouldClick'].ne(0).groupby(df['AuctionId']).apply(sample_true)
+        group_id = df['WouldClick'].ne(0).groupby(df['AuctionId']).groups
+        idx=[]
+        for k,v in group_id.items():
+            idx.append(group_id[k][sample_id[k]])
+        df.loc[idx,'Click']=1
 
+    df['Click'] = df['Click'] * df['WouldClick']
     df.drop(columns=['Uniform', 'WouldClick', 'RankingPClick', 'TruePClick'], inplace=True)
 
+
+
+
     return df, seed+1
 
 def construct_auction_dataset(dataset):
     X = np.hstack((dataset['auctions'][['PClick', 'Position', 'Layout']], dataset['X'][dataset['auctions']['SampleId']]))
-    y = dataset['y'][dataset['auctions']['SampleId']]
+    #y = dataset['y'][dataset['auctions']['SampleId']]
+    y = dataset['auctions']['Click']
+
 
     return X, y
 
diff --git a/ctrf/metrics.py b/ctrf/metrics.py
index 471e210..fd2cc50 100644
--- a/ctrf/metrics.py
+++ b/ctrf/metrics.py
@@ -1,7 +1,8 @@
 import operator
+import numpy as np
+from sklearn.metrics import roc_auc_score,f1_score,log_loss
 
-from sklearn.metrics import roc_auc_score
-
+###AUC Metric
 def compute_auc(preds, ys):
     preds = sorted(zip(preds, ys), key=operator.itemgetter(0), reverse=True)
     pred_p_te, pred_y_te = zip(*preds)
@@ -11,3 +12,29 @@ def compute_model_auc(model, x_te, y_te):
     x_te = x_te.copy()
     y_te = y_te.copy()
     return compute_auc(model.predict_proba(x_te)[:, 1], y_te)
+
+###F1 Score
+def compute_model_f1(model, x_te, y_te):
+    x_te = x_te.copy()
+    y_te = y_te.copy()
+    return f1_score(model.predict(x_te), y_te)
+
+###Bias
+def compute_model_bias(model, x_te,y_te):
+    x_te = x_te.copy()
+    y_te = y_te.copy()
+    pred_ctr = np.mean(model.predict_proba(x_te)[:, 1])
+    real_ctr = np.mean(y_te)
+    return  abs(pred_ctr-real_ctr)/real_ctr
+
+###RIG
+def compute_model_rig(model, x_te, y_te):
+    x_te = x_te.copy()
+    y_te = y_te.copy()
+    real_ctr = np.mean(y_te)
+    real_entropy = - np.log(real_ctr)*real_ctr- np.log(1-real_ctr)*(1-real_ctr)
+    pred_y = model.predict_proba(x_te)[:, 1]
+    l_score = - log_loss(y_te, pred_y)
+
+    return (real_entropy+l_score)/real_entropy
+
diff --git a/ctrf/models.py b/ctrf/models.py
index 24faa5a..88414c2 100644
--- a/ctrf/models.py
+++ b/ctrf/models.py
@@ -1,7 +1,11 @@
 import pandas
 import numpy as np
 import time
-from sklearn.ensemble import RandomForestClassifier
+from sklearn.ensemble import RandomForestClassifier,GradientBoostingClassifier 
+from sklearn.linear_model import LogisticRegression
+# from sklearn import svm
+
+from copy import deepcopy
 
 def train_rf(X, y, seed, **kwargs):
     start = time.time()
@@ -12,15 +16,81 @@ def train_rf(X, y, seed, **kwargs):
     seed += 1
     return model, seed
 
-def train_ctrf(X1, y1, X2, y2, seed, **kwargs):
+def train_combine_rf(X1, y1, X2, y2, seed, **kwargs):
+    start = time.time()
     seed += 1
     model = RandomForestClassifier(criterion='entropy', random_state=seed, **kwargs)
-    model.fit(X1, y1)
+    X_c = np.concatenate([X1, X2], axis=0)
+    y_c = np.concatenate([y1, y2])
+    model.fit(X_c, y_c)
+    print('Runtime:', time.time()-start)
+    seed += 1
+    return model, seed
 
+def train_ctrf(X1, y1, X2, y2, model, combine, seed,**kwargs):
+    seed += 1
+    # model = RandomForestClassifier(criterion='entropy', random_state=seed, **kwargs)
+    # model.fit(X1, y1)
+    #Copy an Oject
+    start = time.time()
+    model = deepcopy(model)
     for e in model.estimators_:
-        df = pandas.DataFrame(zip(e.apply(X2), 1-y2, y2), columns=['LeafId', 'NoClick', 'Click'])
+        if combine==0:
+            df = pandas.DataFrame(zip(e.apply(X2), 1-y2, y2), columns=['LeafId', 'NoClick', 'Click'])
+        else:
+            X_c = np.concatenate([X1,X2],axis=0)
+            y_c = np.concatenate([y1,y2])
+            df = pandas.DataFrame(zip(e.apply(X_c), 1-y_c, y_c), columns=['LeafId', 'NoClick', 'Click'])
         df = df.groupby(['LeafId']).agg(NoClicks=pandas.NamedAgg(column='NoClick', aggfunc='sum'), Clicks=pandas.NamedAgg(column='Click', aggfunc='sum'))
         e.tree_.value[df.index.array] = np.expand_dims(df[['NoClicks', 'Clicks']].values, axis=1)
+    print('Runtime:', time.time()-start)
+    seed += 1
+    return model, seed
+
 
+def calculate_weight(train_X,testing_X):
+    pool_X=np.vstack([train_X,testing_X])
+    pool_Y=np.hstack([np.zeros(train_X.shape[0]),np.ones(testing_X.shape[0])])
+    model=LogisticRegression(solver='liblinear')
+    model.fit(pool_X,pool_Y)
+    pred=model.predict_proba(train_X)
+    weights=pred[:,1]/pred[:,0]
+    weights=weights/np.mean(weights)
+    return weights
+
+def train_lr_model(X,y,seed,**kwargs):
+    start = time.time()
+    seed += 1
+    model = LogisticRegression(solver='liblinear', random_state=seed)
+    model.fit(X,y)
+    print('Runtime:', time.time() - start)
     seed += 1
-    return model, seed
\ No newline at end of file
+    return model,seed
+
+def train_gbdt_model(X,y,seed,**kwargs):
+    start = time.time()
+    seed += 1
+    model = GradientBoostingClassifier(random_state=seed)
+    model.fit(X,y)
+    print('Runtime:', time.time() - start)
+    seed += 1
+    return model,seed
+
+def train_lr_weight_model(X,y,weights,seed,**kwargs):
+    start = time.time()
+    seed += 1
+    model = LogisticRegression(solver='liblinear', random_state=seed)
+    model.fit(X,y,sample_weight=weights)
+    print('Runtime:', time.time() - start)
+    seed += 1
+    return model,seed
+
+def train_gbdt_weight_model(X,y,weights,seed,**kwargs):
+    start = time.time()
+    seed += 1
+    model = GradientBoostingClassifier(random_state=seed)
+    model.fit(X,y,sample_weight=weights)
+    print('Runtime:', time.time() - start)
+    seed += 1
+    return model,seed
+
diff --git a/ctrf/util.py b/ctrf/util.py
new file mode 100644
index 0000000..dbf197a
--- /dev/null
+++ b/ctrf/util.py
@@ -0,0 +1,69 @@
+import numpy as np
+def simu_confounding_data(n=2000,p=20,scenario=2,r=0.65,binary=1):
+    simu_data={'Y':[],'S':[],'V':[],'X':[],'r':r,'scenario':scenario}
+
+    p_s=int(p*0.4)
+    p_v=int(p*0.6)
+
+    i_grid=np.linspace(1,p_s,p_s)
+    alpha=(-1)**i_grid*(i_grid%3+1)*p/3
+    beta=p/2
+    count=0
+    while count<n:
+        if scenario==1:
+            S=np.random.normal(size=p_s)
+            V = np.random.normal(size=p_v)
+        elif scenario==2:
+            #S Causes V
+            S=np.random.normal(size=p_s)
+            V = np.zeros(p_v)
+            for j in range(p_v):
+                V[j]=np.random.normal(loc=int(S[j%p_s]>0)+int(S[(j+1)%p_s]>0))
+        else:
+            #V Causes S
+            V=np.random.normal(size=p_v)
+            S = np.zeros(p_s)
+            for j in range(p_s):
+                S[j]=np.random.normal(loc=int(V[j%p_v]>0)+int(V[(j+1)%p_v]>0))
+
+        S_obs=np.zeros(p_s)
+        V_obs=np.zeros(p_v)
+        S_obs[np.where(S>0)]=1
+        V_obs[np.where(V>0)]=1
+
+        logit=np.sum(alpha*S_obs)+(np.sum(S_obs[1:]*S_obs[:p_s-1]))*beta
+
+        Y=1/(1+np.exp(-logit))+np.random.normal(scale=0.2)
+        if binary==1:
+            Y_obs=0
+            if Y>0.5:
+                Y_obs=1
+        else:
+            Y_obs=Y
+
+        noisy_mean=np.mean(V_obs)
+        inclusion=np.random.uniform()
+        ###Positive Correlation
+        if (Y_obs>0.5 and noisy_mean>0.5) or (Y_obs<0.5 and noisy_mean<0.5):
+            if inclusion<r:
+                simu_data['Y'].append(Y_obs)
+                simu_data['S'].append(S_obs)
+                simu_data['V'].append(V_obs)
+                simu_data['X'].append(np.hstack([S_obs,V_obs]))
+                count+=1
+        ###Negative Correlation
+        else:
+            if inclusion<(1-r):
+                simu_data['Y'].append(Y_obs)
+                simu_data['S'].append(S_obs)
+                simu_data['V'].append(V_obs)
+                simu_data['X'].append(np.hstack([S_obs,V_obs]))
+                count+=1
+    #Into Array
+    simu_data['Y']=np.asarray(simu_data['Y'])
+    simu_data['X']=np.asarray(simu_data['X'])
+    simu_data['V']=np.asarray(simu_data['V'])
+    simu_data['S']=np.asarray(simu_data['S'])
+    
+    return simu_data
+
diff --git a/run_auction.py b/run_auction.py
index 35a6d22..d9277f9 100644
--- a/run_auction.py
+++ b/run_auction.py
@@ -1,10 +1,12 @@
 import argparse
 import numpy as np
+import os
+import pickle
 from sklearn.datasets import make_classification
-
-from ctrf.metrics import compute_auc, compute_model_auc
-from ctrf.models import train_rf, train_ctrf
-from ctrf.auction import run_selection, run_auction, construct_auction_dataset
+#
+from ctrf.metrics import *
+from ctrf.models import *
+from ctrf.auction import *
 
 if __name__ == '__main__':
     parser = argparse.ArgumentParser()
@@ -18,6 +20,7 @@
     parser.add_argument('--n_informative', default=10, type=int, help='Relevant Features')
     parser.add_argument('--n_clusters_per_class', default=10, type=int, help='Clusters per class')
     parser.add_argument('--class_sep', default=1, type=float, help='Class separation')
+    parser.add_argument('--n_experiments', default=100, type=int, help='Number of Experiments')
 
     parser.add_argument('--oracle_n_estimators', default=100, type=int, help='Number of estimators for Oracle')
     parser.add_argument('--oracle_min_samples_leaf', default=100, type=int, help='Minimum number of samples for Oracle to use for labeling')
@@ -33,53 +36,137 @@
     parser.add_argument('--control_reserve', default=.5, type=float, help='Reserve on Control Flight')
     parser.add_argument('--treatment_reserve', default=.7, type=float, help='Reserve on Treatment Flight')
     parser.add_argument('--max_slate', default=5, type=float, help='Maximum slate size')
+    parser.add_argument('--different_slate', default=0, type=float, help='Difference in slate size')
     args = parser.parse_args()
 
     print('Arguments', args)
-    X, y =  make_classification(  n_samples=args.n_tr_samples + args.n_te_cnt_samples + args.n_te_trt_samples + args.n_te_rnd_samples
-                                , n_features=args.n_features
-                                , n_informative=args.n_informative
-                                , n_redundant=0
-                                , n_clusters_per_class=args.n_clusters_per_class
-                                , class_sep=args.class_sep
-                                , random_state=args.seed)
-    args.seed += 1
-
-    print('Splitting Datasets')
-    samples = np.array([0, args.n_tr_samples, args.n_te_cnt_samples, args.n_te_trt_samples, args.n_te_rnd_samples])
-    datasets = {name : {'name':name, 'start':start, 'end':end, 'samples':samples, 'auction_type':auction_type, 'reserve':reserve} for name, start, end, samples, auction_type, reserve in zip(['oracle', 'cnt', 'trt', 'rnd'], np.cumsum(samples), np.cumsum(samples)[1:], samples[1:], [None, 'greedy', 'greedy', 'random'], [None, args.control_reserve, args.treatment_reserve, 0])}
-
-    for dataset, info in datasets.items():
-        print(dataset, info['start'], info['end'], info['samples'])
-        info['X'], info['y'] = X[info['start']:info['end'], :].copy() , y[info['start']:info['end']].copy()
-
-    # Creates an oracle pclick that ignores position and simply observes c/nc
-    oracle, args.seed = train_rf(datasets['oracle']['X'], datasets['oracle']['y'], seed=args.seed, n_estimators=args.oracle_n_estimators, min_samples_leaf=args.oracle_min_samples_leaf)
-    print('Oracle Created -- Test AUC Control:', compute_model_auc(oracle, datasets['cnt']['X'], datasets['cnt']['y']))
-    print('Oracle Created -- Test AUC Treatment:', compute_model_auc(oracle, datasets['trt']['X'], datasets['trt']['y']))
-
-    # Run selection
-    datasets['rnd']['ind'], args.seed = run_selection(args.seed, datasets['rnd']['samples'], args.auction_size, args.n_rnd_auction)
-    datasets['cnt']['ind'], args.seed = run_selection(args.seed, datasets['cnt']['samples'], args.auction_size, args.n_auctions)
-    datasets['trt']['ind'], args.seed = run_selection(args.seed, datasets['trt']['samples'], args.auction_size, args.n_auctions)
-
-    # Run the auction
-    datasets['rnd']['auctions'], args.seed = run_auction(datasets['rnd'], args.seed, oracle, args.epsilon, args.auction_size, args.n_rnd_auction, args.max_slate)
-    datasets['cnt']['auctions'], args.seed = run_auction(datasets['cnt'], args.seed, oracle, args.epsilon, args.auction_size, args.n_auctions, args.max_slate)
-    datasets['trt']['auctions'], args.seed = run_auction(datasets['trt'], args.seed, oracle, args.epsilon, args.auction_size, args.n_auctions, args.max_slate)
-
-    print('Train RF Models')
-    rnd_rf, args.seed = train_rf(*construct_auction_dataset(datasets['rnd']), seed=args.seed, n_estimators=args.auction_n_estimators, max_leaf_nodes=args.auction_max_leaf_nodes)
-    cnt_rf, args.seed = train_rf(*construct_auction_dataset(datasets['cnt']), seed=args.seed, n_estimators=args.auction_n_estimators, max_leaf_nodes=args.auction_max_leaf_nodes)
-    trt_rf, args.seed = train_rf(*construct_auction_dataset(datasets['trt']), seed=args.seed, n_estimators=args.auction_n_estimators, max_leaf_nodes=args.auction_max_leaf_nodes)
-    cnt_ctrf, args.seed = train_ctrf(*construct_auction_dataset(datasets['rnd']), *construct_auction_dataset(datasets['cnt']), seed=args.seed, n_estimators=args.auction_n_estimators, max_leaf_nodes=args.auction_max_leaf_nodes)
-    trt_ctrf, args.seed = train_ctrf(*construct_auction_dataset(datasets['rnd']), *construct_auction_dataset(datasets['trt']), seed=args.seed, n_estimators=args.auction_n_estimators, max_leaf_nodes=args.auction_max_leaf_nodes)
-
-    print('Eval rnd_rf on cnt', compute_model_auc(rnd_rf, *construct_auction_dataset(datasets['cnt'])))
-    print('Eval rnd_rf on trt', compute_model_auc(rnd_rf, *construct_auction_dataset(datasets['trt'])))
-
-    print('Eval cnt_rf on trt', compute_model_auc(cnt_rf, *construct_auction_dataset(datasets['trt'])))
-    print('Eval trt_rf on cnt', compute_model_auc(trt_rf, *construct_auction_dataset(datasets['cnt'])))
-
-    print('Eval cnt_ctrf on trt', compute_model_auc(cnt_ctrf, *construct_auction_dataset(datasets['trt'])))
-    print('Eval trt_ctrf on cnt', compute_model_auc(trt_ctrf, *construct_auction_dataset(datasets['cnt'])))
+
+
+    def results_update(model, results, data_name, method_name):
+        if method_name == 'oracle':
+            results[data_name]['auc'][method_name].append(
+                compute_model_auc(model, datasets[data_name]['X'], datasets[data_name]['y']))
+            results[data_name]['f1_score'][method_name].append(
+                compute_model_f1(model, datasets[data_name]['X'], datasets[data_name]['y']))
+            results[data_name]['bias'][method_name].append(
+                compute_model_bias(model, datasets[data_name]['X'], datasets[data_name]['y']))
+            results[data_name]['rig'][method_name].append(
+                compute_model_rig(model, datasets[data_name]['X'], datasets[data_name]['y']))
+        else:
+            results[data_name]['auc'][method_name].append(
+                compute_model_auc(model, *construct_auction_dataset(datasets[data_name])))
+            results[data_name]['f1_score'][method_name].append(
+                compute_model_f1(model, *construct_auction_dataset(datasets[data_name])))
+            results[data_name]['bias'][method_name].append(
+                compute_model_bias(model, *construct_auction_dataset(datasets[data_name])))
+            results[data_name]['rig'][method_name].append(
+                compute_model_rig(model, *construct_auction_dataset(datasets[data_name])))
+            results[data_name]['feature_imp'][method_name].append(np.argsort(model.feature_importances_))
+
+
+    result_name = '_'.join(
+        ["results", str(args.n_rnd_auction), str(args.n_auctions), str(int(args.treatment_reserve * 100)),
+         str(int(args.control_reserve * 100)), str(args.max_slate),str(args.different_slate)])
+
+    ##Results to Save
+
+    method_list = ['oracle', 'rnd_rf', 'trt_rf', 'cnt_rf', 'trt_ctrf', 'cnt_ctrf', 'combine_rf']
+    results = {name: {'name': name, 'auction_type': auction_type, 'reserve': reserve,
+                      'auc_sample': [],
+                      'auc': {method_name: [] for method_name in method_list},
+                      'f1_score': {method_name: [] for method_name in method_list},
+                      'rig': {method_name: [] for method_name in method_list},
+                      'bias': {method_name: [] for method_name in method_list},
+                      'feature_imp': {method_name: [] for method_name in method_list}}
+               for name, auction_type, reserve in
+               zip(['oracle', 'cnt', 'trt', 'rnd'], [None, 'greedy', 'greedy', 'random'],
+                   [None, args.control_reserve, args.treatment_reserve, 0])}
+    if os.path.isdir('results'):
+         print('result directory existed')
+    else:
+        os.mkdir("results")
+        print("result directory added")
+
+    path = "results/"
+    os.chdir(path)
+
+    for i in range(args.n_experiments):
+
+        print ('------------------------------')
+        print ('Run %d th experiments:' % (i+1))
+        X, y =  make_classification(  n_samples=args.n_tr_samples + args.n_te_cnt_samples + args.n_te_trt_samples + args.n_te_rnd_samples
+                                    , n_features=args.n_features
+                                    , n_informative=args.n_informative
+                                    , n_redundant=0
+                                    , n_clusters_per_class=args.n_clusters_per_class
+                                    , class_sep=args.class_sep
+                                    , random_state=args.seed)
+        args.seed += 1
+
+        #print('Splitting Datasets')
+        samples = np.array([0, args.n_tr_samples, args.n_te_cnt_samples, args.n_te_trt_samples, args.n_te_rnd_samples])
+        datasets = {name : {'name':name, 'start':start, 'end':end, 'samples':samples, 'auction_type':auction_type, 'reserve':reserve} for name, start, end, samples, auction_type, reserve in zip(['oracle', 'cnt', 'trt', 'rnd'], np.cumsum(samples), np.cumsum(samples)[1:], samples[1:], [None, 'greedy', 'greedy', 'random'], [None, args.control_reserve, args.treatment_reserve, 0])}
+
+        for dataset, info in datasets.items():
+        #    print(dataset, info['start'], info['end'], info['samples'])
+            info['X'], info['y'] = X[info['start']:info['end'], :].copy() , y[info['start']:info['end']].copy()
+
+       # Creates an oracle pclick that ignores position and simply observes c/nc
+        oracle, args.seed = train_rf(datasets['oracle']['X'], datasets['oracle']['y'], seed=args.seed, n_estimators=args.oracle_n_estimators, min_samples_leaf=args.oracle_min_samples_leaf)
+        # print('Oracle Created -- Test AUC Control:', compute_model_auc(oracle, datasets['cnt']['X'], datasets['cnt']['y']))
+        # print('Oracle Created -- Test AUC Treatment:', compute_model_auc(oracle, datasets['trt']['X'], datasets['trt']['y']))
+
+        # Run selection
+        datasets['rnd']['ind'], args.seed = run_selection(args.seed, datasets['rnd']['samples'], args.auction_size, args.n_rnd_auction)
+        datasets['cnt']['ind'], args.seed = run_selection(args.seed, datasets['cnt']['samples'], args.auction_size, args.n_auctions)
+        datasets['trt']['ind'], args.seed = run_selection(args.seed, datasets['trt']['samples'], args.auction_size, args.n_auctions)
+
+        # Run the auction
+        datasets['rnd']['auctions'], args.seed = run_auction(datasets['rnd'], args.seed, oracle, args.epsilon, args.auction_size, args.n_rnd_auction, args.max_slate)
+        datasets['cnt']['auctions'], args.seed = run_auction(datasets['cnt'], args.seed, oracle, args.epsilon, args.auction_size, args.n_auctions, args.max_slate)
+        datasets['trt']['auctions'], args.seed = run_auction(datasets['trt'], args.seed, oracle, args.epsilon, args.auction_size, args.n_auctions, args.max_slate+args.different_slate)
+
+        print('Train RF Models')
+        rnd_rf, args.seed = train_rf(*construct_auction_dataset(datasets['rnd']), seed=args.seed, n_estimators=args.auction_n_estimators, max_leaf_nodes=args.auction_max_leaf_nodes)
+        cnt_rf, args.seed = train_rf(*construct_auction_dataset(datasets['cnt']), seed=args.seed, n_estimators=args.auction_n_estimators, max_leaf_nodes=args.auction_max_leaf_nodes)
+        trt_rf, args.seed = train_rf(*construct_auction_dataset(datasets['trt']), seed=args.seed, n_estimators=args.auction_n_estimators, max_leaf_nodes=args.auction_max_leaf_nodes)
+        combine_rf, args.seed = train_combine_rf(*construct_auction_dataset(datasets['rnd']),*construct_auction_dataset(datasets['cnt']), seed=args.seed, n_estimators=args.auction_n_estimators, max_leaf_nodes=args.auction_max_leaf_nodes)
+
+        print('Train CTRF Models')
+        cnt_ctrf, args.seed = train_ctrf(*construct_auction_dataset(datasets['rnd']), *construct_auction_dataset(datasets['cnt']), rnd_rf, combine=1,seed=args.seed, n_estimators=args.auction_n_estimators, max_leaf_nodes=args.auction_max_leaf_nodes)
+        trt_ctrf, args.seed = train_ctrf(*construct_auction_dataset(datasets['rnd']), *construct_auction_dataset(datasets['trt']), rnd_rf,combine=1, seed=args.seed, n_estimators=args.auction_n_estimators, max_leaf_nodes=args.auction_max_leaf_nodes)
+
+        # Results_Update
+        for model_name in ['oracle', 'rnd_rf', 'cnt_rf', 'trt_rf', 'cnt_ctrf', 'trt_ctrf', 'combine_rf']:
+            if model_name == 'oracle':
+                try:
+                    results_update(eval(model_name), results, 'oracle', model_name)
+                except:
+                    continue
+            for data_name in ['rnd', 'cnt', 'trt']:
+                try:
+                    results_update(eval(model_name), results, data_name, model_name)
+                except:
+                    continue
+                # Sample size of impressions
+                results[data_name]['auc_sample'].append(construct_auction_dataset(datasets[data_name])[0].shape[0])
+        print ('Finish %d th experiments:' % (i+1))
+
+
+        f = open(result_name + ".pkl", "wb")
+        pickle.dump(results, f)
+        f.close()
+
+    # # Loading
+    # with open(result_name + '.pkl', 'rb') as pickle_file:
+    #     content = pickle.load(pickle_file)
+
+    # print('Eval rnd_rf on cnt', compute_model_auc(rnd_rf, *construct_auction_dataset(datasets['cnt'])))
+
+    # print('Eval rnd_rf on trt', compute_model_auc(rnd_rf, *construct_auction_dataset(datasets['trt'])))
+    #
+    # print('Eval cnt_rf on trt', compute_model_auc(cnt_rf, *construct_auction_dataset(datasets['trt'])))
+    # print('Eval trt_rf on cnt', compute_model_auc(trt_rf, *construct_auction_dataset(datasets['cnt'])))
+    #
+    # print('Eval cnt_ctrf on trt', compute_model_auc(cnt_ctrf, *construct_auction_dataset(datasets['trt'])))
+    # print('Eval trt_ctrf on cnt', compute_model_auc(trt_ctrf, *construct_auction_dataset(datasets['cnt'])))
diff --git a/run_simulation.py b/run_simulation.py
new file mode 100644
index 0000000..166f414
--- /dev/null
+++ b/run_simulation.py
@@ -0,0 +1,82 @@
+import argparse
+import numpy as np
+from ctrf.util import *
+from ctrf.metrics import *
+from ctrf.models import *
+from ctrf.auction import *
+import os,pickle
+
+###Used for Simulation Data
+def update_results(model_name,results,test_X,test_Y):
+    results[model_name]['auc'].append(compute_model_auc(eval(model_name),test_X,test_Y))
+    results[model_name]['bias'].append(compute_model_bias(eval(model_name),test_X,test_Y))    
+    results[model_name]['rig'].append(compute_model_rig(eval(model_name),test_X,test_Y))
+    results[model_name]['f1_score'].append(compute_model_f1(eval(model_name),test_X,test_Y))
+    
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-s', '--seed', default=16, type=int, help='Random seed for repro')
+    # Settings for true pclick distribution
+    parser.add_argument('--n_rnd', default=1000, type=int, help='Random Data Size')
+    parser.add_argument('--n_log', default=3500, type=int, help='Log Data Size')
+    parser.add_argument('--n_test', default=2000, type=int, help='Random Data Size')
+    parser.add_argument('--p', default=40, type=int, help='Feature Size')
+    parser.add_argument('--scenario_log', default=2, type=int, help='Log data Scenario')
+    parser.add_argument('--scenario_test', default=2, type=int, help='Testing data Scenario')
+    parser.add_argument('--r_log', default=0.7, type=float, help='Sampling Bias for Log data')
+    parser.add_argument('--r_test', default=0.7, type=float, help='Sampling Bias for Testing data')
+    parser.add_argument('--n_experiments', default=200, type=int, help='Number of Experiments')
+    args = parser.parse_args()
+    print('Arguments', args)
+    
+    model_list=['lr_model','lr_weight_model','gbdt_model','gbdt_weight_model','rndrf_model','cntrf_model','trtrf_model','combinerf_model','ctrf_model']
+    results={name:{metric:[] for metric in ['auc','f1_score','bias','rig']} for name in model_list}
+ 
+    if os.path.isdir('simu_results'): 
+        print('')
+    else:
+        os.mkdir("simu_results")
+        print("directory added")
+
+    path="simu_results/"
+    os.chdir(path)
+    result_name='_'.join(["results",str(args.n_log),str(args.p),str(args.scenario_test),str(int(100*args.r_test))])
+
+    for i in range(args.n_experiments):
+        ##Data Generating
+        print ('------------------------------')
+        print ('Run %d th experiments:' % (i+1))
+        random_data=simu_confounding_data(p=args.p,n=args.n_rnd,scenario=1,r=0.5)
+        log_data=simu_confounding_data(p=args.p,n=args.n_log,scenario=args.scenario_log,r=args.r_log)
+        testing_data=simu_confounding_data(p=args.p,n=args.n_test,scenario=args.scenario_test,r=args.r_test)
+        #IPW Calculation
+        ipw_weights=calculate_weight(log_data['X'],testing_data['X'])
+
+        ##Logistic Regression
+        lr_model,args.seed=train_lr_model(log_data['X'],log_data['Y'],seed=args.seed)
+        ##Logistic Regression with IPW adjustment
+        lr_weight_model,args.seed=train_lr_weight_model(log_data['X'],log_data['Y'],ipw_weights,seed=args.seed)
+        ##GBDT
+        gbdt_model,args.seed=train_gbdt_model(log_data['X'],log_data['Y'],seed=args.seed)
+        ##GBDT with IPW adjustment
+        gbdt_weight_model,args.seed=train_gbdt_weight_model(log_data['X'],log_data['Y'],ipw_weights,seed=args.seed)
+        #RF on randomized
+        rndrf_model,args.seed=train_rf(random_data['X'],random_data['Y'],seed=args.seed)
+        #RF on Log
+        cntrf_model,args.seed=train_rf(log_data['X'],log_data['Y'],seed=args.seed)
+        #Oracle RF
+        trtrf_model,args.seed=train_rf(testing_data['X'],testing_data['Y'],seed=args.seed)
+        #Combined RF
+        combinerf_model,args.seed=train_combine_rf(random_data['X'],random_data['Y'],log_data['X'],log_data['Y'],seed=args.seed)       
+        #CTRF
+        ctrf_model,args.seed=train_ctrf(random_data['X'],random_data['Y'],log_data['X'],log_data['Y'],rndrf_model,combine=1,seed=args.seed)
+
+        #Record Results
+        for model_name in model_list:
+            update_results(model_name,results,testing_data['X'],testing_data['Y'])
+
+        print ('Finish %d th experiments:' % (i+1))
+    
+    f = open(result_name+".pkl","wb")
+    pickle.dump(results,f)
+    f.close()
\ No newline at end of file
diff --git a/simulated_auction_reserve_varying.sh b/simulated_auction_reserve_varying.sh
new file mode 100644
index 0000000..7c58a66
--- /dev/null
+++ b/simulated_auction_reserve_varying.sh
@@ -0,0 +1,9 @@
+python run_auction.py --treatment_reserve 0.5 --n_auctions  25000 --n_rnd_auction 10000 --n_experiments 50 &
+python run_auction.py --treatment_reserve 0.55 --n_auctions  25000 --n_rnd_auction 10000 --n_experiments 50 &
+python run_auction.py --treatment_reserve 0.60 --n_auctions  25000 --n_rnd_auction 10000 --n_experiments 50 &
+python run_auction.py --treatment_reserve 0.65 --n_auctions 25000 --n_rnd_auction 10000 --n_experiments 50 &
+python run_auction.py --treatment_reserve 0.70 --n_auctions  25000 --n_rnd_auction 10000 --n_experiments 50 &
+python run_auction.py --treatment_reserve 0.75 --n_auctions  25000 --n_rnd_auction 10000 --n_experiments 50 &
+python run_auction.py --treatment_reserve 0.80 --n_auctions  25000 --n_rnd_auction 10000 --n_experiments 50 &
+python run_auction.py --treatment_reserve 0.85 --n_auctions 25000 --n_rnd_auction 10000 --n_experiments 50 &
+python run_auction.py --treatment_reserve 0.90 --n_auctions  25000 --n_rnd_auction 10000 --n_experiments 50 &
\ No newline at end of file
diff --git a/simulated_data_varying.sh b/simulated_data_varying.sh
new file mode 100644
index 0000000..b566048
--- /dev/null
+++ b/simulated_data_varying.sh
@@ -0,0 +1,27 @@
+python -W ignore run_simulation.py --r_test 0.1 --scenario_test 3 --p 20 &
+python -W ignore run_simulation.py --r_test 0.2 --scenario_test 3 --p 20 &
+python -W ignore run_simulation.py --r_test 0.3 --scenario_test 3 --p 20 &
+python -W ignore run_simulation.py --r_test 0.4 --scenario_test 3 --p 20 &
+python -W ignore run_simulation.py --r_test 0.5 --scenario_test 3 --p 20 &
+python -W ignore run_simulation.py --r_test 0.6 --scenario_test 3 --p 20 &
+python -W ignore run_simulation.py --r_test 0.7 --scenario_test 3 --p 20 &
+python -W ignore run_simulation.py --r_test 0.8 --scenario_test 3 --p 20 &
+python -W ignore run_simulation.py --r_test 0.9 --scenario_test 3 --p 20 & 
+python -W ignore run_simulation.py --r_test 0.1 --scenario_test 3 --p 80 &
+python -W ignore run_simulation.py --r_test 0.2 --scenario_test 3 --p 80 &
+python -W ignore run_simulation.py --r_test 0.3 --scenario_test 3 --p 80 &
+python -W ignore run_simulation.py --r_test 0.4 --scenario_test 3 --p 80 &
+python -W ignore run_simulation.py --r_test 0.5 --scenario_test 3 --p 80 &
+python -W ignore run_simulation.py --r_test 0.6 --scenario_test 3 --p 80 &
+python -W ignore run_simulation.py --r_test 0.7 --scenario_test 3 --p 80 &
+python -W ignore run_simulation.py --r_test 0.8 --scenario_test 3 --p 80 &
+python -W ignore run_simulation.py --r_test 0.9 --scenario_test 3 --p 80 & 
+python -W ignore run_simulation.py --r_test 0.1 --scenario_test 3  &
+python -W ignore run_simulation.py --r_test 0.2 --scenario_test 3  &
+python -W ignore run_simulation.py --r_test 0.3 --scenario_test 3  &
+python -W ignore run_simulation.py --r_test 0.4 --scenario_test 3  &
+python -W ignore run_simulation.py --r_test 0.5 --scenario_test 3  &
+python -W ignore run_simulation.py --r_test 0.6 --scenario_test 3  &
+python -W ignore run_simulation.py --r_test 0.7 --scenario_test 3  &
+python -W ignore run_simulation.py --r_test 0.8 --scenario_test 3  &
+python -W ignore run_simulation.py --r_test 0.9 --scenario_test 3  &