{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom IPython.display import Markdown\nimport gc\nimport os\nclass _TFColor(object):\n    \"\"\"Enum of colors used in TF docs.\"\"\"\n    red = '#F15854'\n    blue = '#5DA5DA'\n    orange = '#FAA43A'\n    green = '#60BD68'\n    pink = '#F17CB0'\n    brown = '#B2912F'\n    purple = '#B276B2'\n    yellow = '#DECF3F'\n    gray = '#4D4D4D'\n    def __getitem__(self, i):\n        return [\n            self.red,\n            self.orange,\n            self.green,\n            self.blue,\n            self.pink,\n            self.brown,\n            self.purple,\n            self.yellow,\n            self.gray,\n        ][i % 9]\nTFColor = _TFColor()\n\n# using and inspired by https://www.kaggle.com/code/zhukovoleksiy/detect-sleep-states-starter-notebook-ensemble","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-09-14T18:37:29.194089Z","iopub.execute_input":"2023-09-14T18:37:29.194622Z","iopub.status.idle":"2023-09-14T18:37:29.216026Z","shell.execute_reply.started":"2023-09-14T18:37:29.194579Z","shell.execute_reply":"2023-09-14T18:37:29.214648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://www.kaggle.com/code/metric/event-detection-ap?scriptVersionId=141456876&cellId=1\n\"\"\"Event Detection Average Precision\n\nAn average precision metric for event detection in time series and\nvideo.\n\n\"\"\"\n\nimport numpy as np\nimport pandas as pd\nimport pandas.api.types\nfrom typing import Dict, List, Tuple\n\n\nclass ParticipantVisibleError(Exception):\n    pass\n\n\n# Set some placeholders for global parameters\nseries_id_column_name = None\ntime_column_name = None\nevent_column_name = None\nscore_column_name = None\nuse_scoring_intervals = None\n\n\ndef score(\n        solution: pd.DataFrame,\n        submission: pd.DataFrame,\n        tolerances: Dict[str, List[float]],\n        series_id_column_name: str,\n        time_column_name: str,\n        event_column_name: str,\n        score_column_name: str,\n        use_scoring_intervals: bool = False,\n) -> float:\n    \"\"\"Event Detection Average Precision, an AUCPR metric for event detection in\n    time series and video.\n\n    This metric is similar to IOU-threshold average precision metrics commonly\n    used in object detection. For events occuring in time series, we replace the\n    IOU threshold with a time tolerance.\n\n    Submissions are evaluated on the average precision of detected events,\n    averaged over timestamp error tolerance thresholds, averaged over event\n    classes.\n\n    Detections are matched to ground-truth events within error tolerances, with\n    ambiguities resolved in order of decreasing confidence.\n\n    Detailed Description\n    --------------------\n    Evaluation proceeds in four steps:\n\n    1. Selection - (optional) Predictions not within a series' scoring\n    intervals are dropped.\n    2. Assignment - Predicted events are matched with ground-truth events.\n    3. Scoring - Each group of predictions is scored against its corresponding\n    group of ground-truth events via Average Precision.\n    4. Reduction - The multiple AP scores are averaged to produce a single\n    overall score.\n\n    Selection\n\n    With each series there may be a defined set of scoring intervals giving the\n    intervals of time over which zero or more ground-truth events might be\n    annotated in that series. A prediction will be evaluated only if it falls\n    within a scoring interval. These scoring intervals can be chosen to improve\n    the fairness of evaluation by, for instance, ignoring edge-cases or\n    ambiguous events.\n\n    It is recommended that, if used, scoring intervals be provided for training\n    data but not test data.\n\n    Assignment\n\n    For each set of predictions and ground-truths within the same `event x\n    tolerance x series_id` group, we match each ground-truth to the\n    highest-confidence unmatched prediction occurring within the allowed\n    tolerance.\n\n    Some ground-truths may not be matched to a prediction and some predictions\n    may not be matched to a ground-truth. They will still be accounted for in\n    the scoring, however.\n\n    Scoring\n\n    Collecting the events within each `series_id`, we compute an Average\n    Precision score for each `event x tolerance` group. The average precision\n    score is the area under the (step-wise) precision-recall curve generated by\n    decreasing confidence score thresholds over the predictions. In this\n    calculation, matched predictions over the threshold are scored as TP and\n    unmatched predictions as FP. Unmatched ground-truths are scored as FN.\n\n    Reduction\n\n    The final score is the average of the above AP scores, first averaged over\n    tolerance, then over event.\n\n    Parameters\n    ----------\n    solution : pd.DataFrame, with columns:\n\n        `series_id_column_name` identifier for each time series\n\n        `time_column_name` the time of occurence for each event as a numeric type\n\n        `event_column_name` class label for each event\n\n        The solution contains the time of occurence of one or more types of\n        event within one or more time series. The metric expects the solution to\n        contain the same event types as those given in `tolerances`.\n\n        When `use_scoring_intervals == True`, you may include `start` and `end`\n        events to delimit intervals within which detections will be scored.\n        Detected events (from the user submission) outside of these events will\n        be ignored.\n\n    submission : pd.DataFrame, with columns as above and in addition:\n\n        `score_column_name` the predicted confidence score for the detected event\n\n    tolerances : Dict[str, List[float]]\n\n        Maps each event class to a list of timestamp tolerances used\n        for matching detections to ground-truth events.\n\n    use_scoring_intervals: bool, default False\n\n        Whether to ignore predicted events outside intervals delimited\n        by `'start'` and `'end'` events in the solution. When `False`,\n        the solution should not include `'start'` and `'end'` events.\n        See the examples for illustration.\n\n    Returns\n    -------\n    event_detection_ap : float\n        The mean average precision of the detected events.\n\n    Examples\n    --------\n    Detecting `'pass'` events in football:\n    >>> column_names = {\n    ...     'series_id_column_name': 'video_id',\n    ...     'time_column_name': 'time',\n    ...     'event_column_name': 'event',\n    ...     'score_column_name': 'score',\n    ... }\n    >>> tolerances = {'pass': [1.0]}\n    >>> solution = pd.DataFrame({\n    ...     'video_id': ['a', 'a'],\n    ...     'event': ['pass', 'pass'],\n    ...     'time': [0, 15],\n    ... })\n    >>> submission = pd.DataFrame({\n    ...     'video_id': ['a', 'a', 'a'],\n    ...     'event': ['pass', 'pass', 'pass'],\n    ...     'score': [1.0, 0.5, 1.0],\n    ...     'time': [0, 10, 14.5],\n    ... })\n    >>> score(solution, submission, tolerances, **column_names)\n    1.0\n\n    Increasing the confidence score of the false detection above the true\n    detections decreases the AP.\n    >>> submission.loc[1, 'score'] = 1.5\n    >>> score(solution, submission, tolerances, **column_names)\n    0.6666666666666666...\n\n    Likewise, decreasing the confidence score of a true detection below the\n    false detection also decreases the AP.\n    >>> submission.loc[1, 'score'] = 0.5  # reset\n    >>> submission.loc[0, 'score'] = 0.0\n    >>> score(solution, submission, tolerances, **column_names)\n    0.8333333333333333...\n\n    We average AP scores over tolerances. Previously, the detection at 14.5\n    would match, but adding smaller tolerances gives AP scores where it does\n    not match. This results in both a FN, since the ground-truth wasn't\n    detected, and a FP, since the detected event matches no ground-truth.\n    >>> tolerances = {'pass': [0.1, 0.2, 1.0]}\n    >>> score(solution, submission, tolerances, **column_names)\n    0.3888888888888888...\n\n    We also average over time series and over event classes.\n    >>> tolerances = {'pass': [0.5, 1.0], 'challenge': [0.25, 0.50]}\n    >>> solution = pd.DataFrame({\n    ...     'video_id': ['a', 'a', 'b'],\n    ...     'event': ['pass', 'challenge', 'pass'],\n    ...     'time': [0, 15, 0],  # restart time for new time series b\n    ... })\n    >>> submission = pd.DataFrame({\n    ...     'video_id': ['a', 'a', 'b'],\n    ...     'event': ['pass', 'challenge', 'pass'],\n    ...     'score': [1.0, 0.5, 1.0],\n    ...     'time': [0, 15, 0],\n    ... })\n    >>> score(solution, submission, tolerances, **column_names)\n    1.0\n\n    By adding scoring intervals to the solution, we may choose to ignore\n    detections outside of those intervals.\n    >>> tolerances = {'pass': [1.0]}\n    >>> solution = pd.DataFrame({\n    ...     'video_id': ['a', 'a', 'a', 'a'],\n    ...     'event': ['start', 'pass', 'pass', 'end'],\n    ...     'time': [0, 10, 20, 30],\n    ... })\n    >>> submission = pd.DataFrame({\n    ...     'video_id': ['a', 'a', 'a'],\n    ...     'event': ['pass', 'pass', 'pass'],\n    ...     'score': [1.0, 1.0, 1.0],\n    ...     'time': [10, 20, 40],\n    ... })\n    >>> score(solution, submission, tolerances, **column_names, use_scoring_intervals=True)\n    1.0\n\n    \"\"\"\n    # Validate metric parameters\n    assert len(tolerances) > 0, \"Events must have defined tolerances.\"\n    assert set(tolerances.keys()) == set(solution[event_column_name]).difference({'start', 'end'}),\\\n        (f\"Solution column {event_column_name} must contain the same events \"\n         \"as defined in tolerances.\")\n    assert pd.api.types.is_numeric_dtype(solution[time_column_name]),\\\n        f\"Solution column {time_column_name} must be of numeric type.\"\n\n    # Validate submission format\n    for column_name in [\n        series_id_column_name,\n        time_column_name,\n        event_column_name,\n        score_column_name,\n    ]:\n        if column_name not in submission.columns:\n            raise ParticipantVisibleError(f\"Submission must have column '{target_name}'.\")\n\n    if not pd.api.types.is_numeric_dtype(submission[time_column_name]):\n        raise ParticipantVisibleError(\n            f\"Submission column '{time_column_name}' must be of numeric type.\"\n        )\n    if not pd.api.types.is_numeric_dtype(submission[score_column_name]):\n        raise ParticipantVisibleError(\n            f\"Submission column '{score_column_name}' must be of numeric type.\"\n        )\n\n    # Set these globally to avoid passing around a bunch of arguments\n    globals()['series_id_column_name'] = series_id_column_name\n    globals()['time_column_name'] = time_column_name\n    globals()['event_column_name'] = event_column_name\n    globals()['score_column_name'] = score_column_name\n    globals()['use_scoring_intervals'] = use_scoring_intervals\n\n    return event_detection_ap(solution, submission, tolerances)\n\n\ndef filter_detections(\n        detections: pd.DataFrame, intervals: pd.DataFrame\n) -> pd.DataFrame:\n    \"\"\"Drop detections not inside a scoring interval.\"\"\"\n    detection_time = detections.loc[:, time_column_name].sort_values().to_numpy()\n    intervals = intervals.to_numpy()\n    is_scored = np.full_like(detection_time, False, dtype=bool)\n\n    i, j = 0, 0\n    while i < len(detection_time) and j < len(intervals):\n        time = detection_time[i]\n        int_ = intervals[j]\n\n        # If the detection is prior in time to the interval, go to the next detection.\n        if time < int_.left:\n            i += 1\n        # If the detection is inside the interval, keep it and go to the next detection.\n        elif time in int_:\n            is_scored[i] = True\n            i += 1\n        # If the detection is later in time, go to the next interval.\n        else:\n            j += 1\n\n    return detections.loc[is_scored].reset_index(drop=True)\n\n\ndef match_detections(\n        tolerance: float, ground_truths: pd.DataFrame, detections: pd.DataFrame\n) -> pd.DataFrame:\n    \"\"\"Match detections to ground truth events. Arguments are taken from a common event x tolerance x series_id evaluation group.\"\"\"\n    detections_sorted = detections.sort_values(score_column_name, ascending=False).dropna()\n    is_matched = np.full_like(detections_sorted[event_column_name], False, dtype=bool)\n    gts_matched = set()\n    for i, det in enumerate(detections_sorted.itertuples(index=False)):\n        best_error = tolerance\n        best_gt = None\n\n        for gt in ground_truths.itertuples(index=False):\n            error = abs(getattr(det, time_column_name) - getattr(gt, time_column_name))\n            if error < best_error and gt not in gts_matched:\n                best_gt = gt\n                best_error = error\n\n        if best_gt is not None:\n            is_matched[i] = True\n            gts_matched.add(best_gt)\n\n    detections_sorted['matched'] = is_matched\n\n    return detections_sorted\n\n\ndef precision_recall_curve(\n        matches: np.ndarray, scores: np.ndarray, p: int\n) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:\n    if len(matches) == 0:\n        return [1], [0], []\n\n    # Sort matches by decreasing confidence\n    idxs = np.argsort(scores, kind='stable')[::-1]\n    scores = scores[idxs]\n    matches = matches[idxs]\n\n    distinct_value_indices = np.where(np.diff(scores))[0]\n    threshold_idxs = np.r_[distinct_value_indices, matches.size - 1]\n    thresholds = scores[threshold_idxs]\n\n    # Matches become TPs and non-matches FPs as confidence threshold decreases\n    tps = np.cumsum(matches)[threshold_idxs]\n    fps = np.cumsum(~matches)[threshold_idxs]\n\n    precision = tps / (tps + fps)\n    precision[np.isnan(precision)] = 0\n    recall = tps / p  # total number of ground truths might be different than total number of matches\n\n    # Stop when full recall attained and reverse the outputs so recall is non-increasing.\n    last_ind = tps.searchsorted(tps[-1])\n    sl = slice(last_ind, None, -1)\n\n    # Final precision is 1 and final recall is 0\n    return np.r_[precision[sl], 1], np.r_[recall[sl], 0], thresholds[sl]\n\n\ndef average_precision_score(matches: np.ndarray, scores: np.ndarray, p: int) -> float:\n    precision, recall, _ = precision_recall_curve(matches, scores, p)\n    # Compute step integral\n    return -np.sum(np.diff(recall) * np.array(precision)[:-1])\n\n\ndef event_detection_ap(\n        solution: pd.DataFrame,\n        submission: pd.DataFrame,\n        tolerances: Dict[str, List[float]],\n) -> float:\n\n    # Ensure solution and submission are sorted properly\n    solution = solution.sort_values([series_id_column_name, time_column_name])\n    submission = submission.sort_values([series_id_column_name, time_column_name])\n\n    # Extract scoring intervals.\n    if use_scoring_intervals:\n        intervals = (\n            solution\n            .query(\"event in ['start', 'end']\")\n            .assign(interval=lambda x: x.groupby([series_id_column_name, event_column_name]).cumcount())\n            .pivot(\n                index='interval',\n                columns=[series_id_column_name, event_column_name],\n                values=time_column_name,\n            )\n            .stack(series_id_column_name)\n            .swaplevel()\n            .sort_index()\n            .loc[:, ['start', 'end']]\n            .apply(lambda x: pd.Interval(*x, closed='both'), axis=1)\n        )\n\n    # Extract ground-truth events.\n    ground_truths = (\n        solution\n        .query(\"event not in ['start', 'end']\")\n        .reset_index(drop=True)\n    )\n\n    # Map each event class to its prevalence (needed for recall calculation)\n    class_counts = ground_truths.value_counts(event_column_name).to_dict()\n\n    # Create table for detections with a column indicating a match to a ground-truth event\n    detections = submission.assign(matched = False)\n\n    # Remove detections outside of scoring intervals\n    if use_scoring_intervals:\n        detections_filtered = []\n        for (det_group, dets), (int_group, ints) in zip(\n            detections.groupby(series_id_column_name), intervals.groupby(series_id_column_name)\n        ):\n            assert det_group == int_group\n            detections_filtered.append(filter_detections(dets, ints))\n        detections_filtered = pd.concat(detections_filtered, ignore_index=True)\n    else:\n        detections_filtered = detections\n\n    # Create table of event-class x tolerance x series_id values\n    aggregation_keys = pd.DataFrame(\n        [(ev, tol, vid)\n         for ev in tolerances.keys()\n         for tol in tolerances[ev]\n         for vid in ground_truths[series_id_column_name].unique()],\n        columns=[event_column_name, 'tolerance', series_id_column_name],\n    )\n\n    # Create match evaluation groups: event-class x tolerance x series_id\n    detections_grouped = (\n        aggregation_keys\n        .merge(detections_filtered, on=[event_column_name, series_id_column_name], how='left')\n        .groupby([event_column_name, 'tolerance', series_id_column_name])\n    )\n    ground_truths_grouped = (\n        aggregation_keys\n        .merge(ground_truths, on=[event_column_name, series_id_column_name], how='left')\n        .groupby([event_column_name, 'tolerance', series_id_column_name])\n    )\n    # Match detections to ground truth events by evaluation group\n    detections_matched = []\n    for key in aggregation_keys.itertuples(index=False):\n        dets = detections_grouped.get_group(key)\n        gts = ground_truths_grouped.get_group(key)\n        detections_matched.append(\n            match_detections(dets['tolerance'].iloc[0], gts, dets)\n        )\n    detections_matched = pd.concat(detections_matched)\n\n    # Compute AP per event x tolerance group\n    event_classes = ground_truths[event_column_name].unique()\n    ap_table = (\n        detections_matched\n        .query(\"event in @event_classes\")\n        .groupby([event_column_name, 'tolerance']).apply(\n            lambda group: average_precision_score(\n                group['matched'].to_numpy(),\n                group[score_column_name].to_numpy(),\n                class_counts[group[event_column_name].iat[0]],\n            )\n        )\n    )\n    # Average over tolerances, then over event classes\n    mean_ap = ap_table.groupby(event_column_name).mean().sum() / len(event_classes)\n\n    return mean_ap","metadata":{"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2023-09-14T18:37:29.330888Z","iopub.execute_input":"2023-09-14T18:37:29.331842Z","iopub.status.idle":"2023-09-14T18:37:29.414371Z","shell.execute_reply.started":"2023-09-14T18:37:29.331790Z","shell.execute_reply":"2023-09-14T18:37:29.412959Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Event Detection Average Precision (EDAP) is a metric used to evaluate the performance of event detection algorithms in time series and video data. It is particularly useful for assessing how well an algorithm can identify specific events at particular time points. The concept is somewhat similar to metrics used in object detection, like Intersection Over Union (IOU), but adapted for time-based data.\n### **Key Components:**\n* Event: An occurrence you want to detect, like a specific action in a video or a particular pattern in a time series.\n* Ground Truth: The actual occurrences of the events, usually annotated manually or obtained through some other reliable method.\n* Prediction: The events that your algorithm detects, along with a confidence score for each detection.\n* Tolerance: A time window within which a predicted event can be considered a match to a ground truth event.\n\nHow EDAP Works:\n1. Selection (Optional) If you choose to use scoring intervals, only predictions within these intervals are considered. This can help to focus the evaluation on specific time ranges and ignore ambiguous or edge-case events.\n2. Assignment Predicted events are matched to ground-truth events based on their time of occurrence. A predicted event is considered a match if it occurs within a specified time tolerance of a ground-truth event. The highest-confidence unmatched prediction within the allowed tolerance is matched to each ground-truth event.\n3. Scoring For each group of matched and unmatched predictions, an Average Precision (AP) score is calculated. This involves generating a precision-recall curve and calculating the area under this curve (AUC). In this context: Precision is the ratio of true positive predictions to the total number of positive predictions (true positives + false positives). Recall is the ratio of true positive predictions to the total number of actual positive events (true positives + false negatives).\n4. Reduction The final EDAP score is an average of the individual AP scores, first averaged over different tolerances and then over different event classes.\n\n### Why Use EDAP?\nFlexibility: It allows for a range of tolerances, making it adaptable to different kinds of time series and video data.\nFairness: By using scoring intervals, you can make the evaluation more fair by ignoring ambiguous or edge-case events.\nComprehensiveness: It takes into account both false positives and false negatives, providing a balanced view of an algorithm's performance.\nIn summary, EDAP provides a robust and flexible way to evaluate the performance of event detection algorithms in time-based data.","metadata":{}},{"cell_type":"code","source":"train = pd.read_parquet(\"/kaggle/input/zzzs-lightweight-training-dataset-target/Zzzs_train.parquet\")","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:37:29.417514Z","iopub.execute_input":"2023-09-14T18:37:29.417998Z","iopub.status.idle":"2023-09-14T18:37:41.527994Z","shell.execute_reply.started":"2023-09-14T18:37:29.417960Z","shell.execute_reply":"2023-09-14T18:37:41.526839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train_series = pd.read_parquet('/kaggle/input/child-mind-institute-detect-sleep-states/train_series.parquet')\ntest = pd.read_parquet('/kaggle/input/child-mind-institute-detect-sleep-states/test_series.parquet')\nsample_submission = pd.read_csv('/kaggle/input/child-mind-institute-detect-sleep-states/sample_submission.csv')\ntrain_events = pd.read_csv('/kaggle/input/child-mind-institute-detect-sleep-states/train_events.csv')","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:37:41.529888Z","iopub.execute_input":"2023-09-14T18:37:41.531116Z","iopub.status.idle":"2023-09-14T18:37:41.612420Z","shell.execute_reply.started":"2023-09-14T18:37:41.531062Z","shell.execute_reply":"2023-09-14T18:37:41.611311Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def create_features(df):\n    # parse the timestamp and create an \"hour\" feature\n    df[\"timestamp\"] = pd.to_datetime(df[\"timestamp\"],utc=True)\n    df[\"hour\"] = (df[\"timestamp\"].dt.hour).astype('int8')\n    df['minute'] = df['timestamp'].dt.minute\n    df['dayofyear'] = df['timestamp'].dt.dayofyear\n    df['year'] = df['timestamp'].dt.year\n\n    # Calculate the half-hour periods\n    df['half_hour'] = (df['hour'] * 2 + (df['minute'] // 30)).astype('int8')\n    \n    df.drop(columns=['minute'], inplace=True)\n    \n    # feature cross\n    df[\"anglez_times_enmo\"] = abs(df[\"anglez\"]) * df[\"enmo\"].astype('float32')\n    # \"rolling\" features\n    periods = 50\n    df[\"anglez_diff\"] = df.groupby('series_id')['anglez'].diff(periods=periods).fillna(method=\"bfill\").astype('float32')\n    df[\"enmo_diff\"]   = df.groupby('series_id')['enmo'].diff(periods=periods).fillna(method=\"bfill\").astype('float32')\n    df[\"anglez_rolling\"] = df[\"anglez\"].rolling(periods,center=True).mean().fillna(method=\"bfill\").fillna(method=\"ffill\").astype('float32')\n    df[\"enmo_rolling\"]   = df[\"enmo\"].rolling(periods,center=True).mean().fillna(method=\"bfill\").fillna(method=\"ffill\").astype('float32')\n    df[\"anglez_diff_rolling\"] = df[\"anglez_diff\"].rolling(periods,center=True).mean().fillna(method=\"bfill\").fillna(method=\"ffill\").astype('float32')\n    df[\"enmo_diff_rolling\"]   = df[\"enmo_diff\"].rolling(periods,center=True).mean().fillna(method=\"bfill\").fillna(method=\"ffill\").astype('float32')\n    \n    return df\n\nfeatures = [\"hour\",\"anglez_times_enmo\", \"half_hour\",\n           \"anglez\",\"anglez_diff\",\"anglez_rolling\",\"anglez_diff_rolling\",\n           \"enmo\",\"enmo_diff\",\"enmo_rolling\",\"enmo_diff_rolling\"]\ntrain = create_features(train)\ntest = create_features(test)","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:37:41.614210Z","iopub.execute_input":"2023-09-14T18:37:41.614780Z","iopub.status.idle":"2023-09-14T18:40:01.525348Z","shell.execute_reply.started":"2023-09-14T18:37:41.614745Z","shell.execute_reply":"2023-09-14T18:40:01.523755Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"series_id = train['series_id'].unique()","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:40:01.528918Z","iopub.execute_input":"2023-09-14T18:40:01.529359Z","iopub.status.idle":"2023-09-14T18:40:02.767455Z","shell.execute_reply.started":"2023-09-14T18:40:01.529323Z","shell.execute_reply":"2023-09-14T18:40:02.766193Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nfor id in series_id[0:3]:\n\n    temp = train[train['series_id'] == id]\n    # vizualize these series\n    display(Markdown('###  anglez for series ' + id))\n    fig, ax = plt.subplots(figsize=(20, 3))\n    sns.lineplot(data=temp, x=\"step\", y=\"anglez\",hue=\"awake\", linewidth = 0.5)\n    plt.show();\n    display(Markdown('###  enmo for series ' + id))\n    fig, ax = plt.subplots(figsize=(20, 3))\n    sns.lineplot(data=temp, x=\"step\", y=\"enmo\",hue=\"awake\", linewidth = 0.5)\n    plt.show();\n    del temp\n    gc.collect();","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:40:02.769339Z","iopub.execute_input":"2023-09-14T18:40:02.769961Z","iopub.status.idle":"2023-09-14T18:40:35.985500Z","shell.execute_reply.started":"2023-09-14T18:40:02.769924Z","shell.execute_reply":"2023-09-14T18:40:35.983900Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"series = train[train['series_id']== np.random.choice(series_id)]\nseries = series[series['dayofyear'] == series.dayofyear.unique()[1]]","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:56:24.789891Z","iopub.execute_input":"2023-09-14T18:56:24.790584Z","iopub.status.idle":"2023-09-14T18:56:27.264754Z","shell.execute_reply.started":"2023-09-14T18:56:24.790547Z","shell.execute_reply":"2023-09-14T18:56:27.263510Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"series['anglez_int'] = (series['anglez'] * 100).astype(int)\nseries['enmo_int'] = (series['enmo'] * 100).astype(int)\nseries.sample()","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:56:27.840526Z","iopub.execute_input":"2023-09-14T18:56:27.841055Z","iopub.status.idle":"2023-09-14T18:56:27.873157Z","shell.execute_reply.started":"2023-09-14T18:56:27.841018Z","shell.execute_reply":"2023-09-14T18:56:27.871973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.barplot(data=series, x=\"hour\", y=\"anglez_int\", hue='awake', linewidth = 0.5)","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:56:32.969792Z","iopub.execute_input":"2023-09-14T18:56:32.970270Z","iopub.status.idle":"2023-09-14T18:56:34.873607Z","shell.execute_reply.started":"2023-09-14T18:56:32.970228Z","shell.execute_reply":"2023-09-14T18:56:34.872280Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(12, 6)) \nsns.barplot(data=series, x=\"half_hour\", y=\"enmo_int\", hue='awake', linewidth = 0.5, )","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:56:56.842514Z","iopub.execute_input":"2023-09-14T18:56:56.844041Z","iopub.status.idle":"2023-09-14T18:57:00.325290Z","shell.execute_reply.started":"2023-09-14T18:56:56.843985Z","shell.execute_reply":"2023-09-14T18:57:00.323982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"temp = series.groupby('half_hour')['enmo_int'].sum().reset_index()","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:57:17.700532Z","iopub.execute_input":"2023-09-14T18:57:17.700951Z","iopub.status.idle":"2023-09-14T18:57:17.711981Z","shell.execute_reply.started":"2023-09-14T18:57:17.700917Z","shell.execute_reply":"2023-09-14T18:57:17.710572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/blob/master/Chapter1_Introduction/Ch1_Introduction_TFP.ipynb\nimport tensorflow as tf\n\nimport tensorflow_probability as tfp\ntfd = tfp.distributions\ntfb = tfp.bijectors\n","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:57:19.678861Z","iopub.execute_input":"2023-09-14T18:57:19.679304Z","iopub.status.idle":"2023-09-14T18:57:19.687318Z","shell.execute_reply.started":"2023-09-14T18:57:19.679244Z","shell.execute_reply":"2023-09-14T18:57:19.685890Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def joint_log_prob(count_data, lambda_1, lambda_2, tau):\n    tfd = tfp.distributions\n \n    alpha = (1. / tf.reduce_mean(count_data))\n    rv_lambda_1 = tfd.Exponential(rate=alpha)\n    rv_lambda_2 = tfd.Exponential(rate=alpha)\n \n    rv_tau = tfd.Uniform()\n \n    lambda_ = tf.gather(\n         [lambda_1, lambda_2],\n         indices=tf.cast(tau * tf.cast(tf.size(count_data), dtype=tf.float32) <= tf.cast(tf.range(tf.size(count_data)), dtype=tf.float32), dtype=tf.int32))\n    rv_observation = tfd.Poisson(rate=lambda_)\n \n    return (\n         rv_lambda_1.log_prob(lambda_1)\n         + rv_lambda_2.log_prob(lambda_2)\n         + rv_tau.log_prob(tau)\n         + tf.reduce_sum(rv_observation.log_prob(count_data))\n    )\n\n\n# Define a closure over our joint_log_prob.\ndef unnormalized_log_posterior(lambda1, lambda2, tau):\n    return joint_log_prob(count_data, lambda1, lambda2, tau)","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:57:21.491241Z","iopub.execute_input":"2023-09-14T18:57:21.491780Z","iopub.status.idle":"2023-09-14T18:57:21.503727Z","shell.execute_reply.started":"2023-09-14T18:57:21.491741Z","shell.execute_reply":"2023-09-14T18:57:21.502235Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Defining our Data and assumptions\ncount_data = tf.constant(temp['enmo_int'], dtype=tf.float32)\nn_count_data = tf.shape(count_data)\ndays = tf.range(n_count_data[0],dtype=tf.int32)\n\n# Visualizing the Results\nplt.figure(figsize=(12.5, 4))\nplt.bar(days.numpy(), count_data, color=\"#5DA5DA\")\nplt.xlabel(\"Time (days)\")\nplt.ylabel(\"count of text-msgs received\")\nplt.title(\"Did the user's texting habits change over time?\")\nplt.xlim(0, n_count_data[0].numpy());","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:57:22.214985Z","iopub.execute_input":"2023-09-14T18:57:22.215452Z","iopub.status.idle":"2023-09-14T18:57:22.715337Z","shell.execute_reply.started":"2023-09-14T18:57:22.215415Z","shell.execute_reply":"2023-09-14T18:57:22.713832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# wrap the mcmc sampling call in a @tf.function to speed it up\n@tf.function(autograph=False)\ndef graph_sample_chain(*args, **kwargs):\n  return tfp.mcmc.sample_chain(*args, **kwargs)\n\nnum_burnin_steps = 5000\nnum_results = 20000\n\n\n# Set the chain's start state.\ninitial_chain_state = [\n    tf.cast(tf.reduce_mean(count_data), tf.float32) * tf.ones([], dtype=tf.float32, name=\"init_lambda1\"),\n    tf.cast(tf.reduce_mean(count_data), tf.float32) * tf.ones([], dtype=tf.float32, name=\"init_lambda2\"),\n    0.5 * tf.ones([], dtype=tf.float32, name=\"init_tau\"),\n]\n\n\n# Since HMC operates over unconstrained space, we need to transform the\n# samples so they live in real-space.\nunconstraining_bijectors = [\n    tfp.bijectors.Exp(),       # Maps a positive real to R.\n    tfp.bijectors.Exp(),       # Maps a positive real to R.\n    tfp.bijectors.Sigmoid(),   # Maps [0,1] to R.  \n]\n\nstep_size = 0.2\n\nkernel=tfp.mcmc.TransformedTransitionKernel(\n        inner_kernel=tfp.mcmc.HamiltonianMonteCarlo(\n            target_log_prob_fn=unnormalized_log_posterior,\n            num_leapfrog_steps=2,\n            step_size=step_size,\n            state_gradients_are_stopped=True),\n        bijector=unconstraining_bijectors)\n\nkernel = tfp.mcmc.SimpleStepSizeAdaptation(\n    inner_kernel=kernel, num_adaptation_steps=int(num_burnin_steps * 0.8))\n\n\n# Sample from the chain.\n[\n    lambda_1_samples,\n    lambda_2_samples,\n    posterior_tau,\n], kernel_results = graph_sample_chain(\n    num_results=num_results,\n    num_burnin_steps=num_burnin_steps,\n    current_state=initial_chain_state,\n    kernel = kernel)\n    \ntau_samples = tf.floor(posterior_tau * tf.cast(tf.size(count_data),dtype=tf.float32))","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:57:28.074935Z","iopub.execute_input":"2023-09-14T18:57:28.075409Z","iopub.status.idle":"2023-09-14T18:57:59.148523Z","shell.execute_reply.started":"2023-09-14T18:57:28.075355Z","shell.execute_reply":"2023-09-14T18:57:59.147352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"acceptance rate: {}\".format(\n    tf.reduce_mean(tf.cast(kernel_results.inner_results.inner_results.is_accepted,dtype=tf.float32))))\nprint(\"final step size: {}\".format(\n    tf.reduce_mean(kernel_results.inner_results.inner_results.accepted_results.step_size[-100:])))\n","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:57:59.150651Z","iopub.execute_input":"2023-09-14T18:57:59.151551Z","iopub.status.idle":"2023-09-14T18:57:59.161727Z","shell.execute_reply.started":"2023-09-14T18:57:59.151509Z","shell.execute_reply":"2023-09-14T18:57:59.160400Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.hist(tau_samples, bins=10)","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:58:06.814473Z","iopub.execute_input":"2023-09-14T18:58:06.815427Z","iopub.status.idle":"2023-09-14T18:58:14.023357Z","shell.execute_reply.started":"2023-09-14T18:58:06.815355Z","shell.execute_reply":"2023-09-14T18:58:14.022464Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# tau_samples, lambda_1_samples, lambda_2_samples contain\n# N samples from the corresponding posterior distribution\n\nN_ = tau_samples.shape[0]\nexpected_texts_per_day = tf.zeros(N_,n_count_data.shape[0]) #(10000,74)\n\nplt.figure(figsize=(12.5, 9))\n\nday_range = tf.range(0,n_count_data[0],delta=1,dtype = tf.int32)\n\n# expand from shape of 74 to (10000,74)\nday_range = tf.expand_dims(day_range,0)\nday_range = tf.tile(day_range,tf.constant([N_,1]))\n\n# expand from shape of 10000 to 10000,74\ntau_samples_per_day = tf.expand_dims(tau_samples,0)\ntau_samples_per_day = tf.transpose(tf.tile(tau_samples_per_day,tf.constant([day_range.shape[1],1])))\n\ntau_samples_per_day = tf.cast(tau_samples_per_day,dtype=tf.int32)\n#ix_day is (10000,74) tensor where axis=0 is number of samples, axis=1 is day. each value is true iff sampleXday value is < tau_sample value\nix_day = day_range < tau_samples_per_day\n\nlambda_1_samples_per_day = tf.expand_dims(lambda_1_samples,0)\nlambda_1_samples_per_day = tf.transpose(tf.tile(lambda_1_samples_per_day,tf.constant([day_range.shape[1],1])))\nlambda_2_samples_per_day = tf.expand_dims(lambda_2_samples,0)\nlambda_2_samples_per_day = tf.transpose(tf.tile(lambda_2_samples_per_day,tf.constant([day_range.shape[1],1])))\n\nexpected_texts_per_day = ((tf.reduce_sum(lambda_1_samples_per_day*tf.cast(ix_day,dtype=tf.float32),axis=0) + tf.reduce_sum(lambda_2_samples_per_day*tf.cast(~ix_day,dtype=tf.float32),axis=0))/N_)\n\nplt.plot(range(n_count_data[0]), expected_texts_per_day, lw=4, color=\"#E24A33\",\n         label=\"expected number of text-messages received\")\nplt.xlim(0, n_count_data.numpy()[0])\nplt.xlabel(\"Day\")\nplt.ylabel(\"Expected # text-messages\")\nplt.title(\"Expected number of text-messages received\")\nplt.ylim(0, 6000)\nplt.bar(np.arange(len(count_data)), count_data, color=\"#5DA5DA\", alpha=0.65,\n        label=\"observed texts per day\")\n\nplt.legend(loc=\"upper left\");","metadata":{"execution":{"iopub.status.busy":"2023-09-14T18:58:18.618981Z","iopub.execute_input":"2023-09-14T18:58:18.619417Z","iopub.status.idle":"2023-09-14T18:58:19.177573Z","shell.execute_reply.started":"2023-09-14T18:58:18.619370Z","shell.execute_reply":"2023-09-14T18:58:19.176298Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}],"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}}