--- title: "Polynomial Library in BsplineQuantReg" author: "Alexandre Abbes" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 3 fig_width: 7 fig_height: 5 vignette: > %\VignetteIndexEntry{Polynomial Library} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, fig.width = 7, fig.height = 5) library(BsplineQuantReg) ``` ## Introduction The `BsplineQuantReg` package provides a comprehensive set of tools for polynomial manipulation, implemented in pure R with a consistent interface. All polynomial functions follow the **decreasing power order convention**, where a polynomial $p(x) = a_0 + a_1 x + a_2 x^2$ is represented as `c(a_2, a_1, a_0)`. This vignette covers: - Basic polynomial operations (addition, multiplication, evaluation, differentiation) - Piecewise polynomial (PP) form - Callable and non-callable PP objects - Extracting parameters from PP objects ## Basic Polynomial Operations ### Polynomial Multiplication ```{r polymul} # (1 + x) * (1 + x) = 1 + 2x + x^2 p1 <- c(1, 1) # 1 + x p2 <- c(1, 1) # 1 + x product <- polymul(p1, p2) print(product) # c(1, 2, 1) → 1 + 2x + x^2 # (x^2 + 2x + 1) * (x - 1) = x^3 + x^2 - x - 1 p3 <- c(1, 2, 1) # x^2 + 2x + 1 p4 <- c(1, -1) # x - 1 polymul(p3, p4) ``` ### Polynomial Addition ```{r polyadd} # (1 + x) + (1 - x) = 2 polyadd(c(1, 1), c(1, -1)) # (x^2 + 1) + (x + 1) = x^2 + x + 2 polyadd(c(1, 0, 1), c(1, 1)) ``` ### Polynomial Evaluation ```{r poly_eval} # P(x) = 1 + x + x^2 p <- c(1, 1, 1) # x^2 + x + 1 poly_eval(p, c(0, 1, 2)) # returns c(1, 3, 7) # Evaluate at many points x <- seq(-2, 2, length.out = 10) y <- poly_eval(c(1, 0, -1), x) # 1 - x^2 plot(x, y, type = "l", main = "1 - x^2") ``` ### Polynomial Differentiation ```{r polyderiv} # P(x) = x^2 → P'(x) = 2x p <- c(1, 0, 0) # x^2 polyderiv(p, der = 1) # c(2, 0) → 2x # Second derivative: P''(x) = 2 polyderiv(p, der = 2) # c(2) # Higher order derivatives p <- c(1, 2, 3, 4) # 4x^3 + 3x^2 + 2x + 1 polyderiv(p, der = 2) # 24x + 6 ``` ### Polynomial Reduction ```{r reduce_pol} # Remove leading zeros reduce_pol(c(0, 0, 1, 2, 1)) # c(1, 2, 1) → x^2 + 2x + 1 reduce_pol(c(0, 0, 0, 5)) # c(5) → constant polynomial ``` ### Taylor Basis Change ```{r change_polynomial_base_taylor} # Convert (x-1)^2 to expansion around 0 # (x-1)^2 = x^2 - 2x + 1 coeffs_a <- c(1, 0, 0) # (x-1)^2 in basis centered at a=1 change_polynomial_base_taylor(coeffs_a, a = 1, b = 0) # Returns c(1, -2, 1) → x^2 - 2x + 1 ``` ## Piecewise Polynomial (PP) Form The PP form represents a function defined piecewise by polynomials on different intervals. ### Creating a PP Object ```{r makpp} # Create a piecewise polynomial with two intervals # Interval 1: x^2 on [0, 1] # Interval 2: 2x - 1 on [1, 2] coeff <- matrix(c( 1, 0, 0, # x^2 0, 2, -1 # 2x - 1 ), nrow = 2, byrow = TRUE) knots <- c(0, 1, 2) pp <- makpp(coeff, knots) print(pp) # Evaluate the PP form x <- seq(0, 2, length.out = 100) y <- evalpp(pp, x) plot(x, y, type = "l", main = "Piecewise Polynomial") abline(v = knots, col = "red", lty = 2) ``` ### Non-Callable PP Objects ```{r non_callable_pp} # A non-callable PP object is a list with components: # - coeff: matrix of polynomial coefficients # - knot: knot positions # - degree: polynomial degree pp <- makpp(coeff, knots, callable = FALSE) class(pp) # "non_callable_pp" # Access components pp$coeff pp$knot pp$degree ``` ### Callable PP Objects A callable PP object can be evaluated directly like a function. ```{r callable_pp} # Create a callable PP pp_call <- makpp(coeff, knots, callable = TRUE) class(pp_call) # "callable_pp" "function" # Evaluate directly x <- seq(0, 2, length.out = 10) y <- pp_call(x) # Print shows information print(pp_call) # Can be used in plots plot(pp_call, xlim = c(0, 2)) ``` ### displaying polynomial in a human readable form ## Introduction The `show_poly()` function displays a polynomial as a human-readable mathematical equation. It supports both canonical and local bases, and can convert between different expansion points using Taylor's formula. A polynomial $P(x) = \sum_{k=0}^d c_k (x-a)^k$ represented by coefficients `c(c_d, c_{d-1}, ..., c_0)` (decreasing power order) can be displayed in any basis $(x-b)^k$. ## Basic Usage ### Canonical Basis (Default) By default, `show_poly()` displays the polynomial in the canonical basis $1, x, x^2, ...$: ```{r } # P(x) = 3x^3 - 2x^2 + x - 5 p <- c(3, -2, 1, -5) show_poly(p) # "3*x^3 - 2*x^2 + 1*x - 5" ``` ## Working with Local Bases ### Display in Local Basis To display a polynomial in the local basis $(x-a)^k$, use the `b` parameter: ```{r local_basis} # P(x) = 3(x-2)^3 - 2(x-2)^2 + (x-2) - 5 p <- c(3, -2, 1, -5) show_poly(p, b = 2) # "3*(x-2)^3 - 2*(x-2)^2 + 1*(x-2) - 5" ``` ### Converting Between Bases The `a` parameter specifies the basis in which the coefficients are given. The function automatically converts using `change_polynomial_base_taylor()`: ```{r convert_bases} # Coefficients are in basis (x-2)^k, display in canonical basis p <- c(3, -2, 1, -5) # In basis (x-2)^k show_poly(p, a = 2, b = 0) # This expands to: 3x^3 - 18x^2 + 37x - 27 ``` ## Examples with Different Bases ### Positive and Negative Centers ```{r centers} # Center at positive value p <- c(1, 0, 2) # (x-3)^2 + 2 show_poly(p, b = 3) # "(x-3)^2 + 2" # Center at negative value p <- c(1, 0, 2) # (x+2)^2 + 2 show_poly(p, b = -2) # "(x+2)^2 + 2" ``` ## Controlling Output Format ### Number of Digits ```{r digits} p <- c(1/3, -2/7, 1/5) show_poly(p, digits = 2) # "0.33*x^2 - 0.29*x + 0.2" show_poly(p, digits = 4) # "0.3333*x^2 - 0.2857*x + 0.2" ``` ### Negative Coefficients The function handles negative coefficients gracefully: ```{r negative} # All coefficients negative p <- c(-3, -2, -1) # -3x^2 - 2x - 1 show_poly(p) # "-3*x^2 - 2*x - 1" # Mixed coefficients p <- c(3, -2, 1, -5) # 3x^3 - 2x^2 + x - 5 show_poly(p) # "3*x^3 - 2*x^2 + 1*x - 5" ``` ### Working with B-splines When working with B-splines, `show_poly()` helps visualize the polynomial pieces: ```{r bspline_example} # Create a B-spline and display its polynomial form sn <- c(0, 0, 0, 0, 0.3, 0.6, 1, 1, 1, 1) basis <- Bspline_base(sn, degree = 3) basis$coeff <- c(1, -0.5, 2, 0.5, -1, 0.75) # Convert to PP pp <- Bsplinetopp(basis, callable = FALSE) # Display each piece cat("Piece 1: ", show_poly(pp$coeff[1, ], b = 0), "\n") cat("Piece 2: ", show_poly(pp$coeff[2, ], b = 0), "\n") cat("Piece 3: ", show_poly(pp$coeff[3, ], b = 0), "\n") ``` ## Integration with show_pp() The `show_poly()` function is used internally by `show_pp()` and `show_pp()` to display piecewise polynomials: ```{r show_pp_integration} # Create a PP with two intervals knot <- c(0, 0.5, 1) coeff <- matrix(c(1, 2, 0.5, 0, 1, -1), nrow = 2, ncol = 3, byrow = TRUE) pp <- makpp(coeff, knot) # Display in local basis show_pp(pp, local = TRUE, digits = 3) ``` ## Summary | Parameter | Description | Default | |-----------|-------------|---------| | `obj` | Polynomial coefficients (decreasing power order) | Required | | `a` | Base of input coefficients $(x-a)^k$ | 0 | | `b` | Base for output display $(x-b)^k$ | 0 | | `digits` | Number of significant digits | 4 | | `verbose` | Print additional information |FALSE | ### Displaying PP Equations The package provides convenient functions to display PP objects as human-readable mathematical equations. The `show_pp()` and `show_pp()` functions format the polynomial expressions on each interval, with support for both canonical and local bases. ```{r show_pp} # Create a PP with two intervals knot <- c(0, 0.5, 1) coeff <- matrix(c(1, 2, 0.5, 0, 1, -1), nrow = 2, ncol = 3, byrow = TRUE) pp <- makpp(coeff, knot) # Display in canonical basis (1, x, x², ...) show_pp(pp, local = FALSE) # [0.000, 0.500] 0.5x^2 + 2x + 1 # [0.500, 1.000] -1x^2 + 1x + 0 # Display in local basis ((x-a)^i) show_pp(pp, local = TRUE) # [0.000, 0.500] 0.5(x-0)^2 + 2(x-0) + 1 # [0.500, 1.000] -1(x-0.5)^2 + 1(x-0.5) + 0 ``` For spline objects, the conversion to PP is automatic: ```{r} # Create a B-spline and display its polynomial form sn <- c(0, 0, 0, 0, 0.3, 0.6, 1, 1, 1, 1) basis <- Bspline_base(sn, degree = 3) basis$coeff <- c(1, -0.5, 2, 0.5, -1, 0.75) show_pp(basis, local = TRUE, verbose = TRUE) # PP Information: # Degree: 3 # Knots: 0, 0.3, 0.6, 1 # Number of intervals: 3 # # [0.000, 0.300] 1(x-0)^3 + 2.5(x-0)^2 + 1.5(x-0) + 1 # [0.300, 0.600] -4.63(x-0.3)^3 + ... # [0.600, 1.000] ... ``` Customization Options The display functions offer several customization options: local: Control the basis (TRUE for local, FALSE for canonical) digits: Number of significant digits (default: 4) verbose: Display additional PP information ### Callable PP Objects A callable PP object can be evaluated directly like a function. # Create a callable PP pp_call <- makpp(coeff, knots, callable = TRUE) class(pp_call) # "callable_pp" "function" # Evaluate directly x <- seq(0, 2, length.out = 10) y <- pp_call(x) # Print shows information print(pp_call) # Can be used in plots plot(pp_call, xlim = c(0, 2)) ### Extrapolating values When evalpp() is called with x values outside the knot range (x \< knots[1] or x \> knots[length(knots)]), it automatically performs extrapolation by extension using the border polynomial piece (the first or last polynomial in the PP structure). For example, if the knot vector is c(0, 1, 2), points with x \< 0 will be evaluated using the polynomial from the first interval [0, 1], while points with x \> 2 will use the polynomial from the last interval [1, 2]. This is equivalent to "extending" the first and last piece of the spline. A warning message is printed to alert the user when extrapolation occurs. # PP with two intervals: x^2 on [0,1], 2x-1 on [1,2] coeff <- matrix(c(1, 0, 0, 0, 2, -1), nrow = 2, byrow = TRUE) knots <- c(0, 1, 2) pp <- makpp(coeff, knots, callable = TRUE) # Evaluate inside and outside the knot range x_test <- c(-0.5, 0.5, 1.5, 2.5) y_test <- pp(x_test) data.frame(x = x_test, y = y_test) # The values at -0.5 and 2.5 are extrapolated ``` ### Extracting Parameters ```{r get_parameters} # For callable PP objects, use get_parameters() params <- get_parameters(pp_call) print(params$degree) print(params$knot) print(params$coeff) # For non-callable PP, access directly pp=makpp(pp,callable=FALSE,verbose=TRUE) pp$degree pp$knot pp$coeff ``` ## Advanced: PP Form Manipulation ### Extracting a Single Interval Polynomial ```{r extract_interval} # Get polynomial for interval i pp=makpp(pp,callable=FALSE) interval_poly <- pp$coeff[1, ] # First interval print(interval_poly) # Evaluate on a subinterval x_sub <- seq(0, 1, length.out = 20) y_sub <- poly_eval(interval_poly, x_sub) ``` ### Combining PP Objects ```{r combine_pp} # Create two PP objects and combine them # This is useful for constructing complex piecewise functions coeff1 <- matrix(c(1, 0), nrow = 1) coeff2 <- matrix(c(0, 1), nrow = 1) pp1 <- makpp(coeff1, c(0, 1)) pp2 <- makpp(coeff2, c(1, 2)) # Combine coefficients and knots combined_coeff <- rbind(coeff1, coeff2) combined_knots <- c(0, 1, 2) pp_combined <- makpp(combined_coeff, combined_knots) # Evaluate combined function x <- seq(0, 2, length.out = 100) y_combined <- evalpp(pp_combined, x) plot(x, y_combined, type = "l") abline(v = c(0, 1, 2), col = "red", lty = 2) ``` ### Visualizing PP with Basis Functions ```{r visualize_pp} # Create a PP with multiple pieces knots <- seq(0, 1, length.out = 5) degree <- 2 n_pieces <- length(knots) - 1 # Generate random coefficients coeff <- matrix(rnorm(n_pieces * (degree + 1)), nrow = n_pieces) pp <- makpp(coeff, knots, callable = TRUE) # Evaluate and plot x <- seq(0, 1, length.out = 200) y <- pp(x) plot(x, y, type = "l", lwd = 2, main = "Random Piecewise Quadratic") abline(v = knots, col = "red", lty = 2) grid() ``` ## Comparison: PP vs B-spline ```{r pp_vs_bspline} # Create a B-spline and convert to PP form sn <- c(0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1) basis <- Bspline_base(sn, degree = 3) basis$coeff <- runif(basis$n_splines) # Convert to PP pp_from_bspline <- Bsplinetopp(basis) # Both represent the same function x <- seq(0, 1, length.out = 100) y_bspline <- spline_eval(basis, x) y_pp <- evalpp(pp_from_bspline, x) # They should match max(abs(y_bspline - y_pp)) ``` ## Summary