R

Ryacas version 1.1.0 publised in Journal of Open Source Software and released to CRAN

It is with great pleasure that I can announce that Ryacas version 1.1.0 has now been accepted into Journal of Open Source Software and same version released to CRAN. (The source code is available at https://github.com/mikldk/ryacas/.) I already wrote about Ryacas many times before. I will refer you to the “Getting started” and “The high-level (symbol) interface” vignettes or one of the others available at the CRAN page or the package’s website.

Variance in reproductive success (VRS) in forensic genetics lineage markers

Back in 2017, David Balding and I published the paper “How convincing is a matching Y-chromosome profile?”. One of the key parameters of the simulation model was the variance in reproductive success (VRS). Here I will discuss and demonstrate this parameter. First note for intuition that in a Wright–Fisher model, all individuals have the same probability of becoming a father, or of having reproductive success. So the VRS here is 0.

Wolfe conditions for deciding step length in inexact line search: An example with the Rosenbrock function

In inexact line search (a numerical optimisation technique) the step length (or learning rate) must be decided. In connection to that the Wolfe conditions are central. Here I will give an example showing why they are useful. More on this topic can be read elsewhere, e.g. in the book of Nocedal and Wright (2006), “Numerical Optimization”, Springer. Wolfe conditions The Wolfe conditions consists of the sufficient decrease condition (SDC) and curvature condition (CC):

Ryacas version 1.0.0 released!

It is with great pleasure that I can announce that Ryacas version 1.0.0 is now released to CRAN (https://cran.r-project.org/package=Ryacas). I wish to thank all co-authors: Rob Goedman, Gabor Grothendieck, Søren Højsgaard, Grzegorz Mazur, Ayal Pinkus. It means that you can install the package by (possible after binaries have been built): install.packages("Ryacas") Followed by: library(Ryacas) (The source code is available at https://github.com/mikldk/ryacas/.) Now you have the yacas computer algebra system fully available!

Approximating small probabilities using importance sampling

Update Oct 14, 2019: Michael Höhle caught a mistake and notified me on Twitter. Thanks! The problem is that I used \(\text{Unif}(-10, 10)\) as importance distribution; this does not have infinite support as the target has. This is required, see e.g. Art B. Owen (2013), “Monte Carlo theory, methods and examples”. I have now updated the post to use a normal distribution instead. Box plots are often used. They are not always the best visualisation (e.

How much pizza and how much frozen yogurt? ...with Gröbner bases

In a recent blog post I tried to get yacas to solve a system of polynomial equations. Unfortunately it could not do that, so I solved it numerically instead. Now it is possible – together with many other systems of polynomial equations thanks to fixing a small error in yacas. It has now been fixed, also in Ryacas (development version), so hurry up and update Ryacas to the latest version 0.

Prediction intervals for Generalized Additive Models (GAMs)

Finding prediction intervals (for future observations) is something different than finding confidence intervals (for unknown population parameters). Here, I demonstrate one approach to doing so. First we load the library and simulate some data: library(mgcv) set.seed(1) dat <- gamSim(eg = 1, n = 400, dist = "normal", scale = 2) ## Gu & Wahba 4 term additive model The simulated in dat contains the “truth” in the f variables:

The cost of evaluating an expression instead of using the language directly

In a recent blog post I used something like this (for use in a call to optim()): obj_fun_expr <- expression((x^4 + 4 * x^2 * y^2 - 12 * x^2 * a + 144 * x^2 - 48 * x * y * a + 144 * x * y - 4320 * x + 36 * y^2 - 2160 * y + 180 * a^2 + 32400)/16) f_expr <- function(par) { x <- par[1] y <- par[2] a <- par[3] val <- eval(obj_fun_expr, list(x = x, y = y, a = a)) return(val) } I have been wondering what the cost of eval(expr, .

Correlation is not transitive, in general at least: A simulation approach

Let \(\rho_{XY}\) be the correlation between the stochastic variables \(X\) and \(Y\) and similarly for \(\rho_{XZ}\) and \(\rho_{YZ}\). If we know two of these, can we say anything about the third? In a recent blog post I dealt with the problem mathematically and I used the concept of a partial correlation coefficient. Here I will take a simulation approach. First z is simulated. Then x and y is simulated based on z in a regression context with a slope between \(-1\) and \(1\).

Automatic 'testthat' test skeletons with new R package 'roxytest' extending 'roxygen2'

It is important to test software. One approach is unit-testing, and for R packages this can e.g. be done using testthat. It is also important to document software. For R packages roxygen2 is really helpful: It enables you to write documentation in the code file in the R/ folder where the function is implemented. And then roxygen2 takes care of handling the Rd files in the man/ folder. I have made a new R package that combines these approaches: roxytest.