predict: return the unsurprising portion of a metric stream by estimating trend and seasonal effects

predict.juttle

// Predict:
// Trend, seasonality, and noise estimation for metric streams.
//
// Predict is a sub which consumes a metric stream and emits a
// prediction stream based on estimated trend, seasonality, and a
// smoothed version of the input. Prediction errors (the difference
// between predict's output and the input) can signal anomalies in the
// stream.
//
// usage: predict  [options]
// options:
// -field ("value") Name of metric field to predict
// -over (:w:)  Period: the length of one repeating cycle in the input metric.
// -every (default is based on period) interval between emitted prediction points
// -pct (0.5) percentile to retain during initial reduction to every-spaced values
// -nonneg (true) do not allow predictions to become negative.
// -detrend (true) if false, do not remove estimated trend
// -deseason (true)  if false, do not remove estimated cyclic effect
// -denoise (true)  if false, do not smooth the detrended/deseasoned value
//
// predict is intended to be used with a historic or superquery that
// includes enough history to initialize its estimators.  predict can
// begin emitting prediction points almost immediately, though the
// early points will simply be de-meaned. After 2 periods have gone by,
// trend and seasonality (at the resolution of -every) will switch on,
// and after 3 periods have passed all estimators have full windows of
// data.
//
// output:
// each output point contains the following fields:
// *field : average value of input field over -every
// T : portion of field predicted by the trend estimator
// S : portion of field predicted by seasonality
// Y : portion of field predicted by smoothing
// P : predicted value of field, T + S + Y
// E : prediction error, P - field
// Z : normalized error based on the sample stddev over trailing periods
//


// Return the sample Z-score of the specified field.
//
// Parameters:
//
//   * field:   fieldname
//
// compute an approximate Z-score for an input value, using sample
// estimates of population mean and stddev:
//     (value - sample mean) / sample stddev
//
reducer z(field) {
    var sum = 0;
    var ssum = 0;
    var n = 0;
    var latest = null;
    function update() {
        var v = *field;
        if (v != null) {
            latest = v;
            sum = sum + v;
            ssum = ssum + v*v;
            n = n + 1;
        }
    }
    function result() {
        if (n < 2 ) {
            return null;
        } else {
            var mean = sum / n;
            var stddev = Math.sqrt(1/(n-1) * (ssum - sum*sum/n));
            return (latest - mean) / stddev;
        }
    }
    function expire() {
        var v = *field;
        if (v != null) {
            sum = sum - v;
            ssum = ssum - v*v;
            n = n - 1;
        }
    }
};

// Exponential weighted moving average
//
// Parameters:
//
//   * field:   fieldname
//   * alpha:   weight of newest value in moving average
//   * initial: initial value of the moving average, defaults to first update
//   * offset:  name of field containing an offset to be applied to the stored value
//
reducer ewma(field, alpha, initial=null, offset=null) {
    var value = initial;

    function update() {
        if (*field != null) {
            if (value == null) {
                value = *field;
            }
            value = value + alpha * (*field - value);
            if (offset != null && *offset != null) {
                value = value - *offset;
            }
        }
    }
    function result() {
        return value;
    }
}

// Return the Nth lagged value of a field
//
// Parameters:
//
//   * field: name of value field
//   * N: (1) lag
//
reducer L(field, N=1, undef=null) {
    var values = [];
    var n = 0, i = 0;
    var current = null;

    function update() {
        if (*field != null) {
            values[i] = *field;
            i = (i + 1) % (N + 1);
            if (n <= N) {
                n = n + 1;
            }
        }
    }
    function result() {
        if (n == N + 1) {
            return values[i];
        } else {
            return undef;
        }
    }
}

// Return a nice subdividing duration, preserving calendarness where possible
//
// Parameters:
//
//   * d: duration to subdivide
//
export function subdivide_duration(d) {
    if (Duration.as(d, "y") == Math.floor(Duration.as(d, "y"))) {
        return :M:;
    } else if (Duration.as(d, "M") == Math.floor(Duration.as(d, "M"))) {
        return :d:;
    } else if (Duration.as(d, "w") == Math.floor(Duration.as(d, "w"))) {
        return :6h:;
    } else if (Duration.as(d, "d") == Math.floor(Duration.as(d, "d"))) {
        return :h:;
    } else if (Duration.as(d, "h") == Math.floor(Duration.as(d, "h"))) {
        return :5m:;
    } else {
        return d / 30;
    }
}

// estimate trailing trend and median as the median
// duration-over-duration change of all samples in a window of -over
// duration (a variant of the Theil-Sen estimator). This can use up to
// 2 periods of historic data per point, but begins outputting a
// windowed median immediately.  trend expects its input to already be
// aggregated -every.
//
// Parameters:
//
//   * field: Name of metric field to predict
//   * over:  Period: length of one repeating cycle in the input metric.
//   * every: interval between emitted prediction points
//   * revise: if true, include the current value in the slope estimate.
//
// Output:
//
//   * T: portion of field predicted by trend
//   * M: slope (change per over), or null if no slope was estimated
//        (T will have a trailing median value if M is null)
//
sub trend(field, over, every, revise=true) {
    const N = Math.round(over / every);
    put __bucket = count() % N
    | put __per_over = delta(field,null) by __bucket // same bucket, successive periods
    | put -over 2 * over
        __t0 = first(time) + (last(time) - first(time)) / 2,
        __y0 = percentile(field,0.5),
        M = __per_over != null ? percentile(__per_over,0.5) : null
    | put M = revise ? M : L("M") ?? M // use current or previous M/y0/t0 to predict trend at time
    | put __y0 = revise ? __y0 : L("__y0") ?? __y0
    | put __t0 = revise ? __t0 : L("__t0") ?? __t0
    | put T = (M != null) ? __y0 + M * (time - __t0) / over : __y0
    | remove __t0, __y0, __bucket, __per_over
}

// estimate cyclic variation as a moving median of the values over
// successive epochs at a given phase.  This can use up to 3 periods
// of historic data per point, and at startup simply replays the most
// recent value.  seasonal expects its input to already be de-trended
// and aggregated -every.
//
// Parameters:
//
//   * field: Name of metric field to predict
//   * over:  Period: length of one repeating cycle in the input metric.
//   * every: interval between emitted prediction points
//   * revise: if true, include the current value in the seasonal estimate.
//
// Output:
//
//   * S: portion of field predicted by seasonality
//
sub seasonal(field, over, every, revise) {
    const N = Math.round(over / every);
    put __bucket = count() % N, __n=count()
    | put -over 3 * over S = (__n > 4) ? percentile(field,0.5) : null by __bucket
    | put S = revise ? S : L("S") by __bucket // use current or previous seasonal at time
    | remove __bucket, __n
}

// windowed median as a less-noisy prediction of mean. the optional
// offset field is examined for a nonzero impulse to be added to the
// windowed values state. This is allows incorporation of a
// discontinuity in the series.
//
// Parameters:
//
//   * field: Name of metric field to predict
//   * every: interval between emitted prediction points
//   * offset: fieldname for an impulse to be added
//   * revise: if false, return previous value as the prediction
//
// Output:
//
//   * Y : portion of field predicted by smoothing
//
sub level(field, every, offset, revise) {
    put -over 3 * every Y = percentile(field, 0.5)
    | put Y = revise ? Y : (L("Y") ?? Y)
    | put -over 3 * every Y = Y + (sum(offset) ?? 0)
}

// predict the mean of the distribution of next field value at the
// current time from prior values (revise=false), optionally including
// the current field value in the trend and seasonal estimates
// (revise=true).  This can use up to 3 periods of historic data per
// point, but can begin producing (noisy) point-to-point estimates
// after a few points.
//
// Parameters:
//
//   * field: ("value") Name of metric field to predict
//   * over: (:w:)  Period: length of one repeating cycle in the input metric.
//   * every: (default is based on value of over) interval between emitted prediction points
//   * pct: (0.5) percentile to retain during initial reduction to every-spaced values
//   * nonneg: (true) do not allow predictions to become negative.
//   * detrend: (true) if false, do not remove estimated trend
//   * deseason: (true)  if false, do not remove estimated cyclic effect
//   * denoise: (true)  if false, do not smooth the detrended/deseasoned value
//   * revise: (true) if false, do not include current value in trend and seasonal estimates
//
// Output:
//
// predict consumes its input stream, and outputs points every -every.
// each output point contains the following fields:
//   * *field: percentile value of input field over -every, from -pct (default median)
//   * T: portion of field predicted by the trend estimator
//   * S: portion of field predicted by seasonality
//   * Y: portion of field predicted by smoothing
//   * P: predicted value of field, T + S + Y
//   * E: prediction error, P - field
//   * Z: normalized error based on the sample stddev over trailing periods
//
export sub predict(field = "value", over = :w:, every=null, pct=0.5, nonneg=true,
                   detrend=true, deseason=true, denoise=true, revise=true) {
    const __every = every ?? subdivide_duration(over);
    reduce -every __every time= Date.quantize(first(time), __every), *field = percentile(field, pct)
    | trend -field field -every __every -over over -revise revise
    | put T = detrend ? T : 0
    | put __detrend = *field  - T
    | put __LT = L("T") ?? 0 // before M becomes valid, T is windowed median.
    | put __offset = ((M != null && L("M") == null) ?  __LT - T : 0)
    | seasonal -field '__detrend' -over over -every __every -revise revise
    | put S = deseason ? S : 0
    | put __offset = __offset - ((S != null && L("S") == null) ? S : 0)
    | put __deseason  = __detrend  - (S ?? 0)
    | level -field '__deseason' -every __every -revise revise -offset '__offset'
    | put Y = denoise ? Y : (L("__deseason") ?? 0) + __offset
    | put P = (S ?? 0) + (T ?? 0) + Y
    | put P = nonneg ? Math.max(P, 0) : P
    | put E = P - *field
    | put -over 3*over Z = z('E')
    | remove __LT, __detrend, __deseason, __offset, M
}

//////////////////////////////////////////////////////////////////////////////

// predict demos: trend and seasonality
//
// the following metric sources all work with the default predict interval of :w:
// all reads are from jx3 or j4qa as indicated.
//

input last: duration -default :30 days: -label 'Show this duration:';

export sub live_juttle_runs(dur1=:3M:) {
    read -from (:now: + :1d: - dur1) -space 'prod'
    event = 'run-juttle' AND
    properties_page_currentPage ='https://app.jut.io/#explorer' AND
    context_ip != '207.141.12.50'
    | reduce -every :h: run_count = count()
    | (@timechart -title "Weekly App Juttle Program Runs" -valueField 'run_count' -display.dataDensity 0; put value=run_count)
}

export sub live_pageviews() {
    read -from :now: - last -space 'prod' type = 'page' properties_url ~ 'http://docs.jut.io/*'
        NOT (context_ip in ['207.141.12.50','54.149.158.203','54.68.207.93','52.10.135.75','54.148.36.156','52.24.101.163'])
    | reduce -every :h: pageviews = count(), users = count_unique(anonymousId)
    | (@timechart -title "live pageviews" -valueField 'pageviews' -display.dataDensity 0; put value=pageviews)
}

export sub live_users() {
    read -from :now: - last -space 'prod' type = 'page' properties_url ~ 'http://docs.jut.io/*'
        NOT (context_ip in ['207.141.12.50','54.149.158.203','54.68.207.93','52.10.135.75','54.148.36.156','52.24.101.163'])
    | reduce -every :h: pageviews = count(), users = count_unique(anonymousId)
    | (@timechart -title "live users" -valueField 'users' -display.dataDensity 0; put value=users)
}

export sub slack_messages() {
    // this runs in j4qa rather than jx3
    read -from :-M: -space "slack"
        type = "message" and user_name != null
    | reduce -every :h: name="slack_messages", value = count()
    | (@timechart -title "slack messages" -display.dataDensity 0; pass)
}

sub pulse(field, value, start, dur) {
  put *field = *field + ((time >= start && time < start + dur) ? value : 0)
}

sub slide(field, value, start, dur) {
    put *field = *field + (
        (time >= start && time < start + dur)
            ? ((time - start) / dur) * value
            : 0)
}

export sub series(from, to, every, over, trend = 1, season = 10, alpha = 0.5, sigma=0.1, initial=0) {
    // generate a series of values having given trend, seasonality, noise, and autocorrelation (ewma alpha)
    emit -from from -to to -every every
    | put n = count(), dy = (time - Date.new(0)) / over, cycle = Math.sin(dy * 2 * Math.PI)
    | put step = sigma * (2 * Math.random() - 1)
    | put step = ewma("step", alpha, initial)
    | put value = n * trend + season * cycle + step
    | slide -field 'value' -value -20 -start from + 5.25 * over -dur 4*every
}

export sub posSeries(from, to, every, over, trend = 1, season = 10, alpha = 0.5, sigma=0.1, initial=0) {
    // like series, but keep it positive
    emit -from from -to to -every every
    | put n = count(), dy = (time - Date.new(0)) / over, cycle = Math.abs(Math.sin(dy * 2 * Math.PI))
    | put step = sigma * (2 * Math.random() - 1)
    | put step = ewma("step", alpha, initial)
    | put value = n * trend + season * cycle + step
    | slide -field 'value' -value -20 -start from + 5.25 * over -dur 4*every
    | put value = Math.abs(value)
}

// null invocation of predict that uses the previous value as the
// prediction of the next value. not bad.
//
export sub dontPredict(field = "value", over = :w:, every=null, pct=0.5, nonneg=true) {
    predict -field field -over over -every every -pct pct -nonneg nonneg
        -detrend false -deseason false -denoise false -revise false
}

const over = :w:;
const every = subdivide_duration(over);
// uncomment any one of the following sources:
posSeries -from Date.new(0) -to Date.new(0)+ 8*over -every every -over over -trend .05 -season 20 -alpha 0.5 -sigma 4//live_juttle_runs
//live_pageviews
//live_users
//slack_messages
//  | dontPredict
    | predict
    | (
        put name="stderr", value = sigma("E") | @tile;
        split value, P, T, S, Y   | @timechart -keyField 'name' -display.dataDensity 0;
        split E   | @timechart -keyField 'name' -display.dataDensity 0;
        split Z   | @timechart -keyField 'name' -display.dataDensity 0;
    )