1+ import numpy as np
2+ import pandas as pd
3+ from sklearn .preprocessing import StandardScaler , PolynomialFeatures
4+ from sklearn .pipeline import make_pipeline
5+ import statsmodels .api as sm
6+ from statsmodels .sandbox .regression .predstd import wls_prediction_std
7+
8+
9+ def check_tolerance (t , y , to_exclude : int = 1 , poly_features : list = [1 , 2 ], alpha : float = 0.05 , predict_all : bool = False ) -> pd .DataFrame :
10+ """
11+ Check that some future values are within a weighted least squares confidence interval.
12+
13+ Parameters
14+ ----------
15+ t : np.array
16+ N explanatory time points of shape (N, 1).
17+ y : np.array
18+ The corresponding response variable values to X, of shape (N, 1).
19+ to_exclude : int, default = 1
20+ How many of the last y values will have their tolerances checked.
21+ poly_features : list, default = [1, 2]
22+ List of degrees of polynomial basis to fit to the data. One model will be
23+ produced for each number in the list, eg. the default will fit a linear and
24+ a second degree polynomial to the data and return both sets of results.
25+ alpha : float, default = 0.05
26+ Alpha parameter for the weighted least squares confidence interval.
27+ predict_all : bool, default = False
28+ Set to true to show predictions for all points of the dataset.
29+
30+
31+ Returns
32+ -------
33+ pd.DataFrame
34+ DataFrame containing:
35+ "yhat_u" : Upper condfidence interval for y
36+ "yobs" : Observed value for y
37+ "yhat" : Predicted value for y
38+ "yhat_l" : Lower confidence interval for y
39+ "polynomial": Max polynomial of model fit to the data
40+
41+
42+ Examples
43+ --------
44+ >>> check_tolerance(
45+ ... t = np.array([1001,1002,1003,1004,1005,1006]),
46+ ... y = np.array([2,3,4,4.5,5,5.1]),
47+ ... to_exclude = 2,
48+ ... )
49+ yhat_u yobs yhat yhat_l polynomial
50+ 0 6.817413 5.0 5.500 4.182587 1
51+ 1 7.952702 5.1 6.350 4.747298 1
52+ 2 9.077182 5.0 4.875 0.672818 2
53+ 3 13.252339 5.1 4.975 -3.302339 2
54+ """
55+
56+ if not isinstance (poly_features , list ):
57+ raise ValueError ("Please input a list of integers from 0 to 4 for poly_features." )
58+ assert all (0 <= degree <= 4 for degree in poly_features ), (
59+ "Please ensure all numbers in poly_features are from 0 to 4."
60+ )
61+ if not isinstance (alpha , float ) or 0 > alpha >= 1 :
62+ raise ValueError ("Please input a float between 0 and 1 for alpha." )
63+ if not isinstance (to_exclude , int ):
64+ raise ValueError ("Please input an integer between 1 and your sample size for to_exclude." )
65+ assert ((len (t ) - to_exclude ) >= 4 ), (
66+ """The sample size for your model is smaller than 4. This will not produce a good
67+ model. Either reduce to_exclude or increase your sample size to continue."""
68+ )
69+ assert np .isfinite (y ).all (), (
70+ f"""Your sample contains missing or infinite values for y at locations
71+ { list (map (tuple , np .where (np .isnan (y ))))} . Exclude these values to continue."""
72+ )
73+ assert np .isfinite (t ).all (), (
74+ f"""Your sample contains missing or infinite values for t at locations
75+ { list (map (tuple , np .where (np .isnan (t ))))} . Exclude these values to continue."""
76+ )
77+
78+
79+ # Sort data by X increasing
80+ idx = np .argsort (t )
81+ t = t [idx ]
82+ y = y [idx ]
83+
84+ results = pd .DataFrame ()
85+ for degree in poly_features :
86+ transforms = make_pipeline (
87+ StandardScaler (),
88+ PolynomialFeatures (degree = degree ),
89+ )
90+
91+ # Fit transforms to training data only, apply to all data.
92+ fitted_transforms = transforms .fit (t [:- to_exclude ].reshape (- 1 , 1 ))
93+ _t = fitted_transforms .transform (t .reshape (- 1 , 1 ))
94+
95+ t_train , y_train = _t [:- to_exclude , :], y [:- to_exclude ]
96+ t_predict , y_predict = (
97+ _t if predict_all else _t [- to_exclude :, :],
98+ y if predict_all else y [- to_exclude :]
99+ )
100+
101+ # Fit ordinary least squares model to the training data, then predict for the
102+ # prediction data.
103+ model = sm .OLS (y_train , t_train ).fit ()
104+ yhat = model .predict (t_predict )
105+
106+ # Calculate prediction intervals of fitted model.
107+ _ , yhat_l , yhat_u = wls_prediction_std (model , t_predict , alpha = alpha )
108+
109+ # Store model results in master frame
110+ results = results .append (
111+ pd .DataFrame ({
112+ "yhat_u" : yhat_u ,
113+ "yobs" : y_predict ,
114+ "yhat" : yhat ,
115+ "yhat_l" : yhat_l ,
116+ "polynomial" : degree
117+ }),
118+ ignore_index = True ,
119+ )
120+
121+ return results
122+
0 commit comments