{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"##################################################################\n#                                                                #\n#             Wonsoek Bae. Chung-Ang Univ. 2019.05.24            #\n#                                                                #\n##################################################################\n\n#======================================================================================================================\n\nprint(\"0. 필요한 모듈 import\")\n\nimport numpy as np\nfrom sklearn.metrics import roc_curve, auc\n\"\"\"fpr, tpr, thresholds = roc_curve(y, model.decision_function(X)) \n입력값은 타켓 y 벡터와 변하는 판결값(함수의 반환값)\n출력값 false positive rate, true positive rate, 판별값들 묶음. 이것들로  plot 그리면 roc curve 그릴 수 있음\nauc : roc_curve 의 넓이\"\"\"\n\nimport pandas as pd\n\"1. pandas.dataframe : 데이터 집합형식-인덱스가 있는 다차원 배열(행렬)\"\nimport lightgbm as lgb\n\"\"\"1. lightgbm.Dataset : lgb 함수에 쓰이는 데이터 집합형식\n 2. lightgbm.cv : 입력값으로 Booster params 와, Data to be trained 를 받아서, 딕셔너리 타입의 Evaluation history 를 반환.\n 3. lightgbm.train : 입력값으로 Parameters for training 와 Data to be trained 를 받아서 학습된 부스터 모델을 반환\"\"\"\n\nfrom bayes_opt import BayesianOptimization\nfrom sklearn.metrics import roc_auc_score\nimport warnings\nimport time\nwarnings.filterwarnings(\"ignore\")\n\n# 필요한 함수들\n\"\"\"여기 함수는 CERN과 함께 이 프로젝트의 스폰서인 yandexdataschool 에서 제공\nhttps://github.com/yandexdataschool/flavours-of-physics-start/blob/master/evaluation.py\"\"\"\n\n\"코릴레이션 검증을 위한 함수\"\ndef __rolling_window(data, window_size):\n    \"\"\"\n    Rolling window: take window with definite size through the array\n\n    :param data: array-like\n    :param window_size: size\n    :return: the sequence of windows\n\n    Example: data = array(1, 2, 3, 4, 5, 6), window_size = 4\n        Then this function return array(array(1, 2, 3, 4), array(2, 3, 4, 5), array(3, 4, 5, 6))\n    \"\"\"\n    shape = data.shape[:-1] + (data.shape[-1] - window_size + 1, window_size)\n    strides = data.strides + (data.strides[-1],)\n    return np.lib.stride_tricks.as_strided(data, shape=shape, strides=strides)\n\ndef __cvm(subindices, total_events):\n    \"\"\"\n    Compute Cramer-von Mises metric.\n    Compared two distributions, where first is subset of second one.\n    Assuming that second is ordered by ascending\n\n    :param subindices: indices of events which will be associated with the first distribution\n    :param total_events: count of events in the second distribution\n    :return: cvm metric\n    \"\"\"\n    target_distribution = np.arange(1, total_events + 1, dtype='float') / total_events\n    subarray_distribution = np.cumsum(np.bincount(subindices, minlength=total_events), dtype='float')\n    subarray_distribution /= 1.0 * subarray_distribution[-1]\n    return np.mean((target_distribution - subarray_distribution) ** 2)\n\ndef compute_cvm(predictions, masses, n_neighbours=100, step=50):\n    \"\"\"\n    Computing Cramer-von Mises (cvm) metric on background events: take average of cvms calculated for each mass bin.\n    In each mass bin global prediction's cdf is compared to prediction's cdf in mass bin.\n\n    :param predictions: array-like, predictions\n    :param masses: array-like, in case of Kaggle tau23mu this is reconstructed mass\n    :param n_neighbours: count of neighbours for event to define mass bin\n    :param step: step through sorted mass-array to define next center of bin\n    :return: average cvm value\n    \"\"\"\n    predictions = np.array(predictions)\n    masses = np.array(masses)\n    assert len(predictions) == len(masses)\n\n    # First, reorder by masses\n    predictions = predictions[np.argsort(masses)]\n\n    # Second, replace probabilities with order of probability among other events\n    predictions = np.argsort(np.argsort(predictions, kind='mergesort'), kind='mergesort')\n\n    # Now, each window forms a group, and we can compute contribution of each group to CvM\n    cvms = []\n    for window in __rolling_window(predictions, window_size=n_neighbours)[::step]:\n        cvms.append(__cvm(subindices=window, total_events=len(predictions)))\n    return np.mean(cvms)\n\n\n\"어그리먼트 검증을 위한 함수\"\ndef __roc_curve_splitted(data_zero, data_one, sample_weights_zero, sample_weights_one):\n    \"\"\"\n    Compute roc curve\n\n    :param data_zero: 0-labeled data\n    :param data_one:  1-labeled data\n    :param sample_weights_zero: weights for 0-labeled data\n    :param sample_weights_one:  weights for 1-labeled data\n    :return: roc curve\n    \"\"\"\n    labels = [0] * len(data_zero) + [1] * len(data_one)\n    weights = np.concatenate([sample_weights_zero, sample_weights_one])\n    data_all = np.concatenate([data_zero, data_one])\n    fpr, tpr, _ = roc_curve(labels, data_all, sample_weight=weights)\n    return fpr, tpr\n\ndef compute_ks(data_prediction, mc_prediction, weights_data, weights_mc):\n    \"\"\"\n    Compute Kolmogorov-Smirnov (ks) distance between real data predictions cdf and Monte Carlo one.\n\n    :param data_prediction: array-like, real data predictions\n    :param mc_prediction: array-like, Monte Carlo data predictions\n    :param weights_data: array-like, real data weights\n    :param weights_mc: array-like, Monte Carlo weights\n    :return: ks value\n    \"\"\"\n    assert len(data_prediction) == len(weights_data), 'Data length and weight one must be the same'\n    assert len(mc_prediction) == len(weights_mc), 'Data length and weight one must be the same'\n\n    data_prediction, mc_prediction = np.array(data_prediction), np.array(mc_prediction)\n    weights_data, weights_mc = np.array(weights_data), np.array(weights_mc)\n\n    assert np.all(data_prediction >= 0.) and np.all(data_prediction <= 1.), 'Data predictions are out of range [0, 1]'\n    assert np.all(mc_prediction >= 0.) and np.all(mc_prediction <= 1.), 'MC predictions are out of range [0, 1]'\n\n    weights_data /= np.sum(weights_data)\n    weights_mc /= np.sum(weights_mc)\n\n    fpr, tpr = __roc_curve_splitted(data_prediction, mc_prediction, weights_data, weights_mc)\n\n    Dnm = np.max(np.abs(fpr - tpr))\n    return Dnm\n\n\"평가를 위한 함수\"\ndef roc_auc_truncated(labels, predictions, tpr_thresholds=(0.2, 0.4, 0.6, 0.8), roc_weights=(4, 3, 2, 1, 0)):\n    \"\"\"\n    These weights were chosen to match the evaluation methodology used by CERN scientists.\n    Note that the weighted AUC is calculated\n    only for events (simulated signal events for tau->µµµ and real background events for tau->µµµ) with min_ANNmuon > 0.4\n    \"\"\"\n\n    \"\"\"\n    Compute weighted area under ROC curve.\n\n    :param labels: array-like, true labels\n    :param predictions: array-like, predictions\n    :param tpr_thresholds: array-like, true positive rate thresholds delimiting the ROC segments\n    :param roc_weights: array-like, weights for true positive rate segments\n    :return: weighted AUC\n    \"\"\"\n    assert np.all(predictions >= 0.) and np.all(predictions <= 1.), 'Data predictions are out of range [0, 1]'\n    assert len(tpr_thresholds) + 1 == len(roc_weights), 'Incompatible lengths of thresholds and weights'\n    fpr, tpr, _ = roc_curve(labels, predictions)\n    area = 0.\n    tpr_thresholds = [0.] + list(tpr_thresholds) + [1.]\n    for index in range(1, len(tpr_thresholds)):\n        tpr_cut = np.minimum(tpr, tpr_thresholds[index])\n        tpr_previous = np.minimum(tpr, tpr_thresholds[index - 1])\n        area += roc_weights[index - 1] * (auc(fpr, tpr_cut, reorder=True) - auc(fpr, tpr_previous, reorder=True))\n    tpr_thresholds = np.array(tpr_thresholds)\n    # roc auc normalization to be 1 for an ideal classifier\n    area /= np.sum((tpr_thresholds[1:] - tpr_thresholds[:-1]) * np.array(roc_weights))\n    return area","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"code","source":"print(\"1. 데이터 가져오기\")\n\ntrain_set = pd.read_csv('../input/training.csv', index_col='id') # [67553 rows x 50 columns]\ntest = pd.read_csv('../input/test.csv', index_col='id') # [855819 rows x 46 columns]\ncheck_agreement = pd.read_csv('../input/check_agreement.csv', index_col='id') # [331147 rows x 48 columns]\ncheck_correlation = pd.read_csv('../input/check_correlation.csv', index_col='id') # [5514 rows x 47 columns]\n# id를 행 인덱스(세로축 레이블)로 .csv 데이터 파일을 끌어온다.","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"2. 데이터 전처리\")\n\n# 2.1 feature(레이블, 대게 물리량) engineering : 특성조합\n\ndef add_features(df):\n    df['flight_dist_sig2'] = (df['FlightDistance'] / df['FlightDistanceError']) ** 2\n    df['p_track_Chi2Dof_MAX'] = df.loc[:, ['p0_track_Chi2Dof', 'p1_track_Chi2Dof', 'p2_track_Chi2Dof']].max(axis=1)\n    df['NEW_FD_SUMP'] =df['FlightDistance']/(df['p0_p']+df['p1_p']+df['p2_p'])\n    df['NEW5_lt']=df['LifeTime']*(df['p0_IP']+df['p1_IP']+df['p2_IP'])/3\n    df['flight_dist_sig'] = df['FlightDistance']/df['FlightDistanceError']\n    df['NEW_IP_dira'] = df['IP']*df['dira']\n    # features from phunter\n    df['p0p2_ip_ratio']=df['IP']/df['IP_p0p2']\n    df['p1p2_ip_ratio']=df['IP']/df['IP_p1p2']\n    df['DCA_MAX'] = df.loc[:, ['DOCAone', 'DOCAtwo', 'DOCAthree']].max(axis=1)\n    df['iso_bdt_min'] = df.loc[:, ['p0_IsoBDT', 'p1_IsoBDT', 'p2_IsoBDT']].min(axis=1)\n    df['iso_min'] = df.loc[:, ['isolationa', 'isolationb', 'isolationc','isolationd', 'isolatione', 'isolationf']].min(axis=1)\n    # features from UGBC GS\n    df['NEW_iso_abc'] = df['isolationa']*df['isolationb']*df['isolationc']\n    df['NEW_iso_def'] = df['isolationd']*df['isolatione']*df['isolationf']\n    df['NEW_pN_IP'] = df['p0_IP']+df['p1_IP']+df['p2_IP']\n    df['NEW_pN_p']  = df['p0_p']+df['p1_p']+df['p2_p']\n    df['NEW_IP_pNpN'] = df['IP_p0p2']*df['IP_p1p2']\n    df['NEW_pN_IPSig'] = df['p0_IPSig']+df['p1_IPSig']+df['p2_IPSig']\n    df['NEW_FD_LT']=df['FlightDistance']/df['LifeTime']\n    return df\n\n# 2.2 새 레이블을 이용해서 데이터 가공\n\ntrain_add = add_features(train_set) # [67553 rows x 68 columns]\ntest_add = add_features(test) # [855819 rows x 64 columns]\nagr_add = add_features(check_agreement) # [331147 rows x 66 columns]\ncor_add = add_features(check_correlation) # [5514 rows x 65 columns]\n\n\"기존 피쳐 지우기\" # idea from UGBC GS\nfilter_out = ['id', 'min_ANNmuon', 'production', 'mass', 'signal', 'SPDhits','CDF1', 'CDF2', 'CDF3',\n              'isolationb', 'isolationc','p0_pt', 'p1_pt', 'p2_pt', 'p0_p', 'p1_p', 'p2_p', 'p0_eta', 'p1_eta', 'p2_eta',\n              'isolationa', 'isolationb', 'isolationc', 'isolationd', 'isolatione', 'isolationf',\n              'p0_IsoBDT', 'p1_IsoBDT', 'p2_IsoBDT', 'p0_IP', 'p1_IP', 'p2_IP', 'IP_p0p2', 'IP_p1p2',\n              'p0_track_Chi2Dof', 'p1_track_Chi2Dof', 'p2_track_Chi2Dof', 'p0_IPSig', 'p1_IPSig', 'p2_IPSig',\n              'DOCAone', 'DOCAtwo', 'DOCAthree']\n\n\"새 레이블에만 있고 옛 레이블에 없는 것만 모아서 col 이라는 레이블 만들기\"\ncol = list(f for f in train_add.columns if f not in filter_out) # 28 columns\n# col에 id없어도 id는 이미 행레이블로 들어갔음.","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"3 Bayesian Optimization으로 최적의 Hyper parameters 찾기\")\n\ndef bayes_parameter_opt_lgb(init_round=20, opt_round=30, n_folds=5, random_seed=42, n_estimators=50000, learning_rate=0.0001):\n    # prepare data\n    train_data = lgb.Dataset(train_add[col], train_set['signal']) # 두 번째 변수가 True인 것만 뽑음\n    # parameters\n    def lgb_eval(num_leaves, feature_fraction, bagging_fraction, max_depth, lambda_l1, lambda_l2, min_split_gain, min_child_weight):\n        params = {'application':'binary','num_iterations': n_estimators, 'learning_rate':learning_rate, 'early_stopping_round':100, 'metric':'auc'}\n        params[\"num_leaves\"] = int(round(num_leaves))\n        params['feature_fraction'] = max(min(feature_fraction, 1), 0)\n        params['bagging_fraction'] = max(min(bagging_fraction, 1), 0)\n        params['max_depth'] = int(round(max_depth))\n        params['lambda_l1'] = max(lambda_l1, 0)\n        params['lambda_l2'] = max(lambda_l2, 0)\n        params['min_split_gain'] = min_split_gain\n        params['min_child_weight'] = min_child_weight\n        cv_result = lgb.cv(params, train_data, nfold=n_folds, seed=random_seed, stratified=True, metrics=['auc'])     \n        \n        max_value_niter = np.argmax(cv_result['auc-mean']) # CV의 리턴값이 딕셔너리로 나올텐데, 그 중에서 auc-mean 키의 밸류들 중에서 최대값에 대응되는 인덱스\n        print('Best number of iterations: {}'.format(max_value_niter)) # Best number of iterations: Best number of iterations: 49999\n        max_value_score = cv_result['auc-mean'][max_value_niter] # cv['auc-mean'] 중에서 그 인덱스의 값 = 즉 최대값\n        print('Best CV score: {}'.format(max_value_score)) # Best CV score: 0.9407012723236857\n        \n        return max(cv_result['auc-mean'])\n    \n    # range \n    lgbBO = BayesianOptimization(lgb_eval, {'num_leaves': (20, 614),# num_leaves = 2^maxmax_depth * 60% \n                                            'feature_fraction': (0.8, 0.9),\n                                            'bagging_fraction': (0.8, 1),\n                                            'max_depth': (5, 10),\n                                            'lambda_l1': (0, 5),\n                                            'lambda_l2': (0, 3),\n                                            'min_split_gain': (0.001, 0.1),\n                                            'min_child_weight': (5, 50)}, random_state=0)\n        \n    # optimize\n    lgbBO.maximize(init_points=init_round, n_iter=opt_round)\n    \n    return lgbBO.res\n\nopt_params = bayes_parameter_opt_lgb(init_round=1, opt_round=1, n_folds=5, random_seed=42, n_estimators=10, learning_rate=0.0001)\nprint(opt_params)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"4. lgb.cv랑 lgb.train 으로 가장 좋은 부스터 모델 만들기\")\n\n# 4.1 Light gradient boosting machine 에 들어갈 3에서 찾은 최적의 하이퍼 파라메터\n\n# Best CV score: 0.934151783702392\n#|  31       |  0.9342   |  1.0      |  0.8      |  4.648    |  0.0      |  9.99     |  5.0      |  0.04201  |  45.0\n# 옵티마이제이션으로 최적의 하이퍼 파라미터를 찾은 뒤, 다른 조건으로 여러 번 옵티마이제이션을 작동 시켰기 때문에\n# 이전에 찾았던 최적의 하이퍼 파라미터를 수동으로 입력해서 사용.\nparams = {'application':'binary','num_iterations': 10, 'learning_rate':0.0001, 'metric':'auc'} # params이라는 딕셔너리 만들기\nparams[\"num_leaves\"] = int(round(45.0 )) # number of leaves in one tree\nparams['feature_fraction'] = 0.8 # 0~1 사이의 값. 0.8이면, 각각의 tree에 80%의 features 사용. 학습 속도 높이고 over-fitting 막는데 씀\nparams['bagging_fraction'] = 1.0 # 0~1 사이의 값. 데이터를 랜덤 추출하게 함. 학습 속도 높이고 over-fitting 막는데 씀\nparams['max_depth'] = int(round(9.99)) # Limit the max depth for tree model. 데이터 작을 때 over-fitting 막는데 씀\nparams['lambda_l1'] = 4.638 # regularization 해서 over-fitting 막는데 씀\nparams['lambda_l2'] = 0.0 # regularization 해서 over-fitting 막는데 씀\nparams['min_split_gain'] = 0.04201 # split 을 만들기 위한 최소 gain. 트리 안에서 useful split 수 조절\nparams['min_child_weight'] = 5.0 # leaf에서 필요한 instance weight (hessian)의 최소 합 .\n\n\"\"\"application':'binary = 참/거짓 이진판별, num_iterations = number of boosting iterations/trees \nmetrics (string or list of strings) – Evaluation metrics to be watched in CV.\nmetrics ex)\n             실제 참   실제 거짓\n  예측 참    10000개     10개\n  예측 거짓     10개    10000개\n\"\"\"\n\n\n\"가공된 훈련 데이터를 lgb 함수에 넣을 수 있게 lgb 데이터 형식으로 만들어 주기\"\ntrain = lgb.Dataset(train_add[col], train_set['signal']) # 두 번째 변수가 True인 것만 뽑음\n\"train = train <lightgbm.basic.Dataset object at 0x000002940095D438> 출력하면 그냥 28행 레이블만 나옴\"\n\n# 4.2 Best boosting number 구하기 : cross validation\ncv = lgb.cv(params, train, nfold=5, seed=42, stratified=True, verbose_eval =100, metrics=['auc'])\n\"\"\"입력값으로 Booster params 와, Data to be trained 를 받아서, 딕셔너리 타입의 Evaluation history 를 반환.\n  하리퍼 파라메터를 가지고, 데이터를 nfold 개수로 나눈 다음, 파라메터를 바꿔가면서 교차 평가하면 평가 (auc) 점수가 오락가락 함 \n  nfold에서 교차해본 평균값이 auc-mean\n  \n  입력값 :\n  params, Dataset, \n  nfold : dataset 을 몇 개의 균등한 크기의 subsamples 로 나눌 건가.(train data가 60000개 안 되서 5개로 나눔)\n  seed - fold 별로 데이터를 구성을 random 하게 나눠주는 것.\n  stratified - sampling of folds should be stratified by the values of outcome labels.)\n  verbose_eval - the eval metric on the valid set is printed at every verbose_eval boosting stage. 100번째마다 출력 \n  ex) [100]\tcv_agg's auc: 0.914896 + 0.0033455\n      [200]\tcv_agg's auc: 0.915306 + 0.00328177\n  \n  반환값 :\n  딕셔너리 타입으로 Evaluation history 를 반환. \n  그 중 최대 auc값과 그 때의 인덱스를 찾을 것.\n  {'auc-mean': [0.900584369149491, 0.9068785796788049, ... 0.9183340848474986, 0.9183348220084587]\n  , 'auc-stdv': [0.003673678624173514, 0.00354709321385367, ... 0.003214837891260049, 0.003214937592912147]}\"\"\"\n\nbest_niter = np.argmax(cv['auc-mean']) # CV의 리턴값이 딕셔너리로 나올텐데, 그 중에서 auc-mean 키의 밸류들 중에서 최대값에 대응되는 인덱스\nprint('Best number of iterations: {}'.format(best_niter)) # Best number of iterations: Best number of iterations: 49999\nbest_score = cv['auc-mean'][best_niter] # cv['auc-mean'] 중에서 그 인덱스의 값 = 즉 최대값\nprint('Best CV score: {}'.format(best_score)) # Best CV score: 0.9407012723236857\n\n# 4.3 분류기 학습시키기\nclf = lgb.train(params, train, num_boost_round=best_niter)\n# 여기가 트레인 돌리기.  이게 우리가 만든 분류기. clf = classifier\n#  입력값 : \n#      params (dict), train_set (Dataset), num_boost_round (int) – Number of boosting iterations\n#      \"num_boost_round = Number of boosting iterations. 이게 몇 학습기를 몇 개 연결시켜줄거냐는 것?\"\n#      \"light bgm 에서 학습기 연결 개수가 tree 개수?\"\n#  반환값 : \n#      booster – The trained Booster model. ex N번 부스터 해서 학습된 부스터 모델을 반환","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"5. 입자물리학 실험에서 요구하는 검증 조건 확인하기\")\n\n# 5.1 어그리먼트 테스트\n\"\"\"분류기가 학습하는 훈련 데이터에는, 실제 데이터와 시뮬레이션(mc)으로 만들어 진 데이터가 섞여 있음.\n 따라서, 시뮬레이션 데이터가 충분히 정확하기 않게 만들어지고, 분류기가 그것에 대해 학습할 수 있음.\n 그래서, τ → 3μ 붕괴 데이터와 유사한 Ds → φπ 붕괴의 실제 데이터와 시뮬레이션 데이터를 분류시켜보고, \n 그것들이 τ → 3μ 붕괴가 아님을 충분히 구별해야, 학습된 분류기를 신뢰할 수 있다. \n 이를 Kolmogorov-Smirnov (KS) 테스트라 한다. * https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test\n \n 여기에서는 KS 값이 0.09보다 작아야 한다. 이것을 검증하기 위한 함수(코드)는 CERN에서 제공 : compute_ks, roc_curve_splitted \"\"\"\n\nagreement_probs = clf.predict(agr_add[col])\n\"가공된 어그리먼트 데이터를 학습된 분류기 모델에 넣어서 (참/거짓)에 대한 prediction 값을 출력 \"\n\"[0.60230961 0.579581   0.60094461 ... 0.62754844 0.6404412  0.6404412 ]\"\n\nks = compute_ks(agreement_probs[check_agreement['signal'].values == 0], agreement_probs[check_agreement['signal'].values == 1],\ncheck_agreement[check_agreement['signal'] == 0]['weight'].values, check_agreement[check_agreement['signal'] == 1]['weight'].values)\n\"\"\" 인자 1: data_prediction = signal 값이 0인 것들의 agreement_probs 정보, \n인자 2: mc_prediction = signal 값이 0인 것들의 agreement_probs 정보\n인자 3 :weights_data, = signal 값이 0인 것들의 weight 정보가 id : 값 딕셔너리 형태로 출력되는데, 거기서 weight 값,\n인자 4 : weights_mc = signal 값이 1인 것들의 weight 정보가 id : 값 딕셔너리 형태로 출력되는데, 거기서 weight 값\"\"\"\nprint ('KS metric : ', ks, ks < 0.09) # KS metric 0.05148947099501236 True\n\n# 5.2 코릴레이션 테스트\n\"\"\"주어진 데이터의 질량은 이미 τ → 3μ 붕괴 후보로서 가능한 값들이다. 따라서 만들어진 모델을 질량과 입자 예측치 간의 \n상관관계가 없어야 한다. (사실 입자 실험에서 질량 값은 적당한 추정치로만 얻어지며, 과학자들이 모델을 만들 때 정확하게 \n신뢰해서는 안 되는 값이다.) Test 데이터에는 잘량 값은 없지만, 다른 물리량들 관계 속에 숨겨진 질량 값을 를 사용하여 \n질량 값과 입자 예측치 간의 상관관계가 있는 지 검증한다. 이를 Cramer-von Mises (cvm) 테스트라 한다. \na) 전체 데이터 집합에 대한 제출에서 예측 된 값\nb) 전체 질량 범위를 따라 롤링 윈도우 방식으로 특정 질량 영역 내의 예측 된 값.\n비교 하고 둘의 평균 값을 반환. 이 CVM 값이 0.02 보다 작아야 한다. \n이것을 검증하기 위한 함수(코드)는 CERN에서 제공 : rolling_window, cvm, compute_cvm\"\"\"\n\n\ncorrelation_probs = clf.predict(cor_add[col])\n\"가공된 코릴레이션 데이터를 학습된 분류기 모델에 넣어서 (참/거짓)에 대한 prediction 값을 출력\"\ncvm = compute_cvm(correlation_probs, check_correlation['mass'])\n\"인자 1: 확률, 인자 2: 질량\"\nprint ('CVM metric : ', cvm, cvm < 0.002) # CvM metric 0.0010729350995806089 True","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"6. 결과 파일 만들기(각각의 id에 대한 참 거짓 prediction 값)\")\n\n#if ks < 0.09 and cvm < 0.002:\nresult = pd.DataFrame({'id': test.index})\nresult['prediction'] = clf.predict(test_add[col])\n\"가공된 test 데이터를 학습된 분류기 모델에 넣어서 (참/거짓)에 대한 prediction 값을 출력\"\nresult.to_csv('result.csv', index=False, sep=',') # 캐글에 올릴 결과 파일.\n\n#result\n#id       index\n#0       14711831\n#1       16316387\n \n#855818  12395271\n    \n#result['prediction'] \n#0         0.542113\n#1         0.542876\n    \n#855818    0.597033\n    \n#만들어서 둘이 합치기","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"7. 테스트 데이터의 정답은 CERN에서 제공하지 않으며, 제출한 파일에 대해 평가만 해준다. \\n 따라서, 위에서 얻은 결과에 대한 roc curve/precision-recall-threshold 그래프는 얻을 수 없다. \\n 따라서, 훈련 데이터를 8:2로 나눠서 80%의 데이터에 대해 다시 훈련시키고, 20%의 데이터에 대해 그것들을 평가해본다.\")\n\nfrom sklearn.model_selection import train_test_split\n%matplotlib inline\nimport matplotlib\nimport matplotlib.pyplot as plt\nfrom sklearn.metrics import precision_recall_curve\nfrom sklearn.metrics import roc_curve\n\n\n# 1. 데이터 나누기\ntrain_set1,train_set2 = train_test_split(train_set, test_size=0.2, random_state=42)\ntrain_add1 = add_features(train_set1)\ntrain_add2 = add_features(train_set2)\n\n# 2. CV 돌리기\n\ntrain = lgb.Dataset(train_add1[col], train_set1['signal']) # 두 번째 변수가 True인 것만 뽑음\ncv = lgb.cv(params, train, nfold=5, seed=42, stratified=True, verbose_eval =100, metrics=['auc'])\n\nbest_niter = np.argmax(cv['auc-mean']) # CV의 리턴값이 딕셔너리로 나올텐데, 그 중에서 auc-mean 키의 밸류들 중에서 최대값에 대응되는 인덱스\nprint('Best number of iterations: {}'.format(best_niter)) # Best number of iterations: Best number of iterations: 49999\nbest_score = cv['auc-mean'][best_niter] # cv['auc-mean'] 중에서 그 인덱스의 값 = 즉 최대값\nprint('Best CV score: {}'.format(best_score)) # Best CV score: 0.9407012723236857\n\n# 3. 분류기 학습시키기\nclf = lgb.train(params, train, num_boost_round=best_niter)\n\n# 4. 분류하기\ntrain_eval = train_add2\ntrain_eval_probs = clf.predict(train_eval[col]) # test_add[col] train_add[col]\n\n# 5. roc 값 구하기\nAUC = roc_auc_score(train_eval['signal'], train_eval_probs)  # 0/1 값(실제 참/거짓) 과 (예측 참/거짓)를 함수에 넣음\nprint('AUC : ', AUC)\n\n# 6. roc curve 구하기\nprecisions, recalls, thresholds = precision_recall_curve(train_eval['signal'], train_eval_probs)\n\ndef plot_precision_recall_vs_threshold(precisions, recalls, thresholds):\n    plt.plot(thresholds, precisions[:-1], \"b--\", label=\"Precision\")\n    plt.plot(thresholds, recalls[:-1], \"g-\", label=\"Recall\")\n    plt.xlabel(\"Threshold\")\n    plt.legend(loc=\"upper left\")\n    plt.ylim([0, 1])\n    \nplot_precision_recall_vs_threshold(precisions, recalls, thresholds)\nplt.show()\n\n#7. precision-recall-threshold 그래프\nfpr, tpr, thresholds = roc_curve(train_eval['signal'], train_eval_probs)\n\ndef plot_roc_curve(fpr, tpr, label=None):\n    plt.plot(fpr, tpr, linewidth=2, label=label)\n    plt.plot([0, 1], [0, 1], 'k--')\n    plt.axis([0, 1, 0, 1])\n    plt.xlabel('False Positive Rate')\n    plt.ylabel('True Positive Rate')\n\nplot_roc_curve(fpr, tpr)\nplt.show()","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}