Trend that saturates on only one end

[bletham]

The logistic growth trend saturates at both lower and upper bounds. It would be useful to have a trend that saturates only on one end, e.g. that saturates at 0 but does not require specifying an upper bound.

This will require adding a new trend, adding to the linear and logistic growth trends that we already have. I think that the logistic loss function $log(1 + exp(k * (t - m)))$ should work very nicely.

If anyone is interested in making a PR for this we can discuss what is involved.

[ibrahimsharaf]

Hi @bletham, I am interested in taking a shot at this, is it a Python or R fix?

[bletham]

This is a little bit of everything - A new Stan model, and then the interfaces in both Python and R. This should be everything that needs to be done:

• Determine the formula for the intercepts of a piecewise logistic-loss function. Basically the "k" in the formula will be offset with deltas, and then "m" will be offset with gammas to maintain continuity. This is described for the logistic loss in Section 3.1.1 of the paper and involves doing the algebra to figure out what gamma needs to be for the segment endpoints to line up.
• Implement a piecewise logistic-loss function that looks like the piecewise logistic function here (In Python and R).
• Implement an initialization for the logistic loss by computing analytically the values for k and m that produce a curve passing through the first and last points in the dataset, like here for logistic growth. (In Python and R).
• Implement the piecewise logistic-loss model in Stan, the one for logistic growth here and here.
• Add it to the list of available models here (Py and R).
• Unit tests and documentation.

This can be done in parts.

[shaidams64]

Hey all, do we have any update on this?

[bletham]

Let's hold off on this until after #501, which will make it much easier and require much less code duplication.

[BernierCR]

What happened on this? #501 doesn't seem to be done yet many months later. So #307 hasn't been started? What is the recommended solution to this problem currently?

[bletham]

@BernierCR the latest Stan a few months ago (2.18) did add a standalone generated quantities functionality that we were planning to use for #501, but it has not yet been piped into rstan or pystan. That piping seems to be going slowly so we are right now looking into handling #501 without waiting for standalone generated quantities. We had tried this last summer and ran into some upstream bugs, hopefully it is successful this time.

In the meantime, you can get that behavior by using the logistic trend but setting the bound on the side that you don't want a bound to just be a a larger number than your forecast ever reaches.

[mmcmm]

@bletham I noticed that the forecast changes based on the upper bond, it seems to scale inversely, the higher the upper bond the smaller the forecast, is this the supposed behaviour?

[bletham]

Could you post the plot so I can see what you mean?

[nartb]

Any update on when this will be available?

[bletham]

#865 is the PR that moves predictions to Stan. There are still some issues there but we're working through them, and once they're resolved this would be next.

[nartb]

Ok, sounds good, thanks! :)

Setting the upper bound to a very high relative value will work though, right?

[germayneng]

hello! just that check the current progress of prophet's feature of not need an upper bound for saturating floor. I currently face the issue of needing to only bound my forecast to 0 yet do not want to enforce an upper bound.

not sure if setting an upper bound will actually alter the intended forecast (even if we set an arbitrary high bound)

[bletham]

No progress, ran into some difficulty in #865 and trying to figure out a good workaround on the rstan side.

[vinodvarma24]

Does anybody have an alternative solution for this meanwhile? because #307 not moving forward.

[bletham]

I was targeting for this to happen after #501, which would have made it easier / more generic, but after the addition of the cmdstanpy backend I've decided to no longer pursue #501 so this should just be done directly.

The new trend function can be added in the same way that the linear and logistic trends are. This is basically what needs to be done, in R and Py:

• Determine the formula for the intercepts of a piecewise logistic-loss function. Basically the "k" in the formula will be offset with deltas, and then "m" will be offset with gammas to maintain continuity. This is described for the logistic loss in Section 3.1.1 of the paper and involves doing the algebra to figure out what gamma needs to be for the segment endpoints to line up.
• Implement a piecewise logistic-loss function that looks like the existing ones:
  

  prophet/python/fbprophet/forecaster.py

  Lines 1197 to 1256 in 46e5611

  	@staticmethod
  	def piecewise_linear(t, deltas, k, m, changepoint_ts):
  	"""Evaluate the piecewise linear function.
  	Parameters
  	----------
  	t: np.array of times on which the function is evaluated.
  	deltas: np.array of rate changes at each changepoint.
  	k: Float initial rate.
  	m: Float initial offset.
  	changepoint_ts: np.array of changepoint times.
  	Returns
  	-------
  	Vector y(t).
  	"""
  	# Intercept changes
  	gammas = -changepoint_ts * deltas
  	# Get cumulative slope and intercept at each t
  	k_t = k * np.ones_like(t)
  	m_t = m * np.ones_like(t)
  	for s, t_s in enumerate(changepoint_ts):
  	indx = t >= t_s
  	k_t[indx] += deltas[s]
  	m_t[indx] += gammas[s]
  	return k_t * t + m_t
  	@staticmethod
  	def piecewise_logistic(t, cap, deltas, k, m, changepoint_ts):
  	"""Evaluate the piecewise logistic function.
  	Parameters
  	----------
  	t: np.array of times on which the function is evaluated.
  	cap: np.array of capacities at each t.
  	deltas: np.array of rate changes at each changepoint.
  	k: Float initial rate.
  	m: Float initial offset.
  	changepoint_ts: np.array of changepoint times.
  	Returns
  	-------
  	Vector y(t).
  	"""
  	# Compute offset changes
  	k_cum = np.concatenate((np.atleast_1d(k), np.cumsum(deltas) + k))
  	gammas = np.zeros(len(changepoint_ts))
  	for i, t_s in enumerate(changepoint_ts):
  	gammas[i] = (
  	(t_s - m - np.sum(gammas))
  	* (1 - k_cum[i] / k_cum[i + 1]) # noqa W503
  	)
  	# Get cumulative rate and offset at each t
  	k_t = k * np.ones_like(t)
  	m_t = m * np.ones_like(t)
  	for s, t_s in enumerate(changepoint_ts):
  	indx = t >= t_s
  	k_t[indx] += deltas[s]
  	m_t[indx] += gammas[s]
  	return cap / (1 + np.exp(-k_t * (t - m_t)))

  
  

  prophet/R/R/prophet.R

  Lines 1354 to 1410 in 46e5611

  	#' Evaluate the piecewise linear function.
  	#'
  	#' @param t Vector of times on which the function is evaluated.
  	#' @param deltas Vector of rate changes at each changepoint.
  	#' @param k Float initial rate.
  	#' @param m Float initial offset.
  	#' @param changepoint.ts Vector of changepoint times.
  	#'
  	#' @return Vector y(t).
  	#'
  	#' @keywords internal
  	piecewise_linear <- function(t, deltas, k, m, changepoint.ts) {
  	# Intercept changes
  	gammas <- -changepoint.ts * deltas
  	# Get cumulative slope and intercept at each t
  	k_t <- rep(k, length(t))
  	m_t <- rep(m, length(t))
  	for (s in 1:length(changepoint.ts)) {
  	indx <- t >= changepoint.ts[s]
  	k_t[indx] <- k_t[indx] + deltas[s]
  	m_t[indx] <- m_t[indx] + gammas[s]
  	}
  	y <- k_t * t + m_t
  	return(y)
  	}
  	#' Evaluate the piecewise logistic function.
  	#'
  	#' @param t Vector of times on which the function is evaluated.
  	#' @param cap Vector of capacities at each t.
  	#' @param deltas Vector of rate changes at each changepoint.
  	#' @param k Float initial rate.
  	#' @param m Float initial offset.
  	#' @param changepoint.ts Vector of changepoint times.
  	#'
  	#' @return Vector y(t).
  	#'
  	#' @keywords internal
  	piecewise_logistic <- function(t, cap, deltas, k, m, changepoint.ts) {
  	# Compute offset changes
  	k.cum <- c(k, cumsum(deltas) + k)
  	gammas <- rep(0, length(changepoint.ts))
  	for (i in 1:length(changepoint.ts)) {
  	gammas[i] <- ((changepoint.ts[i] - m - sum(gammas))
  	* (1 - k.cum[i] / k.cum[i + 1]))
  	}
  	# Get cumulative rate and offset at each t
  	k_t <- rep(k, length(t))
  	m_t <- rep(m, length(t))
  	for (s in 1:length(changepoint.ts)) {
  	indx <- t >= changepoint.ts[s]
  	k_t[indx] <- k_t[indx] + deltas[s]
  	m_t[indx] <- m_t[indx] + gammas[s]
  	}
  	y <- cap / (1 + exp(-k_t * (t - m_t)))
  	return(y)
  	}

• Implement an initialization by computing analytically the values for k and m that produce a curve passing through the first and last points in the dataset, like here:
  

  prophet/python/fbprophet/forecaster.py

  Lines 986 to 1050 in 46e5611

  	@staticmethod
  	def linear_growth_init(df):
  	"""Initialize linear growth.
  	Provides a strong initialization for linear growth by calculating the
  	growth and offset parameters that pass the function through the first
  	and last points in the time series.
  	Parameters
  	----------
  	df: pd.DataFrame with columns ds (date), y_scaled (scaled time series),
  	and t (scaled time).
  	Returns
  	-------
  	A tuple (k, m) with the rate (k) and offset (m) of the linear growth
  	function.
  	"""
  	i0, i1 = df['ds'].idxmin(), df['ds'].idxmax()
  	T = df['t'].iloc[i1] - df['t'].iloc[i0]
  	k = (df['y_scaled'].iloc[i1] - df['y_scaled'].iloc[i0]) / T
  	m = df['y_scaled'].iloc[i0] - k * df['t'].iloc[i0]
  	return (k, m)
  	@staticmethod
  	def logistic_growth_init(df):
  	"""Initialize logistic growth.
  	Provides a strong initialization for logistic growth by calculating the
  	growth and offset parameters that pass the function through the first
  	and last points in the time series.
  	Parameters
  	----------
  	df: pd.DataFrame with columns ds (date), cap_scaled (scaled capacity),
  	y_scaled (scaled time series), and t (scaled time).
  	Returns
  	-------
  	A tuple (k, m) with the rate (k) and offset (m) of the logistic growth
  	function.
  	"""
  	i0, i1 = df['ds'].idxmin(), df['ds'].idxmax()
  	T = df['t'].iloc[i1] - df['t'].iloc[i0]
  	# Force valid values, in case y > cap or y < 0
  	C0 = df['cap_scaled'].iloc[i0]
  	C1 = df['cap_scaled'].iloc[i1]
  	y0 = max(0.01 * C0, min(0.99 * C0, df['y_scaled'].iloc[i0]))
  	y1 = max(0.01 * C1, min(0.99 * C1, df['y_scaled'].iloc[i1]))
  	r0 = C0 / y0
  	r1 = C1 / y1
  	if abs(r0 - r1) <= 0.01:
  	r0 = 1.05 * r0
  	L0 = np.log(r0 - 1)
  	L1 = np.log(r1 - 1)
  	# Initialize the offset
  	m = L0 * T / (L0 - L1)
  	# And the rate
  	k = (L0 - L1) / T
  	return (k, m)

  
  

  prophet/R/R/prophet.R

  Lines 1089 to 1150 in 46e5611

  	#' Initialize linear growth.
  	#'
  	#' Provides a strong initialization for linear growth by calculating the
  	#' growth and offset parameters that pass the function through the first and
  	#' last points in the time series.
  	#'
  	#' @param df Data frame with columns ds (date), y_scaled (scaled time series),
  	#' and t (scaled time).
  	#'
  	#' @return A vector (k, m) with the rate (k) and offset (m) of the linear
  	#' growth function.
  	#'
  	#' @keywords internal
  	linear_growth_init <- function(df) {
  	i0 <- which.min(df$ds)
  	i1 <- which.max(df$ds)
  	T <- df$t[i1] - df$t[i0]
  	# Initialize the rate
  	k <- (df$y_scaled[i1] - df$y_scaled[i0]) / T
  	# And the offset
  	m <- df$y_scaled[i0] - k * df$t[i0]
  	return(c(k, m))
  	}
  	#' Initialize logistic growth.
  	#'
  	#' Provides a strong initialization for logistic growth by calculating the
  	#' growth and offset parameters that pass the function through the first and
  	#' last points in the time series.
  	#'
  	#' @param df Data frame with columns ds (date), cap_scaled (scaled capacity),
  	#' y_scaled (scaled time series), and t (scaled time).
  	#'
  	#' @return A vector (k, m) with the rate (k) and offset (m) of the logistic
  	#' growth function.
  	#'
  	#' @keywords internal
  	logistic_growth_init <- function(df) {
  	i0 <- which.min(df$ds)
  	i1 <- which.max(df$ds)
  	T <- df$t[i1] - df$t[i0]
  	# Force valid values, in case y > cap or y < 0
  	C0 <- df$cap_scaled[i0]
  	C1 <- df$cap_scaled[i1]
  	y0 <- max(0.01 * C0, min(0.99 * C0, df$y_scaled[i0]))
  	y1 <- max(0.01 * C1, min(0.99 * C1, df$y_scaled[i1]))
  	r0 <- C0 / y0
  	r1 <- C1 / y1
  	if (abs(r0 - r1) <= 0.01) {
  	r0 <- 1.05 * r0
  	}
  	L0 <- log(r0 - 1)
  	L1 <- log(r1 - 1)
  	# Initialize the offset
  	m <- L0 * T / (L0 - L1)
  	# And the rate
  	k <- (L0 - L1) / T
  	return(c(k, m))
  	}

• Implement the piecewise logistic-loss model in Stan, like here:
  

  prophet/python/stan/unix/prophet.stan

  Lines 48 to 76 in 46e5611

  	vector logistic_trend(
  	real k,
  	real m,
  	vector delta,
  	vector t,
  	vector cap,
  	matrix A,
  	vector t_change,
  	int S
  	) {
  	vector[S] gamma;
  	gamma = logistic_gamma(k, m, delta, t_change, S);
  	return cap .* inv_logit((k + A * delta) .* (t - (m + A * gamma)));
  	}
  	// Linear trend function
  	vector linear_trend(
  	real k,
  	real m,
  	vector delta,
  	vector t,
  	matrix A,
  	vector t_change
  	) {
  	return (k + A * delta) .* t + (m + A * (-t_change .* delta));
  	}
  	}

• Add it to the places where we switch the trend (these are the Py places, corresponding places are in R):
  

  prophet/python/stan/unix/prophet.stan

  Line 116 in 46e5611

  if (trend_indicator == 0) {

  
  

  prophet/python/fbprophet/forecaster.py

  Line 1111 in 46e5611

  'trend_indicator': int(self.growth == 'logistic'),

  
  

  prophet/python/fbprophet/forecaster.py

  Lines 1116 to 1121 in 46e5611

  	if self.growth == 'linear':
  	dat['cap'] = np.zeros(self.history.shape[0])
  	kinit = self.linear_growth_init(history)
  	else:
  	dat['cap'] = history['cap_scaled']
  	kinit = self.logistic_growth_init(history)

  
  

  prophet/python/fbprophet/forecaster.py

  Lines 1274 to 1279 in 46e5611

  	if self.growth == 'linear':
  	trend = self.piecewise_linear(t, deltas, k, m, self.changepoints_t)
  	else:
  	cap = df['cap_scaled']
  	trend = self.piecewise_logistic(
  	t, cap, deltas, k, m, self.changepoints_t)

  
  

  prophet/python/fbprophet/forecaster.py

  Lines 1471 to 1476 in 46e5611

  	if self.growth == 'linear':
  	trend = self.piecewise_linear(t, deltas, k, m, changepoint_ts)
  	else:
  	cap = df['cap_scaled']
  	trend = self.piecewise_logistic(t, cap, deltas, k, m,
  	changepoint_ts)

• Unit tests and documentation.

[domarcon]

How does inputing an arbitrary large value on cap would affect the results?

[bletham]

I think one of the main applications of this is data that has to be positive but isn't necessarily saturating at an upper limit. There's discussion of that setting in #1668 along with some strategies so I just wanted to point people to that issue in case anyone wants to try the approaches described there. I'm especially interested to hear if the ProphetPos class that is described in that issue provides a satisfactory solution to the use cases here.

[bletham]

Following up on my last comment: The documentation now provides examples of how to add custom trends in both R and Py (#1466, #1794, #1778). For now we are not going to add any new trend functions; there is a long tail of trends that one might be interested in, and we don't want to add that much complexity/weight to the package. So the PRs listed above will provide a guide for those who would like to try something totally custom, including a trend that saturates on only one end.

I suspect the main source of interest in this issue is to produce forecasts that saturate at 0 and thus stay positive. There is discussion on that problem in #1668. A trend that saturates at 0 would not be sufficient by itself to produce positive forecasts, so I think that problem should be discussed separately from this.