Как вы создаете прогноз линейной регрессии по данным временных рядов в python

Ответ 1

В тексте вашего вопроса вы четко заявляете, что хотите верхняя и нижняя границы вашего регрессионного выхода, а также выход прогнозирование. Вы также упоминаете использование алгоритмов Holt-Winters для прогнозирование в частности.

Пакеты, предлагаемые другими респондентами, полезны, но вы можете заметить что sklearn LinearRegression не дает вам границ ошибок "из box ", statsmodels делает не предоставляет Holt-Winters прямо сейчас.

Поэтому я предлагаю использовать эту реализацию Holt-Winters. К сожалению, ее лицензия неясна, поэтому я не могу воспроизвести ее здесь полный. Теперь я не уверен, действительно ли вы хотите, чтобы Holt-Winters (сезонное) предсказание или линейное экспоненциальное сглаживание Холта алгоритм. Я предполагаю, что последний дал название должности. Таким образом, вы можете использовать функцию linear() связанной библиотеки. техника подробно описана здесь для заинтересованных читателей.

В интересах не предоставления ссылки только ответ - я опишу основные функции здесь. Определяется функция, которая принимает данные i.e.

 def linear(x, fc, alpha = None, beta = None):

x - это данные, которые нужно поместить, fc - количество требуемых временных меток для прогнозирования, альфа и бета принимают свои обычные значения Холт-Винтерса: примерно параметр, чтобы контролировать количество сглаживания на "уровень", и к "тренду" соответственно. Если alpha или beta не они оцениваются с помощью scipy.optimize.fmin_l_bfgs_b для минимизации RMSE.

Функция просто применяет алгоритм Холта путем циклического существующих точек данных, а затем возвращает прогноз следующим образом:

 return Y[-fc:], alpha, beta, rmse

где Y[-fc:] - точки прогноза, alpha и beta - это фактически используемые значения, а rmse - среднеквадратичная ошибка. К сожалению, как вы можете видеть, нет никакой верхней или меньшей уверенности интервалы. Кстати, мы должны, вероятно, называть их предсказанием интервалы.

Интервалы предсказания maths

Алгоритм Холта и алгоритм Холта-Винтера являются экспоненциальным сглаживанием методы и поиск доверительных интервалов для генерируемых прогнозов из них сложный вопрос. Они называются " правило thumb " и, в случае мультипликативной формы Холта-Винтера алгоритм без базовая статистическая модель. Однако окончательная сноска к этой странице утверждает, что:

Можно рассчитать доверительные интервалы вокруг долгосрочных прогнозы, производимые экспоненциальными сглаживающими моделями, путем их учета в качестве особых случаев моделей ARIMA. (Остерегайтесь: не все программное обеспечение вычисляет доверительные интервалы для этих моделей правильно.) Ширина доверительные интервалы зависят от (i) ошибки RMS модели, (ii) тип сглаживания (простой или линейный); (iii) значение (и) константа сглаживания; и (iv) количество периодов, в течение которых вы находитесь прогнозирования. В общем, интервалы распространяются быстрее, когда α получает больше в модели SES, и они распространяются намного быстрее, когда линейные вместо простого сглаживания.

Мы видим здесь, что модель ARIMA (0,2,2) эквивалентна модели Holt линейная модель с аддитивными ошибками

Код интервалов предсказания (т.е. как действовать)

В комментариях вы указываете, что вы "можете легко сделать это в R" . я предположим, что вы можете использовать функцию holt(), предоставляемую forecast в R и, следовательно, ожидая таких интервалов. В в этом случае - вы можете адаптировать библиотеку python, чтобы предоставить их вам на тот же базис.

Глядя на R holt code, мы видим, что он возвращает объект на основе forecast(ets(...). Под капотом - это в конечном итоге вызывает эта функция class1, которая возвращает среднее значение mu и дисперсию var (а также cj, который я должен признаться, я не понимаю). Дисперсия используется для вычисления верхней и нижней границ здесь.

Чтобы сделать что-то подобное в Python - нам нужно было что-то сделать аналогично функции class1 R, которая оценивает дисперсию каждого прогнозирование. Эта функция принимает остатки, найденные в моделировании и умножает их на один раз на каждом временном шаге, чтобы получить дисперсию в это timestep. В частном случае линейного алгоритма Холта фактор является суммой сумм alpha + k*beta где k - это число предсказаний временных меток. Когда у вас есть это дисперсии в каждой точке прогнозирования, обрабатывать ошибки как обычно распределяется и получает значение X% из нормального распределения.

Вот идея, как это сделать в Python (используя код, который я связал как ваша линейная функция)

#Copy, import or reimplement the RMSE and linear function from
#https://gist.github.com/andrequeiroz/5888967

#factor in case there are not 1 timestep per day - in your case
#assuming the timesteps are UTC epoch - I think they're 5 min
# spaced i.e. 288 per day
timesteps_per_day = 288

# Note - big assumption here - your known data will be regular in time
# i.e. timesteps_per_day observations per day.  From the timestamps this seems valid.
# if you can't guarantee that - you'll need to interpolate the data
def holt_predict(data, timestamps, forecast_days, pred_error_level = 0.95):
    forecast_timesteps = forecast_days*timesteps_per_day
    middle_predictions, alpha, beta, rmse = linear(data,int(forecast_timesteps))
    cum_error = [beta+alpha]
    for k in range(1,forecast_timesteps):
        cum_error.append(cum_error[k-1] + k*beta + alpha)

    cum_error = np.array(cum_error)
    #Use some numpy multiplication to get the intervals
    var = cum_error * rmse**2
    # find the correct ppf on the normal distribution (two-sided)
    p = abs(scipy.stats.norm.ppf((1-pred_error_level)/2))
    interval = np.sqrt(var) * p
    upper = middle_predictions + interval
    lower = middle_predictions - interval
    fcast_timestamps = [timestamps[-1] + i * 86400 / timesteps_per_day for i in range(forecast_timesteps)]

    ret_value = []

    ret_value.append({'target':'Forecast','datapoints': zip(middle_predictions, fcast_timestamps)})
    ret_value.append({'target':'Upper','datapoints':zip(upper,fcast_timestamps)})
    ret_value.append({'target':'Lower','datapoints':zip(lower,fcast_timestamps)})
    return ret_value

if __name__=='__main__':
    import numpy as np
    import scipy.stats
    from math import sqrt

    null = None
    data_in = [{"target": "average", "datapoints": [[null, 1435688679],
    [34.870499801635745, 1435688694], [null, 1435688709], [null,
    1435688724], [null, 1435688739], [null, 1435688754], [null, 1435688769],
    [null, 1435688784], [null, 1435688799], [null, 1435688814], [null,
    1435688829], [null, 1435688844], [null, 1435688859], [null, 1435688874],
    [null, 1435688889], [null, 1435688904], [null, 1435688919], [null,
    1435688934], [null, 1435688949], [null, 1435688964], [null, 1435688979],
    [38.180000209808348, 1435688994], [null, 1435689009], [null,
    1435689024], [null, 1435689039], [null, 1435689054], [null, 1435689069],
    [null, 1435689084], [null, 1435689099], [null, 1435689114], [null,
    1435689129], [null, 1435689144], [null, 1435689159], [null, 1435689174],
    [null, 1435689189], [null, 1435689204], [null, 1435689219], [null,
    1435689234], [null, 1435689249], [null, 1435689264], [null, 1435689279],
    [30.79849989414215, 1435689294], [null, 1435689309], [null, 1435689324],
    [null, 1435689339], [null, 1435689354], [null, 1435689369], [null,
    1435689384], [null, 1435689399], [null, 1435689414], [null, 1435689429],
    [null, 1435689444], [null, 1435689459], [null, 1435689474], [null,
    1435689489], [null, 1435689504], [null, 1435689519], [null, 1435689534],
    [null, 1435689549], [null, 1435689564]]}]

    #translate the data.  There may be better ways if you're
    #prepared to use pandas / input data is proper json
    time_series = data_in[0]["datapoints"]

    epoch_in = []
    Y_observed = []

    for (y,x) in time_series:
        if y and x:
            epoch_in.append(x)
            Y_observed.append(y)

    #Pass in the number of days to forecast
    fcast_days = 30
    res = holt_predict(Y_observed,epoch_in,fcast_days)
    data_out = data_in + res
    #data_out now holds the data as you wanted it.

    #Optional plot of results
    import matplotlib.pyplot as plt
    plt.plot(epoch_in,Y_observed)
    m,tstamps = zip(*res[0]['datapoints'])
    u,tstamps = zip(*res[1]['datapoints'])
    l,tstamps = zip(*res[2]['datapoints'])
    plt.plot(tstamps,u, label='upper')
    plt.plot(tstamps,l, label='lower')
    plt.plot(tstamps,m, label='mean')
    plt.show()

N.B. Результат, который я дал, добавляет в ваш объект точки типа tuple. Если вам действительно нужно list, замените zip(upper,fcast_timestamps) на map(list,zip(upper,fcast_timestamps)), где код добавляет upper, lower и forecast dicts к результату.

Этот код предназначен для частного случая линейного алгоритма Холта - это не общий способ вычисления правильных интервалов прогнозирования.

Важное примечание

В ваших примерах входных данных, похоже, много null и только 3 подлинных данных. Это маловероятно, чтобы быть хорошей основой для предсказание таймсерий - особенно, поскольку все они кажутся с 15 минутами, и вы пытаетесь прогнозировать до 3 месяцев!. Действительно - если вы подаете эти данные в R holt(), он скажет:

You've got to be joking. I need more data!

Я предполагаю, что у вас есть больший набор данных для тестирования. Я попробовал код выше на цены открытия фондового рынка на 2015 год и, как представляется, дал разумные результаты (см. Ниже).

Вы можете подумать, что интервалы прогнозирования выглядят немного шире. Этот блог от автора модуля прогноза R подразумевает, что это преднамеренно:

"доверительные интервалы для среднего значения намного более узкие, чем интервалы прогнозирования"

Ответ 2

Scikit - отличный модуль для python

Классы X и Y должны быть разделены на два массива, где, если бы вы их нарисовали (X, Y), индекс из одного совпадал бы с другим массивом.

Итак, я думаю, что в ваших данных временных рядов вы разделили бы значения X и Y следующим образом:

null = None
time_series = [{"target": "average", "datapoints": [[null, 1435688679], [34.870499801635745, 1435688694], [null, 1435688709], [null, 1435688724], [null, 1435688739], [null, 1435688754], [null, 1435688769], [null, 1435688784], [null, 1435688799], [null, 1435688814], [null, 1435688829], [null, 1435688844], [null, 1435688859], [null, 1435688874], [null, 1435688889], [null, 1435688904], [null, 1435688919], [null, 1435688934], [null, 1435688949], [null, 1435688964], [null, 1435688979], [38.180000209808348, 1435688994], [null, 1435689009], [null, 1435689024], [null, 1435689039], [null, 1435689054], [null, 1435689069], [null, 1435689084], [null, 1435689099], [null, 1435689114], [null, 1435689129], [null, 1435689144], [null, 1435689159], [null, 1435689174], [null, 1435689189], [null, 1435689204], [null, 1435689219], [null, 1435689234], [null, 1435689249], [null, 1435689264], [null, 1435689279], [30.79849989414215, 1435689294], [null, 1435689309], [null, 1435689324], [null, 1435689339], [null, 1435689354], [null, 1435689369], [null, 1435689384], [null, 1435689399], [null, 1435689414], [null, 1435689429], [null, 1435689444], [null, 1435689459], [null, 1435689474], [null, 1435689489], [null, 1435689504], [null, 1435689519], [null, 1435689534], [null, 1435689549], [null, 1435689564]]}]

 # assuming the time series is the X axis value

input_X_vars = []
input_Y_vars = []

for pair in time_series[0]["datapoints"]:
    if (pair[0] != None and pair[1] != None):
        input_X_vars.append(pair[1])
        input_Y_vars.append(pair[0])



import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model

regr = linear_model.LinearRegression()
regr.fit(input_X_vars, input_Y_vars)

test_X_vars = [1435688681, 1435688685]

results = regr.predict(test_X_vars)

forecast_append = {"target": "Lower", "datapoints": results}

time_series.append(forecast_append)

Я устанавливаю значение null как None, поскольку ключевое слово 'null' анализируется как просто переменная в python.

Ответ 3

Примечание: это не подробный канонический ответ, а ссылки на доступные библиотеки, которые применяются к вашему домену (статистические модели).

Для python вы можете использовать:

• scipy: scipy.stats.linregress
  • например. http://glowingpython.blogspot.com/2012/03/linear-regression-with-numpy.html (использует scipy.stats.linregress)
• statsmodels: здесь несколько примеров.
• sklearn: здесь документация.

Есть несколько хороших статей:

• http://davidcoallier.com/blog/linear-regression-from-r-to-python/, в частности, он содержит раздел о прогнозировании.
• http://connor-johnson.com/2014/02/18/linear-regression-with-python/, другая статья, показывает как с помощью statsmodels, так и sklearn.

Ответ 4

Я сделал линейную регрессию с использованием градиентного спуска в Python, возможно, это может вам помочь. После инициализации данных создается функция потерь, которая является обычным наименьшим квадратом (OLS). Затем создается функция градиента, и она используется в итеративной процедуре для поиска оптимальных параметров. Затем их можно использовать для прогнозирования. Вот рабочий пример о прогнозе полученного класса, учитывая количество книг и количество часов в классе.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import *
from random import random, seed
from matplotlib import cm

# Initialize datapoints [# Books, hours in class], on the basis of which we want to make a prediction:
#x = np.array([[2104,3],[1600,3],[2400,3],[1416,2],[300,4],[2001,3]])
x = np.array([[0, 9],
[1, 15],
[0, 10],
[2, 16],
[4, 10],
[4, 20],
[1, 11],
[4, 20],
[3, 15],
[0, 15],
[2, 8],
[1, 13],
[4, 18],
[1, 10],
[0, 8],
[1, 10],
[3, 16],
[0, 11],
[1, 19],
[4, 12],
[4, 11],
[0, 19],
[2, 15],
[3, 15],
[1, 20],
[0, 6],
[3, 15],
[3, 19],
[2, 14],
[2, 13],
[3, 17],
[2, 20],
[2, 11],
[3, 20],
[4, 20],
[4, 20],
[3, 9],
[1, 8],
[2, 16],
[0, 10]])

# Initialize y, the grade obtained
y = np.array([45,
57,
45,
51,
65,
88,
44,
87,
89,
59,
66,
65,
56,
47,
66,
41,
56,
37,
45,
58,
47,
64,
97,
55,
51,
61,
69,
79,
71,
62,
87,
54,
43,
92,
83,
94,
60,
56,
88,
62])

#Initialize starting point for theta, the parameter we want to optimize:
#theta = np.array([90,0.1,-8])
theta = np.array([1,1,1])

# Initialize h_theta(x):
def h(x, theta):
    return theta[0] + theta[1]*x[:,0] + theta[2]*x[:,1]

# Initialize the loss function (ordinary least squares), J(theta):
def lossf(theta):
    loss = np.sum((h(x,theta)-y)*(h(x,theta)-y))
    return loss
    '''
    for i in range(x.shape[0]):
        lossIncrement = (h(x[[i]],theta) - y[i])*(h(x[[i]],theta) - y[i])
        loss = loss + lossIncrement

    return 0.5*loss
    '''

# Define the gradient of the loss function for any parameters theta, use vector calculations here instead of for loops:
def gradlossf(theta):
    # dJ/dTheta0:
    d1 = np.sum(h(x,theta)-y)
    d2 = np.sum((h(x,theta)-y)*x[:,0])
    d3 = np.sum((h(x,theta)-y)*x[:,1])
    return np.array([d1,d2,d3])

epsilon = 1e-2
max_iter = 1000
norm = np.linalg.norm(gradlossf(theta))
alpha0 = 1
iter = 0
beta = 1e-4
tolerance = 1e-6
tau = 0.5

# Start iterative procedure:

while iter < max_iter and norm > epsilon:
    print theta

    alpha = alpha0
    alpha_found = False

    while alpha_found is False:
        xpoint = theta - alpha*gradlossf(theta)
        val1 = lossf(xpoint)
        val2 = lossf(theta)

        if val1<val2+tolerance:
            alpha_found = True
            print 'We choose alpha =',alpha
        else:
            alpha = tau*alpha

    iter = iter + 1
    theta = xpoint
    print 'New value of loss function is: ', val1


# THE CORRECT VALUES: [37.379, 4.037, 1.283] (from SPSS) Residual: 7306.230

# Plot datapoints
fig = plt.figure()
plt.hold(True)
x1_surf=np.arange(0, 5, 0.1)                # generate a mesh
x2_surf=np.arange(5, 25, 0.1)
x1_surf, x2_surf = np.meshgrid(x1_surf, x2_surf)
z_surf = theta[0] + theta[1]*x1_surf + theta[2]*x2_surf             # ex. function, which depends on x and y

ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x1_surf, x2_surf, z_surf, cmap=cm.hot);    # plot a 3d surface plot
ax.scatter(x[:,0],x[:,1],y)

ax.set_xlabel('# Books')
ax.set_ylabel('Hours in class')
ax.set_zlabel('Grade')

plt.show()