The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
mpoly is a simple collection of tools to help deal
with multivariate polynomials symbolically and functionally in
R. Polynomials are defined with the mp() function:
library("mpoly")
# Registered S3 methods overwritten by 'ggplot2':
# method from
# [.quosures rlang
# c.quosures rlang
# print.quosures rlang
mp("x + y")
# x + y
mp("(x + 4 y)^2 (x - .25)")
# x^3 - 0.25 x^2 + 8 x^2 y - 2 x y + 16 x y^2 - 4 y^2Term orders are available with the reorder function:
(p <- mp("(x + y)^2 (1 + x)"))
# x^3 + x^2 + 2 x^2 y + 2 x y + x y^2 + y^2
reorder(p, varorder = c("y","x"), order = "lex")
# y^2 x + y^2 + 2 y x^2 + 2 y x + x^3 + x^2
reorder(p, varorder = c("x","y"), order = "glex")
# x^3 + 2 x^2 y + x y^2 + x^2 + 2 x y + y^2Vectors of polynomials (mpolyList’s) can be specified in
the same way:
mp(c("(x+y)^2", "z"))
# x^2 + 2 x y + y^2
# zYou can extract pieces of polynoimals using the standard
[ operator, which works on its terms:
p[1]
# x^3
p[1:3]
# x^3 + x^2 + 2 x^2 y
p[-1]
# x^2 + 2 x^2 y + 2 x y + x y^2 + y^2There are also many other functions that can be used to piece apart polynomials, for example the leading term (default lex order):
LT(p)
# x^3
LC(p)
# [1] 1
LM(p)
# x^3You can also extract information about exponents:
exponents(p)
# [[1]]
# x y
# 3 0
#
# [[2]]
# x y
# 2 0
#
# [[3]]
# x y
# 2 1
#
# [[4]]
# x y
# 1 1
#
# [[5]]
# x y
# 1 2
#
# [[6]]
# x y
# 0 2
multideg(p)
# x y
# 3 0
totaldeg(p)
# [1] 3
monomials(p)
# x^3
# x^2
# 2 x^2 y
# 2 x y
# x y^2
# y^2Arithmetic is defined for both polynomials (+,
-, * and ^)…
p1 <- mp("x + y")
p2 <- mp("x - y")
p1 + p2
# 2 x
p1 - p2
# 2 y
p1 * p2
# x^2 - y^2
p1^2
# x^2 + 2 x y + y^2… and vectors of polynomials:
(ps1 <- mp(c("x", "y")))
# x
# y
(ps2 <- mp(c("2 x", "y + z")))
# 2 x
# y + z
ps1 + ps2
# 3 x
# 2 y + z
ps1 - ps2
# -1 x
# -1 z
ps1 * ps2
# 2 x^2
# y^2 + y zYou can compute derivatives easily:
p <- mp("x + x y + x y^2")
deriv(p, "y")
# x + 2 x y
gradient(p)
# y^2 + y + 1
# 2 y x + xYou can turn polynomials and vectors of polynomials into functions
you can evaluate with as.function(). Here’s a basic example
using a single multivariate polynomial:
f <- as.function(mp("x + 2 y")) # makes a function with a vector argument
# f(.) with . = (x, y)
f(c(1,1))
# [1] 3
f <- as.function(mp("x + 2 y"), vector = FALSE) # makes a function with all arguments
# f(x, y)
f(1, 1)
# [1] 3Here’s a basic example with a vector of multivariate polynomials:
(p <- mp(c("x", "2 y")))
# x
# 2 y
f <- as.function(p)
# f(.) with . = (x, y)
f(c(1,1))
# [1] 1 2
f <- as.function(p, vector = FALSE)
# f(x, y)
f(1, 1)
# [1] 1 2Whether you’re working with a single multivariate polynomial or a
vector of them (mpolyList), if it/they are actually
univariate polynomial(s), the resulting function is vectorized. Here’s
an example with a single univariate polynomial.
f <- as.function(mp("x^2"))
# f(.) with . = x
f(1:3)
# [1] 1 4 9
(mat <- matrix(1:4, 2))
# [,1] [,2]
# [1,] 1 3
# [2,] 2 4
f(mat) # it's vectorized properly over arrays
# [,1] [,2]
# [1,] 1 9
# [2,] 4 16Here’s an example with a vector of univariate polynomials:
(p <- mp(c("t", "t^2")))
# t
# t^2
f <- as.function(p)
f(1)
# [1] 1 1
f(1:3)
# [,1] [,2]
# [1,] 1 1
# [2,] 2 4
# [3,] 3 9You can use this to visualize a univariate polynomials like this:
library("tidyverse"); theme_set(theme_classic())f <- as.function(mp("(x-2) x (x+2)"))
# f(.) with . = x
x <- seq(-2.5, 2.5, .1)
qplot(x, f(x), geom = "line")
For multivariate polynomials, it’s a little more complicated:
f <- as.function(mp("x^2 - y^2"))
# f(.) with . = (x, y)
s <- seq(-2.5, 2.5, .1)
df <- expand.grid(x = s, y = s)
df$f <- apply(df, 1, f)
qplot(x, y, data = df, geom = "raster", fill = f)
Using tidyverse-style coding
(install tidyverse packages with
install.packages("tidyverse")), this looks a bit
cleaner:
f <- as.function(mp("x^2 - y^2"), vector = FALSE)
# f(x, y)
seq(-2.5, 2.5, .1) %>%
list("x" = ., "y" = .) %>%
cross_df() %>%
mutate(f = f(x, y)) %>%
ggplot(aes(x, y, fill = f)) +
geom_raster()
Grobner bases are no longer implemented in mpoly; they’re now in m2r.
# polys <- mp(c("t^4 - x", "t^3 - y", "t^2 - z"))
# grobner(polys)Homogenization and dehomogenization:
(p <- mp("x + 2 x y + y - z^3"))
# x + 2 x y + y - z^3
(hp <- homogenize(p))
# x t^2 + 2 x y t + y t^2 - z^3
dehomogenize(hp, "t")
# x + 2 x y + y - z^3
homogeneous_components(p)
# x + y
# 2 x y
# -1 z^3mpoly can make use of many pieces of the
polynom and orthopolynom packages with
as.mpoly() methods. In particular, many special polynomials
are available.
You can construct Chebyshev polynomials as follows:
chebyshev(1)
# Loading required package: polynom
#
# Attaching package: 'polynom'
# The following object is masked from 'package:mpoly':
#
# LCM
# x
chebyshev(2)
# -1 + 2 x^2
chebyshev(0:5)
# 1
# x
# 2 x^2 - 1
# 4 x^3 - 3 x
# 8 x^4 - 8 x^2 + 1
# 16 x^5 - 20 x^3 + 5 xAnd you can visualize them:
s <- seq(-1, 1, length.out = 201); N <- 5
(chebPolys <- chebyshev(0:N))
# 1
# x
# 2 x^2 - 1
# 4 x^3 - 3 x
# 8 x^4 - 8 x^2 + 1
# 16 x^5 - 20 x^3 + 5 x
df <- as.function(chebPolys)(s) %>% cbind(s, .) %>% as.data.frame()
names(df) <- c("x", paste0("T_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree)
s <- seq(-1, 1, length.out = 201); N <- 5
(jacPolys <- jacobi(0:N, 2, 2))
# 1
# 5 x
# 17.5 x^2 - 2.5
# 52.5 x^3 - 17.5 x
# 144.375 x^4 - 78.75 x^2 + 4.375
# 375.375 x^5 - 288.75 x^3 + 39.375 x
df <- as.function(jacPolys)(s) %>% cbind(s, .) %>% as.data.frame
names(df) <- c("x", paste0("P_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree) +
coord_cartesian(ylim = c(-25, 25))
s <- seq(-1, 1, length.out = 201); N <- 5
(legPolys <- legendre(0:N))
# 1
# x
# 1.5 x^2 - 0.5
# 2.5 x^3 - 1.5 x
# 4.375 x^4 - 3.75 x^2 + 0.375
# 7.875 x^5 - 8.75 x^3 + 1.875 x
df <- as.function(legPolys)(s) %>% cbind(s, .) %>% as.data.frame
names(df) <- c("x", paste0("P_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree)
s <- seq(-3, 3, length.out = 201); N <- 5
(hermPolys <- hermite(0:N))
# 1
# x
# x^2 - 1
# x^3 - 3 x
# x^4 - 6 x^2 + 3
# x^5 - 10 x^3 + 15 x
df <- as.function(hermPolys)(s) %>% cbind(s, .) %>% as.data.frame
names(df) <- c("x", paste0("He_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree)
s <- seq(-5, 20, length.out = 201); N <- 5
(lagPolys <- laguerre(0:N))
# 1
# -1 x + 1
# 0.5 x^2 - 2 x + 1
# -0.1666667 x^3 + 1.5 x^2 - 3 x + 1
# 0.04166667 x^4 - 0.6666667 x^3 + 3 x^2 - 4 x + 1
# -0.008333333 x^5 + 0.2083333 x^4 - 1.666667 x^3 + 5 x^2 - 5 x + 1
df <- as.function(lagPolys)(s) %>% cbind(s, .) %>% as.data.frame
names(df) <- c("x", paste0("L_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree) +
coord_cartesian(ylim = c(-25, 25))
Bernstein
polynomials are not in polynom or
orthopolynom but are available in
mpoly with bernstein():
bernstein(0:4, 4)
# x^4 - 4 x^3 + 6 x^2 - 4 x + 1
# -4 x^4 + 12 x^3 - 12 x^2 + 4 x
# 6 x^4 - 12 x^3 + 6 x^2
# -4 x^4 + 4 x^3
# x^4
s <- seq(0, 1, length.out = 101)
N <- 5 # number of bernstein polynomials to plot
(bernPolys <- bernstein(0:N, N))
# -1 x^5 + 5 x^4 - 10 x^3 + 10 x^2 - 5 x + 1
# 5 x^5 - 20 x^4 + 30 x^3 - 20 x^2 + 5 x
# -10 x^5 + 30 x^4 - 30 x^3 + 10 x^2
# 10 x^5 - 20 x^4 + 10 x^3
# -5 x^5 + 5 x^4
# x^5
df <- as.function(bernPolys)(s) %>% cbind(s, .) %>% as.data.frame
names(df) <- c("x", paste0("B_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree)
You can use the bernstein_approx() function to compute
the Bernstein polynomial approximation to a function. Here’s an
approximation to the standard normal density:
p <- bernstein_approx(dnorm, 15, -1.25, 1.25)
round(p, 4)
# -0.1624 x^2 + 0.0262 x^4 - 0.002 x^6 + 0.0001 x^8 + 0.3796
x <- seq(-3, 3, length.out = 101)
df <- data.frame(
x = rep(x, 2),
y = c(dnorm(x), as.function(p)(x)),
which = rep(c("actual", "approx"), each = 101)
)
# f(.) with . = x
qplot(x, y, data = df, geom = "path", color = which)
You can construct Bezier polynomials
for a given collection of points with bezier():
points <- data.frame(x = c(-1,-2,2,1), y = c(0,1,1,0))
(bezPolys <- bezier(points))
# -10 t^3 + 15 t^2 - 3 t - 1
# -3 t^2 + 3 tAnd viewing them is just as easy:
df <- as.function(bezPolys)(s) %>% as.data.frame
ggplot(aes(x = x, y = y), data = df) +
geom_point(data = points, color = "red", size = 4) +
geom_path(data = points, color = "red", linetype = 2) +
geom_path(size = 2)
Weighting is available also:
points <- data.frame(x = c(1,-2,2,-1), y = c(0,1,1,0))
(bezPolys <- bezier(points))
# -14 t^3 + 21 t^2 - 9 t + 1
# -3 t^2 + 3 t
df <- as.function(bezPolys, weights = c(1,5,5,1))(s) %>% as.data.frame
ggplot(aes(x = x, y = y), data = df) +
geom_point(data = points, color = "red", size = 4) +
geom_path(data = points, color = "red", linetype = 2) +
geom_path(size = 2)
To make the evaluation of the Bezier polynomials stable,
as.function() has a special method for Bezier polynomials
that makes use of de
Casteljau’s algorithm. This allows bezier() to be used
as a smoother:
s <- seq(0, 1, length.out = 201)
df <- as.function(bezier(cars))(s) %>% as.data.frame
qplot(speed, dist, data = cars) +
geom_path(data = df, color = "red")
I’m starting to put in methods for some other R functions:
set.seed(1)
n <- 101
df <- data.frame(x = seq(-5, 5, length.out = n))
df$y <- with(df, -x^2 + 2*x - 3 + rnorm(n, 0, 2))
mod <- lm(y ~ x + I(x^2), data = df)
(p <- mod %>% as.mpoly %>% round)
# 1.983 x - 1.01 x^2 - 2.709
qplot(x, y, data = df) +
stat_function(fun = as.function(p), colour = 'red')
# f(.) with . = x
s <- seq(-5, 5, length.out = n)
df <- expand.grid(x = s, y = s) %>%
mutate(z = x^2 - y^2 + 3*x*y + rnorm(n^2, 0, 3))
(mod <- lm(z ~ poly(x, y, degree = 2, raw = TRUE), data = df))
#
# Call:
# lm(formula = z ~ poly(x, y, degree = 2, raw = TRUE), data = df)
#
# Coefficients:
# (Intercept)
# -0.070512
# poly(x, y, degree = 2, raw = TRUE)1.0
# -0.004841
# poly(x, y, degree = 2, raw = TRUE)2.0
# 1.005307
# poly(x, y, degree = 2, raw = TRUE)0.1
# 0.001334
# poly(x, y, degree = 2, raw = TRUE)1.1
# 3.003755
# poly(x, y, degree = 2, raw = TRUE)0.2
# -0.999536
as.mpoly(mod)
# -0.004840798 x + 1.005307 x^2 + 0.001334122 y + 3.003755 x y - 0.9995356 y^2 - 0.07051218From CRAN: install.packages("mpoly")
From Github (dev version):
# install.packages("devtools")
devtools::install_github("dkahle/mpoly")This material is based upon work partially supported by the National Science Foundation under Grant No. 1622449.
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.