पायथन + आर का उपयोग कर एक SARIMA मॉडल का निर्माण

परिचय

शुभ दोपहर, प्रिय पाठकों।
पायथन में टाइम सीरीज़ के विश्लेषण के बारे में पिछली पोस्ट लिखने के बाद, मैंने उन टिप्पणियों को ठीक करने का फैसला किया, जो टिप्पणियों में इंगित की गई थीं, लेकिन जब उन्हें ठीक किया गया, तो मैं कई समस्याओं में भाग गया, उदाहरण के लिए, जब एक मौसमी ARIMA मॉडल का निर्माण किया गया, क्योंकि मुझे स्टैटमॉडल पैकेज में एक समान फ़ंक्शन नहीं मिला। अंत में, मैंने इसके लिए R से फ़ंक्शन का उपयोग करने का निर्णय लिया, और खोजों ने मुझे rpy2 लाइब्रेरी में ले गया, जो आपको उल्लिखित भाषा के लाइब्रेरी से फ़ंक्शंस का उपयोग करने की अनुमति देता है।
कई पूछ सकते हैं, "यह आवश्यक क्यों है?", क्योंकि आर को लेना और उसमें सभी काम करना आसान है। मैं इस कथन से पूरी तरह सहमत हूं, लेकिन यह मुझे प्रतीत होता है कि अगर डेटा को प्रारंभिक प्रसंस्करण की आवश्यकता होती है, तो इसे पायथन में उत्पादन करना आसान होता है, और विश्लेषण के लिए आवश्यक होने पर आर क्षमताओं का उपयोग किया जा सकता है।
इसके अलावा, यह दिखाएगा कि फ़ंक्शन आर के आउटपुट के परिणामों को आईपीथॉन नोटबुक में कैसे एकीकृत किया जाए ।

Rpy2 को स्थापित और कॉन्फ़िगर करें

आरंभ करने के लिए, आपको rpy2 स्थापित करना होगा। आप कमांड का उपयोग करके ऐसा कर सकते हैं:

pip install rpy2 

यह ध्यान दिया जाना चाहिए कि, इस पुस्तकालय को काम करने के लिए, स्थापित आर भाषा की आवश्यकता है। आप इसे बंद से डाउनलोड कर सकते हैं । साइट ।
अगली बात यह है कि आवश्यक सिस्टम चर जोड़ें। विंडोज के लिए, निम्नलिखित चर जोड़ें:

• पथ - R.dll और R.exe के लिए पथ
  
• R_HOME - उस फ़ोल्डर का पथ जिसमें R स्थापित है
  
• R_USER - उपयोगकर्ता नाम जिसके तहत विंडो लोड की गई है
  

मैक ओएस एक्स पर इंस्टॉलेशन भी किया गया था, जहां कोई हेरफेर डेटा की आवश्यकता नहीं थी।

शुरुआत हो रही है

इसलिए, यदि आप आईपीथॉन नोटबुक में काम कर रहे हैं, तो आपको एक निर्देश जोड़ने की आवश्यकता है:

 %load_ext rmagic 

यह एक्सटेंशन आपको rpy2 के माध्यम से कुछ आर फ़ंक्शन को कॉल करने की अनुमति देता है, और परिणाम को सीधे IPython नोटबुक कंसोल में प्रदर्शित करता है, जो बहुत सुविधाजनक है (नीचे हम यह दिखाएंगे कि यह कैसे करना है)। अधिक विवरण यहां लिखे गए हैं ।
अब आवश्यक लाइब्रेरी डाउनलोड करें:

 from pandas import read_csv, DataFrame, Series import statsmodels.api as sm import rpy2.robjects as R from rpy2.robjects.packages import importr import pandas.rpy.common as com from pandas import date_range 

अब, पिछले लेख की तरह, डेटा लोड करें और साप्ताहिक अंतराल पर आगे बढ़ें:

 dataset = read_csv('DataSets/tovar_moving.csv',';', index_col=['date'], parse_dates=['date'], dayfirst=True) otg = dataset.qty w = otg.resample('w', how='sum') w.plot(figsize=(12,6)) 

तो, ग्राफ से आप वार्षिक मौसमी (52 सप्ताह) और एक स्पष्ट प्रवृत्ति देख सकते हैं। इसलिए, मॉडल के निर्माण से पहले, हमें प्रवृत्ति और मौसमी से छुटकारा पाने की आवश्यकता है।

प्रारंभिक डेटा विश्लेषण

इसलिए, शुरुआत के लिए, हम मूल्यों को संरेखित करने के लिए मूल श्रृंखला का प्रस्ताव करेंगे:

 w_log = log(w) w_log.plot(figsize=(12,6)) 

जैसा कि आप देख सकते हैं, ग्राफ में मौसमी है और तदनुसार, श्रृंखला स्थिर नहीं है। हम डिक्की-फुलर परीक्षण का उपयोग करके इसे सत्यापित करेंगे, जो इकाई जड़ों की उपस्थिति की परिकल्पना की जांच करता है और तदनुसार, यदि वे हैं, तो श्रृंखला स्थिर नहीं होगी। आँकड़ों की लाइब्रेरी का उपयोग करके इस परीक्षा का संचालन कैसे करें, मैंने पिछली बार दिखाया था। अब मैं प्रदर्शित करता हूं कि यह कैसे आर से adf.test () फ़ंक्शन का उपयोग करके किया जा सकता है ।
तो, यह फ़ंक्शन ट्रेसिस आर-लाइब्रेरी में है। यह समय श्रृंखला विश्लेषण के लिए डिज़ाइन किया गया है और वैकल्पिक है। आप आयातकर्ता () फ़ंक्शन का उपयोग करके वांछित लाइब्रेरी लोड कर सकते हैं।

 stats = importr('stats') tseries = importr('tseries') 

आप देख सकते हैं कि ट्रेसियों के अलावा, हमने सांख्यिकी पुस्तकालय भी डाउनलोड किया है। टाइप रूपांतरण के लिए हमें इसकी आवश्यकता है।
अब आपको पायथन प्रकार से डेटा को समझने योग्य R प्रकार में कनवर्ट करने की आवश्यकता है। आप Data_rame को खिलाए गए इनपुट के लिए Convert_to_r_dataframe () फ़ंक्शन का उपयोग करके ऐसा कर सकते हैं, और आउटपुट R के लिए एक वेक्टर है।

 r_df = com.convert_to_r_dataframe(DataFrame(w_log)) 

तो, वेक्टर अगला चरण है, आपको इसे समय श्रृंखला के प्रारूप में अनुवाद करने की आवश्यकता है। इसके लिए, आर में एक फ़ंक्शन ts () है, इसकी कॉल इस तरह दिखाई देगी:

 y = stats.ts(r_df) 

प्रारंभिक डेटा तैयारी पूरी हो गई है और हम उस फ़ंक्शन को कॉल कर सकते हैं जिसकी हमें आवश्यकता है:

 ad = tseries.adf_test(y, alternative="stationary", k=52) 

मापदंडों के रूप में, यह समय श्रृंखला और अंतराल की संख्या है जिसके लिए परीक्षण की गणना की जाएगी। आर्थिक मॉडल में, इस मूल्य को वर्ष के बराबर लेने की प्रथा है, और चूंकि हमारे पास साप्ताहिक डेटा और एक वर्ष में 52 सप्ताह हैं, इसलिए पैरामीटर का ऐसा मूल्य है।
विज्ञापन चर में अब एक R ऑब्जेक्ट है। इसकी संरचना एक सूची के रूप में वर्णित है, जिसका वर्णन मुझे नहीं मिला। इसलिए, दृश्य विश्लेषण का उपयोग करते हुए, मैंने एक कोड लिखा जो फ़ंक्शन के परिणाम को समझ में आता है:

 a = ad.names[:5] {ad.names[i]:ad[i][0] for i in xrange(len(a))} 

{'वैकल्पिक': 'स्थिर',
'विधि': 'संवर्धित डिकी-फुलर टेस्ट',
'p.value': 0.23867869477446427,
'पैरामीटर': 52.0,
'आँकड़ा': -2.8030060277420006}

परीक्षण के परिणामों के आधार पर, प्रारंभिक श्रृंखला स्थिर नहीं है। क्योंकि इकाई जड़ों की उपस्थिति की परिकल्पना को कम संभावना के साथ स्वीकार किया जाता है, और, तदनुसार, श्रृंखला स्थिर नहीं होती है। अब स्टेशनरी के लिए पहले कई अंतरों की जाँच करें।
पहले, उन्हें पायथन का उपयोग करके प्राप्त करें, और फिर ADF परीक्षण लागू करें:

 diff1lev = w.diff(periods=1).dropna() print 'p.value: %f' % sm.tsa.adfuller(diff1lev, maxlag=52)[1] 

p.value: 0.000000

 diff1lev.plot(figsize=(12,6)) 

परीक्षण के परिणामों के अनुसार, पहले अंतर की एक संख्या स्थिर हो गई। और ग्राफ हमें यह सुनिश्चित करने में मदद करता है कि कोई प्रवृत्ति नहीं है। यह मौसमी से छुटकारा पाने के लिए बनी हुई है।
ऐसा करने के लिए, आपको हमारी परिणामी श्रृंखला से मौसमी अंतर को लेना होगा। यहाँ और पढ़ें यदि परिणामस्वरूप पंक्ति स्थिर है, तो इसके पहले अंतर को बुनना और इसकी जांच करना आवश्यक होगा।

 diff1lev_season = diff1lev.diff(52).dropna() 

आइए R से ADF परीक्षण का उपयोग करते हुए इसे स्टेशनरी के लिए जांचें:

 r_df = com.convert_to_r_dataframe(DataFrame(diff1lev_season)) y = stats.ts(r_df) ad = tseries.adf_test(y, alternative="stationary", k=52) a = ad.names[:5] {ad.names[i]:ad[i][0] for i in xrange(len(a))} 

{'वैकल्पिक': 'स्थिर',
'विधि': 'संवर्धित डिकी-फुलर टेस्ट',
'p.value': 0.551977997289418,
'पैरामीटर': 52.0,
'आँकड़ा': -2.0581183466564776}

तो, श्रृंखला स्थिर नहीं है, इसका पहला अंतर लें:

 diff1lev_season1lev = diff1lev_season.diff().dropna() print 'p.value: %f' % sm.tsa.adfuller(diff1lev_season1lev, maxlag=52)[1] 

p.value: 0.000395

परिणामस्वरूप पंक्ति स्थिर है। अब आप मॉडल के निर्माण के लिए आगे बढ़ सकते हैं

एक मॉडल का निर्माण।

तो, प्रारंभिक विश्लेषण पूरा हो गया है, और हम मौसमी मॉडल ARIMA (SARIMA) के निर्माण के लिए आगे बढ़ सकते हैं।
इस मॉडल का सामान्य दृश्य
इस मॉडल में, पैरामीटर निम्नलिखित दर्शाते हैं:

• - मॉडल ऑर्डर
• - स्रोत डेटा के एकीकरण का क्रम
• - मॉडल ऑर्डर
• - मौसमी घटक का क्रम
• - मौसमी घटक को एकीकृत करने की प्रक्रिया
• - मौसमी घटक का क्रम
• - मौसमी का आयाम (माह, तिमाही, आदि)

पी , डी , क्यू का निर्धारण कैसे करें, मैंने पिछली बार दिखाया था। अब मैं पी , डी , क्यू के मौसमी घटकों के क्रम को निर्धारित करने का वर्णन करूंगा ।
हम पैरामीटर डी को परिभाषित करके शुरू करते हैं । यह मौसमी अंतर के एकीकरण के क्रम को निर्धारित करेगा, अर्थात। हमारे मामले में, यह 1 के बराबर है। P और Q को निर्धारित करने के लिए , हमें अभी भी ACR और PACF correlograms के निर्माण की आवश्यकता है।

 fig = plt.figure(figsize=(12,8)) ax1 = fig.add_subplot(211) fig = sm.graphics.tsa.plot_acf(diff1lev_season1lev.values.squeeze(), lags=150, ax=ax1) ax2 = fig.add_subplot(212) fig = sm.graphics.tsa.plot_pacf(diff1lev_season1lev, lags=150, ax=ax2) 

पीएसीएफ गारफिक से यह देखा जा सकता है कि एआर ऑर्डर पी = 4 होगा , और एसीएफ के अनुसार, ऑर्डर एमए क्यू = 13 , क्योंकि 13 लैग 0 के अलावा अन्य अंतिम अंतराल है।
अब हम मौसमी घटकों पर चलते हैं। उनका मूल्यांकन करने के लिए, आपको उन लैग्स को देखने की जरूरत है जो मौसमी के आकार के गुणक हैं, यानी, यदि, हमारे उदाहरण के लिए, मौसमी 52 है, तो लैग 52, 104, 156, ...
हमारे मामले में, पैरामीटर पी और क्यू 0 के बराबर होगा (यह देखा जा सकता है यदि आप उपरोक्त लैग पर एसीएफ और पीएसीएफ को देखते हैं)।
हमारे शोध के परिणामस्वरूप, हमें एक मॉडल मिला
जैसा कि इस लेख की शुरुआत में संकेत दिया गया था, मुझे पायथन में इस मॉडल को बनाने के तरीके नहीं मिले, इसलिए मैंने आर से एरीमा () फ़ंक्शन का उपयोग करने का फैसला किया। इसके लिए, एआरआईएमए मॉडल का क्रम और, यदि आवश्यक हो, तो मौसमी घटक के आदेश को पारित कर दिया जाता है। । लेकिन इसे कॉल करने से पहले, आपको कुछ डेटा तैयार करने की आवश्यकता है।
आरंभ करने के लिए, हम अपने स्रोत को R प्रारूप में सेट करते हैं और समय श्रृंखला प्रारूप में अनुवाद करते हैं।

 r_df = com.convert_to_r_dataframe(DataFrame(w)) y = stats.ts(r_df) 

मॉडल का क्रम वेक्टर आर के रूप में पारित किया गया है, तो चलिए इसे बनाते हैं:

 order = R.IntVector((4,1,13)) 

इसके अलावा, मौसमी घटक के एक पैरामीटर के रूप में, एक सूची पारित की जाती है जिसमें इसका क्रम और अवधि का आकार होता है:

 season = R.ListVector({'order': R.IntVector((0,1,0)), 'period' : 52}) 

अब हम एक मॉडल बनाने के लिए तैयार हैं:

 model = stats.arima(y, order = order, seasonal=season) 

इसलिए हमारा मॉडल तैयार है और हम इसके आधार पर पूर्वानुमान के निर्माण के लिए आगे बढ़ सकते हैं।

मॉडल की पर्याप्तता की जाँच करना।

इसलिए, मॉडल की पर्याप्तता की जांच करने के लिए, आपको यह जांचने की आवश्यकता है कि क्या मॉडल के अवशेष "सफेद शोर" के अनुरूप हैं। हम फेफड़े - बॉक्स क्यू परीक्षण करके इसे सत्यापित करते हैं और अवशेषों के सहसंबंध को सत्यापित करते हैं। ऐसा करने के लिए, आर में एक tsdiag () फ़ंक्शन है, परीक्षण के लिए मॉडल और लैग्स की संख्या एक पैरामीटर के रूप में पारित की जाती है।
आप इस फ़ंक्शन को इस तरह से कॉल कर सकते हैं:

 %Rpush model %R tsdiag(model, 100) 

पहली पंक्ति में, % Rpush निर्देश आर में उपयोग के लिए वस्तुओं को लोड करता है। दूसरी पंक्ति में % R कथन R भाषा प्रारूप में कोड को कॉल करता है। यह निर्माण IPython नोटबुक में काम करता है।
ऊपर दिए गए ग्राफ़ से, आप देख सकते हैं कि अवशिष्ट स्वतंत्र हैं (इसे ACF से देखा जा सकता है)। इसके अलावा, क्यू-सांख्यिकी के ग्राफ से, आप देख सकते हैं कि सभी बिंदुओं पर पी-वैल्यू महत्व स्तर से अधिक है, इससे हम निष्कर्ष निकाल सकते हैं कि अवशेष "सफेद शोर" होने की सबसे अधिक संभावना है।

भविष्यवाणी

इसे चलाने के लिए, आपको पूर्वानुमान पुस्तकालय को लोड करने की आवश्यकता है

 forecast = importr('forecast') 

पूर्वानुमान परिणाम प्राप्त करने के दो तरीके हैं।

विधि 1. यह IPython और R. की एकीकरण क्षमताओं का उपयोग करना है, जो पिछले भाग में दिखाया गया था:

 %R plot(forecast(model)) 

विधि 2. दूसरी विधि पूर्वानुमान लाइब्रेरी का उपयोग करके एक पूर्वानुमान बनाने के लिए है, और फिर परिणाम को एक अस्थायी पांडा श्रृंखला में अनुवाद करें और उन्हें स्क्रीन पर प्रदर्शित करें। जिस कोड को निष्पादित किया जाएगा, वह नीचे लिखा गया है:

 f = forecast.forecast(model) #  pred = [i[1] for i in f[3].iteritems()] #  dt = date_range(w.index[-1], periods=len(pred)+1,freq='W')[1:] #    pr = Series(pred, index = dt) #   TimeSeries 

अब परिणाम प्रदर्शित करें

 w.plot(figsize=(12,6)) pr.plot(color = 'red') 

निष्कर्ष

अंत में, मैं यह नोट करना चाहूंगा कि मौसमी के अधिक सटीक विश्लेषण के लिए, आपको 7-10 मौसमों के साथ डेटा होना चाहिए। इसके अलावा, मैं लेख तैयार करने की प्रक्रिया में सहायता के लिए वेयरवोल्फ को एक बड़ा धन्यवाद कहना चाहूंगा।
इस लेख में मैंने यह दिखाने की कोशिश की कि कैसे एक मौसमी ARIMA मॉडल का निर्माण किया जाता है, और यह भी दिखाया गया है कि आप डेटा विश्लेषण के लिए R और पायथन भाषाओं के एक समूह का उपयोग कैसे कर सकते हैं। मुझे लगता है कि एक और प्रेमिका भाषा के विशेषज्ञ पाएंगे कि वर्णित गुच्छा को प्रभावी ढंग से कैसे लागू किया जाए।