Is there a piece of software that can help me compute the expectation of many piecewise-linear functions?

Question

Take a piecewise-linear function

$$f(x,y)=\begin{cases} ax+b \quad x\leq y\\ cx+d \quad x > y \end{cases}$$

I need to obtain an expression of $\mathbb{E}_{y}f(x,y)$, knowing that there are only finitely many $y$ values, say $y_1,\ldots,y_N$ and their own probabilities $p_1,\ldots,p_N$.

The expectation will itself be a piecewise-linear function. If $N$ is fairly small, say $2$ or $3$, I can manually do the calculations and obtain an expression of the expectation. However, as $N$ grows, this task soon becomes impractical. I could probably write my own code in e.g., Python, to do that, but I feel there is a chance I am reinventing the wheel. So, here is the question: Is there a piece of (possibly free) software that can give me the expression of $g(x) = \mathbb{E}_{y}f(x,y) = \sum_{i=1}^Np_if(x,y_i)$?

Answer 1

Here's a function in R that I wrote for $g(x)$:

# Expectation of many piecewise function
g <- function(x){
  a <- 2
  b <- 3
  c <- -5
  d <- 7
  y <- c(-8,-4,1,9) # should be ascending
  p <- c(0.1,0.2,0.3,0.4) # should sum to 1

  if ( !identical(order(y),1:length(y)) ) {
    stop("y vector in definition of g is not ascending")
  }

  if (sum(p)!=1) {
    stop("probabilities do not sum to one")
  }

  if ( length(p) != length(y) ) {
    stop("length of p does not equal length of y")
  }

  # return value
  (1-sum((y<=x)*p))*(a*x+b) + sum((y<=x)*p)*(c*x+d)
}

# Plot of g(x)
ordinate <- numeric(length = length((-150:150)/10))
for (i in -150:150) {
  ordinate[[i+151]] <- g(i/10)
}
plot((-150:150)/10, ordinate, xlab = "x", ylab = "g(x)")