kloppenborg.ca

Explaining A- and B-Basis Values

2022-11-21T00:00:00-05:00

I’ve been thinking about how to explain an A- or B-Basis value to people without much statistical knowledge. These are the names used in aircraft certification for the lower tolerance bounds on material strength. The definition used for transport category aircraft is given in 14 CFR 25.613 (the same definition is used for other categories of aircraft too). This definition is precise, but not easy to understand.

… (b) Material design values must be chosen to minimize the probability of structural failures due to material variability. … compliance must be shown by selecting material design values which assure material strength with the following probability:

(1) Where applied loads are eventually distributed through a single member within an assembly, the failure of which would result in loss of structural integrity of the component, 99 percent probability with 95 percent confidence.

(2) For redundant structure, in which the failure of individual elements would result in applied loads being safely distributed to other load carrying members, 90 percent probability with 95 percent confidence. …

Another way of stating this definition is the lower 95% confidence bound on the 1-st or 10-th percentile of the population, respectively. But, describing it that way doesn’t help to explain the concept to a person who’s not well versed in statistics.

The Explanation

There’s some random variation in all material properties. Some pieces of any material will be a little bit different than other pieces of the same material. To account for this variation in the material properties when we design aircraft structure, we design it so that there’s at least a 90% chance that redundant structure is stronger than it needs to be (or 99% of non-redundant structure).

When we test a material property, we get a sample. This sample is not a perfect representation of the material property. A good analogy is that a sample is like a low-resolution photo: it gives us an idea of what we’re seeing, but we don’t get all the detail. We can get a better idea of what we’re seeing by taking a higher resolution photo: this is akin to testing more and getting a larger sample size.

We choose a statistical distribution that fits the data, then find the tenth (or first) percentile of that distribution. But since we only have a sample of the material property (a “low-resolution photo”, in the analogy), we’re not sure if the distribution that we chose is correct. To account for that uncertainty, we try out many possible distributions for the material property and determine how likely each is to be true based on the sample (the data). Distributions that look a lot like the data are highly likely; distributions that look different than the data are less likely, but depending on how “low-resolution” our data is, they could be correct. For each of these possible distributions, we find the 10th percentile (for B-Basis; it would be the 1st percentile for A-Basis). Next, we weight each of those individual 10th percentiles based on the likelihood that the corresponding distribution is true, and we find a lower bound where 95% of those weighted 10th percentiles are above that lower bound.

Or in graphical form:

I hope that this explanation helps explain this complicated topic. If you think you have a better explanation, please connect with me on LinkedIn and message me there.

Developing the Graph

Let’s look at how I developed this graph. The graph was developed using the R language. As for most R code, we start by loading the required packages.

library(tidyverse)
library(cmstatr)
library(stats4)

In this example, we’ll use some of the sample data that comes with the cmstatr package. We’ll be using the room-temperature warp-tension example data.

dat <- carbon.fabric.2 %>% 
  filter(test == "WT" & condition == "RTD")

Let’s start by plotting this data. We’re plotting 1-D data, so we only need one axis. But in order to make sure that none of the data points overlap, we’ll add some “jitter” (random vertical position). We’ll also hide the vertical axis since this axis is not meaningful.

dat %>% 
  ggplot(aes(x = strength)) +
  geom_jitter(aes(y = 0.01), height = 0.005) +
  scale_y_continuous(name = NULL,
                    breaks = NULL)

We can fit a normal distribution to this data. The sample mean and standard deviation are point-estimates of the mean and standard deviation of the distribution. We’ll use those point-estimates and draw the PDF superimposed over the data, assuming that the distribution is normal. We can also add the 10-th percentile of this distribution to the plot.

dat %>% 
  ggplot(aes(x = strength)) +
  geom_jitter(aes(y = 0.01), height = 0.005, color = "magenta") +
  stat_function(fun = function(x) dnorm(x, mean(dat$strength),
                                        sd(dat$strength))) +
  geom_vline(xintercept = qnorm(0.1, mean(dat$strength), sd(dat$strength)),
             color = "blue") +
  scale_y_continuous(name = NULL,
                    breaks = NULL)

But, the distribution that we’ve drawn is just a point-estimate from the data. There is uncertainty in our estimate. Based on the data, we’ve concluded that this estimate is the most likely, but we shouldn’t be surprised if the true population distribution is a bit different. This point-estimate is actually the Maximum Likelihood Estimate (MLE), based on this particular data. We can calculate the likelihood of various potential estimates of distribution (or rather, the parameters of the distribution) using the following equation:

$$ L\left(\mu, \sigma\right) = \prod_{i=1}^{n} f\left(X_i;\,\mu, \sigma\right) $$

Here $X_i$ are the data and there are two parameters for a normal distribution are $\mu$ and $\sigma$. The function $f()$ is the probability density function.

It turns out that computers have trouble multiplying a bunch of small numbers together and coming up with an accurate result. We can avoid this problem by using a log transform:

$$ \mathcal{L}\left(\mu, \sigma\right) = \sum_{i=1}^{n} \log f\left(X_i;\,\mu, \sigma\right) $$

Implementing this in R:

log_likelihood <- function(x, mu, sigma) {
  sum(dnorm(x, mu, sigma, log = TRUE))
}

To make sure that this function works, we can find the log likelihood of our MLE of the parameters. The actual numerical value of the likelihood doesn’t mean very much to us, but we’ll be interested in the distribution of the likelihoods as we change the parameters.

log_likelihood(dat$strength, mean(dat$strength), sd(dat$strength))

## [1] -92.55627

It’s going to make our life a lot easier if we can work with a single parameter instead of two ($\mu$ and $\sigma$). We’ll treat $\sigma$ as a nuisance parameter and find the value of $\sigma$ that produces the greatest likelihood for any given values of $\mu$. To avoid working with very tiny numbers, we’ll calculate the relative likelihood (the likelihood divided by the maximum likelihood). We can do this in R as follows:

rel_likelihood_mu <- function(x, mu) {
  ll_hat <- log_likelihood(x, mean(x), sd(x))

  opt <- optimize(
    function(sigma) exp(log_likelihood(x, mu, sigma) - ll_hat),
    lower = 0,
    upper = 20 * sd(x),  # pick an upper bound that's big
    maximum = TRUE
  )
  # We'll return a list of the sigma and the relative likelihood:
  list(
    sigma = opt$maximum,
    rel_likelihood = opt$objective
  )
}

We can also do the same thing to calculate the relative likelihood of a a particular 10th percentile ($x_p$). We use the transformation $\mu = x_p - \sigma \Phi(0.1)$.

rel_likelihood_xp <- function(x, xp) {
  ll_hat <- log_likelihood(x, mean(x), sd(x))

  opt <- optimize(
    function(sigma) exp(log_likelihood(x, xp - sigma * qnorm(0.1), sigma) - ll_hat),
    lower = 0,
    upper = 20 * sd(x),  # pick an upper bound that's big
    maximum = TRUE
  )
  # We'll return a list of the sigma and the relative likelihood:
  list(
    sigma = opt$maximum,
    rel_likelihood = opt$objective
  )
}

Now, we can draw our same plot again, but this time, we’ll draw a bunch of potential distributions (and 10-th percentiles), coloring them according to their likelihood:

p <- dat %>% 
  ggplot(aes(x = strength))

walk(
  seq(from = 0.95 * mean(dat$strength),
      to = 1.05 * mean(dat$strength),
      length.out = 55),
  function(mu) {
    rl <- rel_likelihood_mu(dat$strength, mu)
    p <<- p + stat_function(aes(alpha = rl$rel_likelihood),
                            fun = function(x) dnorm(x, mu, rl$sigma))
  }
)

xp_dist <- imap_dfr(
  seq(from = 0.9 * mean(dat$strength),
      to = 0.98 * mean(dat$strength),
      length.out = 55),
  function(xp, ii) {
    rl <- rel_likelihood_xp(dat$strength, xp)
    data.frame(xp = xp, rl = rl$rel_likelihood)
  })

p +
  geom_vline(aes(xintercept = xp, alpha = rl), data = xp_dist, color = "blue") +
  geom_jitter(aes(y = 0.01), height = 0.005, color = "magenta") +
  scale_y_continuous(name = NULL,
                    breaks = NULL)

Next, we can add a plot of the distribution of the 10th percentiles. We’ll also plot the B-Basis, as calculated by the package cmstatr.

p <- dat %>% 
  mutate(f = 1) %>% 
  ggplot(aes(x = strength))

walk(
  seq(from = 0.95 * mean(dat$strength),
      to = 1.05 * mean(dat$strength),
      length.out = 55),
  function(mu) {
    rl <- rel_likelihood_mu(dat$strength, mu)
    p <<- p + stat_function(aes(alpha = rl$rel_likelihood),
                            fun = function(x) dnorm(x, mu, rl$sigma))
  }
)

xp_dist <- imap_dfr(
  seq(from = 0.9 * mean(dat$strength),
      to = 0.98 * mean(dat$strength),
      length.out = 55),
  function(xp, ii) {
    rl <- rel_likelihood_xp(dat$strength, xp)
    data.frame(xp = xp, rl = rl$rel_likelihood)
  })

p +
  geom_vline(aes(xintercept = xp, alpha = rl), data = xp_dist, color = "blue") +
  geom_vline(
    aes(xintercept = basis_normal(dat, strength)$basis),
    color = "red",
    data = xp_dist %>% mutate(f = 2),
    inherit.aes = FALSE
  ) +
  geom_line(aes(x = xp, y = rl),
            data = xp_dist %>% 
              mutate(f = 2),
            inherit.aes = FALSE,
            color = "black") +
  geom_jitter(aes(y = 0.01), height = 0.005, color = "magenta") +
  facet_grid(f ~ .,
             scales = "free_y") +
  scale_y_continuous(name = NULL,
                    breaks = NULL) +
  theme(strip.text = element_blank())

The astute reader might recognize the lower curve as a non-central t-distribution. Since it’s a relative likelihood and not a probability, the vertical scale (which is hidden) won’t match the non-central t-distribution, but it’s the same shape. Just for fun, we can plot the lower curve shown above and a non-central t-distribution:

bind_rows(
  xp_dist %>% 
    mutate(f = "Relative Likelihood"),
  xp_dist %>% 
    mutate(n = length(dat$strength),
           rl = dt(sqrt(n) * (mean(dat$strength) - xp) / sd(dat$strength), df = n - 1, ncp = qnorm(0.9) * sqrt(n)),
           f = "t-Distribution") %>% 
    select(-c(n))
) %>% 
  ggplot(aes(x = xp, y = rl)) +
  geom_line() +
  facet_grid(f ~ ., scales = "free_y") +
  ylab("")

Turning our attention back to creating the graph for this blog post, we’ll improve the aesthetics of the graph and also add the annotations:

p <- dat %>% 
  mutate(f = 1) %>% 
  ggplot(aes(x = strength))

walk(
  seq(from = 0.95 * mean(dat$strength),
      to = 1.05 * mean(dat$strength),
      length.out = 55),
  function(mu) {
    rl <- rel_likelihood_mu(dat$strength, mu)
    p <<- p + stat_function(aes(alpha = rl$rel_likelihood),
                            fun = function(x) dnorm(x, mu, rl$sigma))
  }
)

xp_dist <- imap_dfr(
  seq(from = 0.9 * mean(dat$strength),
      to = 0.98 * mean(dat$strength),
      length.out = 55),
  function(xp, ii) {
    rl <- rel_likelihood_xp(dat$strength, xp)
    data.frame(xp = xp, rl = rl$rel_likelihood)
  })

p +
  geom_vline(aes(xintercept = xp, alpha = rl), data = xp_dist, color = "blue") +
  geom_vline(
    aes(xintercept = basis_normal(dat, strength)$basis),
    color = "red",
    data = xp_dist %>% mutate(f = 2),
    inherit.aes = FALSE
  ) +
  geom_line(aes(x = xp, y = rl),
            data = xp_dist %>% 
              mutate(f = 2),
            inherit.aes = FALSE,
            color = "black") +
  geom_jitter(aes(y = 0.01), height = 0.005, color = "magenta") +
  facet_grid(f ~ .,
             scales = "free_y") +
  theme_bw() +
  guides(alpha = guide_none()) +
  scale_y_continuous(name = NULL,
                    breaks = NULL) +
  theme(strip.text = element_blank()) +
  geom_text(
    aes(x = x, y = y, f = f, label = label),
    data = tribble(
      ~x, ~y, ~f, ~label,
      140, 0.025, 1, "(1) The data tells us which\ndistributions are most likely.",
      140, 0.00, 1,
      "...but we don't know the true distribution.",
      140, 0.7, 2,
      "(2) The data also tells us which\n10th percentiles are likely.",
      123, 0.6, 2,
      "(3) Considering the\nlikelihood of all the\npossible 10th\npercentiles, there is\n95% confidence that\nthe true values is above\nthe the B-Basis."
    ),
    color = c("black", "black", "blue", "red")
  ) +
  xlim(c(120, 150))

And that’s the graph at the beginning of this post.

Blogging with Quarto

2022-07-18T00:00:00-04:00

I’ve recently started using Quarto, which is a new open source project backed by RStudio. Quarto is a system for producing reports, presentations, books and blog posts. It takes text formatted with markdown and code written in Python or R and produces PDFs, HTML or several other formats that contain the formatted text, the code (optionally) and the outputs from that code. In a lot of ways, Quarto is like R Markdown or Jupyter Notebooks. Quarto uses Pandoc to do actually convert document formats, and Quarto actually works quite well with version control software like git, unlike Jupyter Notebooks.

I’ve written about using R Markdown or Jupyter for reproducibility in engineering reports, and I’ve written about creating custom document templates for reports written using Pandoc. Much of what I’ve written in those posts should be applicable to Quarto.

To date, I’ve written two posts on this blog using Quarto: Violin Bow Stiffness and Shear of Adhesive Bonded Joints. This post describes some of my experiences using Quarto. Overall, I’ve been quite happy with it.

Blogging with Quarto and Pelican

Quarto has built-in support for blogging with Hugo. However, this blog is using Pelican, not Hugo. A few tweaks are needed.

Since the Quarto posts are written in Markdown, they have YAML headers. Here is the header for one of my blog posts:

---
title: Shear of Adhesive Bonded Joints
date: 2022-06-25
format: commonmark_x
keep-yaml: true
Tags: Engineering, Python, Adhesive Bonding
Category: Posts
filters:
    - attach-filter.lua
---

Let’s go through the lines in this header one at a time:

title: Self-explanatory — the title of the post
date: Also self-explanatory — the date of the post
format: There are several output formats for Quarto. I’ve found that commonmark_x works the best for Pelican. This output format produces a .md file in a format that (mostly) works with Pelican.
keep-yaml: Setting this option to true tells Quarto to copy the present YAML header to the output .md file.
Tags: This is an option used by Pelican. Since we’ve set keep-yaml = true, this gets copied to the .md file that Pelican will process.
Category: Another option used by Pelican.
filters: We’ll talk about this next.

Lua Filters

Pandoc uses something called a filter to alter the output. These filters are written in a language called Lua. In order for Pelican to include an image, the filename of the image needs to start with {attach}. This tells Pelican to include the image file in the website output.

The following filter edits each image element when it’s being processed. The name of the filter (Image) means that it applies to images. This filter concatenates the string {attach} with the src attribute of the image and stores the result in the src attribute of the resulting element.

function Image (elem)
    elem.src = "{attach}" .. elem.src
    return elem
end

Similarly, for links, Pelican requires that links to internal files on the blog start with {filename}. External links are used as-is. To do this, I use the following filter, which applies to Link elements. It checks if the target of the link starts with http. If so, it uses the link as-is. Otherwise, I assume that the link is an internal link, so the filter pre-pends the target with {filename}.

function Link (elem)
    if( string.find(elem.target, "http", 0) )
    then
        return elem
    else
        elem.target = "{filename}" .. elem.target
        return elem
    end
end

I’ve created a file called attach-filter.lua containing both of the filters above. The filters line in the YAML header tells Quarto to use these filters when processing the file.

Shear of Adhesive Bonded Joints

2022-06-25T00:00:00-04:00

There are a lot of misconceptions about bonded joints. One of the misconceptions that I’ve seen most often is that people think that the average shear stress in a lap joint is predictive of the strength. This same misconception is usually phrased as either:

Doubling the overlap length of a lap joint doubles the strength (wrong!)
Calculate $P/A$ for the joint and make sure that the value is less than the lap shear strength on the adhesive data-sheet (wrong!)

In this post, I’m going to explain why these statements are incorrect. I’m going to try to give you a understanding of how load transfer works in an adhesive joint, and I’m going to share some Python code that produces a first approximation of the stress distribution.

For simplicity, we’re going to ignore the effects of peel. Peel is the tendency for the ends of a lap joint to separate. This can cause the joint to fail in some cases, but considering the effects of peel complicates the analysis of the joint — since the purpose of this post is to give a basic understanding of the mechanics of the joint, I’m going to ignore this complicating factor.

The Joint

In this post, I’m going to focus on a simple lap joint. In this type of joint, two adherends overlap each other by a certain amount and there is adhesive connecting the two adherends over the area in which they overlap. These two adherends are then pulled apart. In this post, we’re going to assume that the two adherends are homogeneous isotropic materials (for example, sheet metal) and are uniform thickness. This joint is shown in the figure below. In the top part of the figure, we see the unstressed joint, and in the bottom, we see the joint under load.

Obviously, the deformation of the joint is exaggerated, but it allows us to see what’s happening.

First, let’s look at the lower adherend. We see that the left edge of the adherend is built-in (i.e. it can’t move). When load is applied, the left portion of the lower adherend stretches a lot because it is carrying the entirety of the reaction load.

As we move our gaze further to the right, but still focusing on the lower adherend, we see that the further right we go, the less the adherend is stretching. This is because the adhesive is transferring load along the length of the joint. When we look at the right portion of the lower adherend, it’ hardly stretched at all. Sure, it moved because the rest of the adherend has stretched, but the right part of the lower adherend has hardly stretched at all since it’s carrying no load.

If the two adherends are the same thickness, then symmetry will tell us that the upper adherend behaves in the same way — but now is the right end of the upper adherend that stretches a lot and the left end that doesn’t stretch.

Now, let’s turn our attention to the adhesive. Shear strain can be though of as an angle. At the very left edge of adhesive, the shear strain is quite large, in the middle of the adhesive, the shear strain is moderate and at the right edge of the adhesive, the shear strain is quite large again. The relationship between shear stress and shear strain of an adhesive is not linear, but nonetheless, a large strain produces a large stress and a small strain produces a small stress. So, the shear stress distribution in this joint is “U” shaped — there’s a lot of stress at the ends and a smaller stress in the middle.

This “U” shaped shear stress distribution should be the first clue about why using the average shear stress in the joint to predict failure might not be the best idea.

A Linear Model of the Joint

The actual shear stress-shear strain relationship for most adhesives is non-linear, but we’ll start our analysis of a lap joint by making the assumption that the adhesive is linear-elastic.

Let’s start with defining the variables that we’ll need. The variables are shown in the following figure.

A force balance on the two adherends gives us:

$$ \frac{dN_1}{dx} - \tau = 0 $$

$$ \frac{dN_2}{dx} + \tau = 0 $$

We can find the deformation of the two adherends as follows:

$$ \frac{du_1}{dx} = \frac{N_1}{E_1^\prime t_1} $$

$$ \frac{du_2}{dx} = \frac{N_2}{E_2^\prime t_2} $$

Where $E_1^\prime$ and $E_2^\prime$ are the adherend plane-strain elastic moduli.

The shear strain of the adhesive layer is given by:

$$ \gamma = \frac{1}{t_A} \left(u_1 - u_2\right) $$

We can differentiate this with respect to $x$ and then substitute in the previous equations to get:

$$ \frac{d\gamma}{dx} = \frac{1}{t_A}\left( \frac{du_1}{dx} - \frac{du_2}{dx} \right) \ {} = \frac{1}{t_A}\left( \frac{N_1}{E_1^\prime t_1} - \frac{N_2}{E_2^\prime t_2} \right) $$

We can then differentiate this again with respect to $x$ and substituting in the first equations, we get:

$$ \frac{d^2\gamma}{dx^2} = \frac{1}{t_A}\left( \frac{dN_1}{dx}\frac{1}{E_1^\prime t_1} - \frac{dN_2}{dx}\frac{1}{E_2^\prime t_2} \right) \ {} = \tau\frac{1}{t_A}\left( \frac{1}{E_1^\prime t_1} + \frac{1}{E_2^\prime t_2} \right) $$

Remember that for now, we’re assuming that the adhesive is linear-elastic. Thus:

$$ \tau = G_A \gamma $$

We can solve the second-order differential equation above, but we need two boundary conditions. The boundary conditions that we choose are the loads at the ends of the adherends. At the left end ($x=0$), the unit load on the lower adherend ($N_2$) must be equal to the applied load ($P$) divided by the width ($w$) and the load on the upper adherend ($N_1$) must be zero. The opposite is true at the other end ($x=L$). Thus:

$$ \left.N_1\right|{x=0} = 0 \ \left.N_2\right| = P/w $$

$$ \left.N_1\right|{x=L} = P/w \ \left.N_2\right| = 0 $$

We can plug these into the equation for $\frac{d\gamma}{dx}$ at the two ends of the joint and get the following boundary conditions that we will enforce for the solution.

$$ \left.\frac{d\gamma}{dx}\right|_{x=0} = \frac{1}{t_A}\left( \frac{-P / w}{E_2^\prime t_2} \right) $$

$$ \left.\frac{d\gamma}{dx}\right|_{x=L} = \frac{1}{t_A}\left( \frac{P / w}{E_1^\prime t_1} \right) $$

There is a closed-form solution to this boundary value problem, which we could find, but I think it’s more instructive to just find a numerical solution — plus it’s easier to extend the numerical solution to the case where the adhesive is non-linear. In order to find the numerical solution, we’re going to use the Python package scipy, which includes the function solve_bvp() for solving boundary-value problems. We’ll start by importing the packages that we’ll use.

import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt

Next, we’ll set the parameters for our solution. These include the elastic moduli, thicknesses, overlap length and load.

E1 = 10.5e6 / (1 - 0.33**2)
E2 = 10.5e6 / (1 - 0.33**2)
t1 = 0.063
t2 = 0.063
Ga = 65500
ta = 0.005

L = 0.5
w = 1.

P = 2700

The function solve_bvp requires two arguments: (i) a function that returns the derivatives of the variables, and (ii) a function that returns the residuals for the boundary conditions. We’re going to reduce the second-order differential equation to a system of two first-order differential equations by defining $y$ as follows. Based on this definition, we can implement the two functions required by solve_bvp.

$$ y = \left[ \begin{matrix} \frac{d\tau}{dx} & \tau \end{matrix} \right]^T $$

def func1(x, y):
  D = np.matrix([
    [0, Ga / ta * (1. / (E1 * t1) + 1. / (E2 * t2))],
    [1, 0]
  ])
  return D @ y

def bc1(ya, yb):
  return np.array([
    ya[0] - 1. / ta * (-P / w / (E2 * t2)),
    yb[0] - 1. / ta * (P / w / (E1 * t1))
  ])

res1 = scipy.integrate.solve_bvp(
  func1,
  bc1,
  x = np.linspace(0, L, num=50),
  y = np.zeros((2, 50))
)

The variable res1 now contains the solution to our differential equation. We can plot the shear strain ($\gamma$) over the length of the joint as follows:

plt.plot(res1.x, res1.y[1,:])
plt.title("Linear Elastic Adhesive")
plt.ylabel("Shear Strain, $\\gamma$")
plt.xlabel("$x$")
plt.show()

Because we’re assuming that the adhesive is linear-elastic, we can find the shear stress by simply multiplying the elastic modulus $G_A$ by the shear strain. The shear stress in the adhesive over the length of the joint is thus:

plt.plot(res1.x, Ga * res1.y[1,:])
plt.title("Linear Elastic Adhesive")
plt.ylabel("Shear Stress, $\\tau$")
plt.xlabel("$x$")
plt.show()

Stress-Strain Curve

The shear stress-strain curve for most adhesives is linear at low strain, but highly nonlinear above a certain value of strain. It’s common to idealize the stress-strain curve for an adhesive as elastic-perfectly plastic. The important parameters for the adhesive stress-strain curve are the initial shear modulus ($G_A$), the strain at yield ($\gamma_y$), from which you can calculate a shear stress at yield. The other important parameter is the ultimate strain, which we’ll talk about later. The idealized stress-strain curve therefore looks like this:

The ordinary approach to solving the stress distribution within a bonded joint involves finding the points along the length of the joint at which the adhesive transitions from elastic to plastic and then solving the elastic and plastic portions of the joint separately. If you try to naively solve the equations above with an elastic-perfectly plastic adhesive model, you’ll get errors since the Jacobian become singular. For the purpose of keeping this blog post simple, we’ll cheat a little bit and give the the stress-strain curve a very small slope above the yield stress. This will eliminate the numerical issues, and as as long as this slope is small enough, it won’t affect the results very much.

With this in mind, and considering that the strain could be positive or negative, we implement a function to find the stress based on the strain as follows:

def calc_tau(gamma):
  gamma_y = 0.09  # the yield strain
  G_final = 1  # a very small slope for the upper part of the curve
  sign = np.sign(gamma)
  if np.abs(gamma) <= gamma_y:
    tau_unsigned = Ga * np.abs(gamma)
  else:
    tau_unsigned = Ga * gamma_y + \
      G_final * (np.abs(gamma) - gamma_y)
  return sign * tau_unsigned

We’ll vectorize this function so that we can calculate an array of stress values based on an array of strain values:

calc_tau_vec = np.vectorize(calc_tau)

A Nonlinear Model of the Joint

Now that we have a function to describe the way in which the adhesive creates shear stress depending on its shear strain, we can implement the solution to the differential equation again. Since the boundary conditions don’t depend on the behavior of the adhesive, we can re-use the same function for calculating the residuals of the boundary condition.

def func2(x, y):
  b = (1. / (E1 * t1) + 1. / (E2 * t2)) / ta
  return np.row_stack((
    calc_tau_vec(y[1, :]) * b,
    y[0, :]
  ))


res2 = scipy.integrate.solve_bvp(
  func2,
  bc1,
  x = np.linspace(0, L, num=50),
  y = np.zeros((2, 50))
)

Here is the strain solution that we get:

plt.plot(res2.x, res2.y[1,:])
plt.title("Elastic-Plastic Adhesive")
plt.ylabel("Shear Strain, $\\gamma$")
plt.xlabel("$x$")
plt.show()

And the corresponding adhesive shear stress solution is as follows:

plt.plot(res2.x, calc_tau_vec(res2.y[1,:]))
plt.title("Elastic-Plastic Adhesive")
plt.ylabel("Shear Stress, $\\tau$")
plt.xlabel("$x$")
plt.show()

We’ll overlay the linear and the elastic-plastic models on top of each other to clearly show the differences between the two models. First, we notice that the elastic-plastic model has flat spots in the stress distribution where the adhesive has yielded. These occur near the ends of the joint. Next, we notice that the middle of the two stress distributions look similar, but shifted: for the elastic-plastic model, the stress in the “trough” is higher because the ends of this joint take a smaller proportion of the entire load.

plt.plot(res1.x, Ga * res1.y[1,:], label="Elastic Adhesive")
plt.plot(res2.x, calc_tau_vec(res2.y[1,:]), label="Elastic-Plastic Adhesive")
plt.title("Comparison of Shear Stress for Both Models")
plt.ylabel("Shear Stress, $\\tau$")
plt.xlabel("$x$")
plt.legend()
plt.show()

Exploration

We’ll create a function that takes several of the joint parameters as arguments and returns the stress and strain distributions. We’ll use this function to explore the effect of some of the joint parameters. We’re only going to implement this for the elastic-plastic model.

def model(t1, t2, ta, L, P):
  def ode(x, y):
    b = (1. / (E1 * t1) + 1. / (E2 * t2)) / ta
    return np.row_stack((
      calc_tau_vec(y[1, :]) * b,
      y[0, :]
    ))

  def bc(ya, yb):
    return np.array([
      ya[0] - 1. / ta * (-P / w / (E2 * t2)),
      yb[0] - 1. / ta * (P / w / (E1 * t1))
    ])

  res = scipy.integrate.solve_bvp(
    ode,
    bc,
    x = np.linspace(0, L, num=50),
    y = np.zeros((2, 50))
  )

  x = res.x
  gamma = res.y[1,:]
  tau = calc_tau_vec(gamma)

  return x, gamma, tau

First, we’ll keep all of the parameters constant except that we’ll vary the load. This will show us how the stress distribution changes as we increase the load. The results aren’t surprising. At low loads, the joint is fully elastic. As the load is increased, the adhesive at the ends of the overlap start to yield. As load is increased further, the yielded area grows and the “trough” gets shallower. Finally, the joint becomes fully plastic. At this point, the joint would surely fail, but since our model doesn’t check for failure, we don’t see this.

for Pi in np.linspace(1750, 2950, num=4):
  x_i, gamma_i, tau_i = model(
    t1=t1, t2=t2, ta=ta, L=L, P=Pi
  )
  plt.plot(x_i, tau_i, label=f"P={Pi}")
plt.title("Shear Stress With Various Loads")
plt.ylabel("Shear Stress, $\\tau$")
plt.xlabel("$x$")
plt.legend()
plt.show()

Next, we’ll see what happens when we change the thickness of the upper adherend. In this example, the lower adherend has a thickness of $t_2=0.063$ and we vary the thickness of the upper adherend ($t_1$) from half this thickness to four times this thickness. As we can see, this changes the length of the two plastic zones: in the extreme case of $t_1=0.250$, there is no plastic zone on the right because the adherend carrying the load at the right end of the joint is so stiff.

for t1_i in [0.032, 0.063, .125, .250]:
  x_i, gamma_i, tau_i = model(
    t1=t1_i, t2=t2, ta=ta, L=L, P=P
  )
  plt.plot(x_i, tau_i, label=f"t1={t1_i}")
plt.title("Various Upper Adherend Thicknesses")
plt.ylabel("Shear Stress, $\\tau$")
plt.xlabel("$x$")
plt.legend()
plt.show()

Finally, we’ll see the effect of changing the overlap length. This time, we’re going to vary the overlap length $L$ and keep the average shear stress ($P/A$) constant.

for Li in np.linspace(0.5, 1.5, num=3):
  x_i, gamma_i, tau_i = model(
    t1=t1, t2=t2, ta=ta, L=Li, P=5600*Li
  )
  plt.plot(x_i, tau_i, label=f"L={Li}, P={5600*Li}")
plt.title("Various Lap Lengths, Constant $P/A$")
plt.ylabel("Shear Stress, $\\tau$")
plt.xlabel("$x$")
plt.legend()
plt.show()

Here, we see that for all three overlap lengths considered, the adhesive at the ends of the lap is plastic and that there’s an elastic “trough” in the middle of each joint. At this point, we might be tempted to declare that all of the joints are able to carry at least the same average shear stress ($P/A$), but before we do so, let’s look at the shear strain in the adhesive layer for each of these cases.

for Li in np.linspace(0.5, 1.5, num=3):
  x_i, gamma_i, tau_i = model(
    t1=t1, t2=t2, ta=ta, L=Li, P=5600*Li
  )
  plt.plot(x_i, gamma_i, label=f"L={Li}, P={5600*Li}")
plt.title("Various Lap Lengths, Constant $P/A$")
plt.ylabel("Shear Strain, $\\gamma$")
plt.xlabel("$x$")
plt.legend()
plt.show()

Here we see that the shear strain in the adhesive at the ends of the longest joint is almost 0.9. Think about what that means: the “top” of the adhesive layer has moved sideways relative to the “bottom” of the layer by an amount almost equal to the thickness of the layer. In other words, that a huge amount of strain.

The ultimate shear strain is going to depend on the type of adhesive we’re using, as well as the environmental conditions (temperature, moisture content, etc.). For a lot of adhesives, the ultimate strain is going to be somewhere in the range of $0.2$ to $0.6$. So, in the three examples shown here, the first overlap length ($L=0.5$) can probably carry this value of $P/A$, the second overlap length ($L=1.0$) might be able to carry it, but the third overlap length ($L=1.5$) almost certainly will fail. This is the reason that you can’t use the average shear stress ($P/A$) to size lap joints.

If you want to play around with this model, I’ve created a widget that implements this model.

Violin Bow Stiffness

2022-06-15T00:00:00-04:00

I’ve made a few violin bows and a couple cello bows. I’m very much a novice bow maker, but I’m learning. As I’m an engineer, I’m naturally trying to apply engineering principles to bow making, which isn’t necessarily easy since violin bows are actually very complex, despite looking quite simple.

The stiffness of a bow affects what the player is able to do with it. If a bow is too stiff, it becomes nearly unplayable; if it’s too soft, they player can’t apply much force to the string before the stick bottoms out and contacts the string (normally the hair of the bow contacts the string). The stiffness affects how much camber the bow maker must add to the stick. The wrong combination of stiffness and camber can lead to a torsional-bending buckling mode, which will make the bow unplayable. The mass and mass distribution of the bow has a large effect on playability. Plus, the aesthetics of the bow are of importance. As I said, a bow is quite complex.

The “standard” wood for making violin bows has been pernambuco for the past 250 years. However, the tree that produces this wood is endangered and hence this wood is difficult to obtain. I’ve been making bows out of other types of wood — mostly ipe and snakewood. In order for a bow made from ipe to have the same stiffness as a bow made from pernambuco, the dimensions need to be altered. Hence, having a good understanding between the taper of the stick and the resulting stiffness is important.

Taper

Henry Saint-George provides a procedure for calculating the taper of a bow based on measurements of Tourte bows ( SaintGeorge (1896) ). In this procedure, the bow is divided into 12 (unequal) segments. Referring to the figure below (reproduced from Saint-George’s book), line AC is constructed perpendicular to the bow with a length of 110 mm. A second line BD is constructed perpendicular to the stick at the other end. Saint-George indicates that the line BD is 22 mm when the total length (AB) is 700 mm. A compass is used to draw the arc Ce. A line perpendicular to the stick is then constructed starting from point e and terminating at the line CD. The compass is re-set to draw the arc fg and the process is repeated. The points A, e, g, i, k, etc. are the points at which the diameter of the bow is set. At points A and e, the diameter are set equal to one another. At points y and B, they are equal to another fixed value. The diameter at the remaining points are each decremented by a fixed value. But, since those points are not uniformly spaced, the taper is not linear, but instead accelerates along the length of the stick.

This procedure seems quite complicated. However, the keen reader might recognize that the points along the stick form a geometric series. The keen reader may also recognize that the values 22 mm and 700 mm cannot both be taken as fixed: if you change the length of the bow (which affects the slope of the line CD), you also need to change the length of line BD, otherwise the procedure described above will not produce the correct overall length.

The sum of each of these segments is given by:

$$ L = \sum_{k=0}^{12} C r^k = C\left(\frac{r^{12}-1}{r-1}\right) $$

Here, the value C is selected as 110 mm and the value of $r$ needs to be found based on the value of $L$ chosen. This can be done numerically in Python. The following code does that, then computes the points and the diameters of the bow:

import scipy.optimize

length = 700.
length_constant = 110.
d_butt = 8.6
d_head = 5.6

r = scipy.optimize.root(
    lambda r: length_constant * (r**12 - 1) / (r - 1) - length,
    22
).x[0]
print(f"Found r = {r}\n")

x_points = [0.] * 13
d_points = [0.] * 13

for i in range(13):
    if i == 0:
        d_points[i] = d_butt
    else:
        x_points[i] = length_constant * (r**i - 1) / (r - 1)
        d_points[i] = d_butt + (d_head - d_butt) * (i - 1.) / 10.
    if i == 12:
        d_points[i] = d_head

    print(f"x = {x_points[i]:.1f}, d = {d_points[i]:.2f}")

Found r = 0.8741349707251251

x = 0.0, d = 8.60
x = 110.0, d = 8.60
x = 206.2, d = 8.30
x = 290.2, d = 8.00
x = 363.7, d = 7.70
x = 427.9, d = 7.40
x = 484.0, d = 7.10
x = 533.1, d = 6.80
x = 576.0, d = 6.50
x = 613.5, d = 6.20
x = 646.3, d = 5.90
x = 675.0, d = 5.60
x = 700.0, d = 5.60

We can plot the diameter of the stick:

import matplotlib.pyplot as plt

plt.plot(x_points, d_points)
plt.title("Bow Diameter")
plt.xlabel("x")
plt.show()

Stiffness

Section Properties

Bows are either (approximately) round or octagonal in cross-section. The area moment of inertia of each of these are as follows ( Oberg et al. (2000) ):

Shape	Area Moment of Inertia
Circle	$\frac{\pi d^4}{64} = 0.0490874 d^4$
Octagon	$\frac{2 d^2 \tan\frac{\pi}{8}}{12}\left[\frac{d^2 \left(1 + 2 \cos^2\frac{\pi}{8}\right)}{4\cos^2\frac{\pi}{8}}\right] = 0.0547379 d^4$

Of course, when determining the stiffness of the bow, the modulus of elasticity also needs to be known. From my research, the modulus of elasticity of pernambuco is about 30 GPa. From my measurements, the modulus of elasticity of ipe is about 20 GPa.

Finite Element Method

In order to determine the stiffness of the stick, we’ll use the finite element method with tapered beam elements. This analysis will be done in two dimensions. We’ll define a node at each of the x points found in the previous calculation of bow taper with a tapered beam element connecting adjacent nodes. The diameter (or width across flats in the case of an octagonal cross-section) is known at each of the nodes. Our model will assume that the variation in the diameter is linear between nodes.

The following derivation is based on Chapter 3 from Cook et al. (2001) , but differs since the elements are tapered beams instead of constant section beams.

Each node will have two degrees of freedom: a transverse displacement and a rotation. The degrees of freedom associated with a single element (which connects two nodes) is thus:

$$ [d] = \left[ \matrix{ \nu_1 & \theta_1 & \nu_2 & \theta_2 } \right] $$

Some of the algebra that we’ll use in the following derivation gets a bit tedious, so we’ll use the symbolic mathematics package sympy to help us:

import sympy

# Due to the way that my blogging platform works, we need to
# define a new function for printing symbolic math:
def sym_print(x):
    print('$${}$$'.format(sympy.printing.latex(x)))

The shape function for our element is a function of the element length $L$ and the position along the element $x$ and is given by:

L = sympy.var("L")
x = sympy.var("x")

B = sympy.Matrix([[
    -6 / L**2 + 12 * x / L**3,
    -4 / L + 6 * x / L**2,
    6 / L**2 - 12 * x / L**3,
    -2 / L + 6 * x / L**2
]])
sym_print(B)

$$\left[\begin{matrix}- \frac{6}{L^{2}} + \frac{12 x}{L^{3}} & - \frac{4}{L} + \frac{6 x}{L^{2}} & \frac{6}{L^{2}} - \frac{12 x}{L^{3}} & - \frac{2}{L} + \frac{6 x}{L^{2}}\end{matrix}\right]$$

For the purpose of stiffness calculations, we’re idealizing the taper of the bow so that within each element the taper is linear. This means that the diameter of the stick at the point $x$ is given by the following. Note that in this section, $x$ and $L$ refer to the distance along the length of the element dn the length of the element, respectively, rather than the dimensions of the bow.

$$ d = d_1 + \frac{x}{L}\left(d_2 - d_1\right) $$

where $d_1$ and $d_2$ are the diameters at nodes 1 and 2, respectively. So that we don’t have to carry around so many variables, we’ll define the variable $\beta$ such that:

$$ d = d_1 + \beta x $$

As we found earlier, for both circular sections and octagonal sections, the moment of inertia ($I$) is a function of $d^4$. We’ll define a new variable $\alpha$ such that:

$$ EI = \alpha d^4 $$

Combining the previous two equations and entering this into sympy, we get:

alpha = sympy.var("\\alpha")
d1 = sympy.var("d_1")
beta = sympy.var("\\beta")
EI = alpha * (d1 + beta * x)**4
sym_print(EI)

$$\alpha \left(\beta x + d_{1}\right)^{4}$$

The stiffness matrix for the element is given by:

$$ [k] = \int_0^L \left[B\right]^T EI \left[B\right] dx $$

Solving and simplifying this using sympy, we get the following. The stiffness matrix is a 4x4 matrix that is quite complex, so we’ll show one column at a time in this post:

k = sympy.simplify(
    sympy.integrate(B.T * EI * B, (x, 0, L))
)

# The first column
sym_print(k[:,0])

$$\left[\begin{matrix}\frac{12 \alpha \left(11 L^{4} \beta^{4} + 49 L^{3} \beta^{3} d_{1} + 84 L^{2} \beta^{2} d_{1}^{2} + 70 L \beta d_{1}^{3} + 35 d_{1}^{4}\right)}{35 L^{3}}\frac{2 \alpha \left(19 L^{4} \beta^{4} + 84 L^{3} \beta^{3} d_{1} + 147 L^{2} \beta^{2} d_{1}^{2} + 140 L \beta d_{1}^{3} + 105 d_{1}^{4}\right)}{35 L^{2}}\frac{12 \alpha \left(- 11 L^{4} \beta^{4} - 49 L^{3} \beta^{3} d_{1} - 84 L^{2} \beta^{2} d_{1}^{2} - 70 L \beta d_{1}^{3} - 35 d_{1}^{4}\right)}{35 L^{3}}\frac{2 \alpha \left(47 L^{4} \beta^{4} + 210 L^{3} \beta^{3} d_{1} + 357 L^{2} \beta^{2} d_{1}^{2} + 280 L \beta d_{1}^{3} + 105 d_{1}^{4}\right)}{35 L^{2}}\end{matrix}\right]$$

# The second column
sym_print(k[:,1])

$$\left[\begin{matrix}\frac{2 \alpha \left(19 L^{4} \beta^{4} + 84 L^{3} \beta^{3} d_{1} + 147 L^{2} \beta^{2} d_{1}^{2} + 140 L \beta d_{1}^{3} + 105 d_{1}^{4}\right)}{35 L^{2}}\frac{4 \alpha \left(3 L^{4} \beta^{4} + 14 L^{3} \beta^{3} d_{1} + 28 L^{2} \beta^{2} d_{1}^{2} + 35 L \beta d_{1}^{3} + 35 d_{1}^{4}\right)}{35 L}\frac{2 \alpha \left(- 19 L^{4} \beta^{4} - 84 L^{3} \beta^{3} d_{1} - 147 L^{2} \beta^{2} d_{1}^{2} - 140 L \beta d_{1}^{3} - 105 d_{1}^{4}\right)}{35 L^{2}}\frac{2 \alpha \left(13 L^{4} \beta^{4} + 56 L^{3} \beta^{3} d_{1} + 91 L^{2} \beta^{2} d_{1}^{2} + 70 L \beta d_{1}^{3} + 35 d_{1}^{4}\right)}{35 L}\end{matrix}\right]$$

# The third column
sym_print(k[:,2])

$$\left[\begin{matrix}\frac{12 \alpha \left(- 11 L^{4} \beta^{4} - 49 L^{3} \beta^{3} d_{1} - 84 L^{2} \beta^{2} d_{1}^{2} - 70 L \beta d_{1}^{3} - 35 d_{1}^{4}\right)}{35 L^{3}}\frac{2 \alpha \left(- 19 L^{4} \beta^{4} - 84 L^{3} \beta^{3} d_{1} - 147 L^{2} \beta^{2} d_{1}^{2} - 140 L \beta d_{1}^{3} - 105 d_{1}^{4}\right)}{35 L^{2}}\frac{12 \alpha \left(11 L^{4} \beta^{4} + 49 L^{3} \beta^{3} d_{1} + 84 L^{2} \beta^{2} d_{1}^{2} + 70 L \beta d_{1}^{3} + 35 d_{1}^{4}\right)}{35 L^{3}}\frac{2 \alpha \left(- 47 L^{4} \beta^{4} - 210 L^{3} \beta^{3} d_{1} - 357 L^{2} \beta^{2} d_{1}^{2} - 280 L \beta d_{1}^{3} - 105 d_{1}^{4}\right)}{35 L^{2}}\end{matrix}\right]$$

# The fourth column
sym_print(k[:,3])

$$\left[\begin{matrix}\frac{2 \alpha \left(47 L^{4} \beta^{4} + 210 L^{3} \beta^{3} d_{1} + 357 L^{2} \beta^{2} d_{1}^{2} + 280 L \beta d_{1}^{3} + 105 d_{1}^{4}\right)}{35 L^{2}}\frac{2 \alpha \left(13 L^{4} \beta^{4} + 56 L^{3} \beta^{3} d_{1} + 91 L^{2} \beta^{2} d_{1}^{2} + 70 L \beta d_{1}^{3} + 35 d_{1}^{4}\right)}{35 L}\frac{2 \alpha \left(- 47 L^{4} \beta^{4} - 210 L^{3} \beta^{3} d_{1} - 357 L^{2} \beta^{2} d_{1}^{2} - 280 L \beta d_{1}^{3} - 105 d_{1}^{4}\right)}{35 L^{2}}\frac{4 \alpha \left(17 L^{4} \beta^{4} + 77 L^{3} \beta^{3} d_{1} + 133 L^{2} \beta^{2} d_{1}^{2} + 105 L \beta d_{1}^{3} + 35 d_{1}^{4}\right)}{35 L}\end{matrix}\right]$$

We can now write a function that outputs the stiffness matrix for a tapered beam element:

import numpy as np

def elm_k(L, d1, d2, alpha):
    b = (d2 - d1) / L
    return 2 * alpha / (35*L**3) * np.array(
        [
            [
                6*(11*L**4*b**4 + 49*L**3*b**3*d1 + 84*L**2*b**2*d1**2 + 70*L*b*d1**3 + 35*d1**4),
                L*(19*L**4*b**4 + 84*L**3*b**3*d1 + 147*L**2*b**2*d1**2 + 140*L*b*d1**3 + 105*d1**4),
                6*(-11*L**4*b**4 - 49*L**3*b**3*d1 - 84*L**2*b**2*d1**2 - 70*L*b*d1**3 - 35*d1**4),
                L*(47*L**4*b**4 + 210*L**3*b**3*d1 + 357*L**2*b**2*d1**2 + 280*L*b*d1**3 + 105*d1**4)
            ],
            [
                L*(19*L**4*b**4 + 84*L**3*b**3*d1 + 147*L**2*b**2*d1**2 + 140*L*b*d1**3 + 105*d1**4),
                2*L**2*(3*L**4*b**4 + 14*L**3*b**3*d1 + 28*L**2*b**2*d1**2 + 35*L*b*d1**3 + 35*d1**4),
                L*(-19*L**4*b**4 - 84*L**3*b**3*d1 - 147*L**2*b**2*d1**2 - 140*L*b*d1**3 - 105*d1**4), 
                L**2*(13*L**4*b**4 + 56*L**3*b**3*d1 + 91*L**2*b**2*d1**2 + 70*L*b*d1**3 + 35*d1**4)
            ],
            [
                6*(-11*L**4*b**4 - 49*L**3*b**3*d1 - 84*L**2*b**2*d1**2 - 70*L*b*d1**3 - 35*d1**4),
                L*(-19*L**4*b**4 - 84*L**3*b**3*d1 - 147*L**2*b**2*d1**2 - 140*L*b*d1**3 - 105*d1**4),
                6*(11*L**4*b**4 + 49*L**3*b**3*d1 + 84*L**2*b**2*d1**2 + 70*L*b*d1**3 + 35*d1**4),
                L*(-47*L**4*b**4 - 210*L**3*b**3*d1 - 357*L**2*b**2*d1**2 - 280*L*b*d1**3 - 105*d1**4)
            ],
            [
                L*(47*L**4*b**4 + 210*L**3*b**3*d1 + 357*L**2*b**2*d1**2 + 280*L*b*d1**3 + 105*d1**4),
                L**2*(13*L**4*b**4 + 56*L**3*b**3*d1 + 91*L**2*b**2*d1**2 + 70*L*b*d1**3 + 35*d1**4),
                L*(-47*L**4*b**4 - 210*L**3*b**3*d1 - 357*L**2*b**2*d1**2 - 280*L*b*d1**3 - 105*d1**4),
                2*L**2*(17*L**4*b**4 + 77*L**3*b**3*d1 + 133*L**2*b**2*d1**2 + 105*L*b*d1**3 + 35*d1**4)
            ]
        ]
    )

Stroup Test

The Stroup Test is a way of testing the stiffness of the stick of a bow. In this test, the bow is mounted in a jig that supports the stick on two rollers that are 575 mm apart. A transverse force of 2 lb is applied mid-way between the two rollers and the deflection at the force application point is measured. From what I can tell, there were a small number of people advocating this test some time ago, but it has since become quite uncommon — most makers will assess the stiffness of a stick by feel. However, the Stroup Test can be easily implemented using the finite element method for the purpose of assessing relative stiffness of sticks made from different materials with different dimensions.

Implementing the Stroup Test

We already have a list of nodal locations. We’ll choose one of these nodes as the location of one of the supports (we’ll use the second last node for this). We’ll need to ensure that there are two other nodes for the load application point and the other support in the correct location. We’ll likely need to create these nodes and sub-divide the existing elements. We can do this in Python as follows:

x_nodes = []
d_nodes = []

x_s2 = x_points[11]
x_s1 = x_s2 - 575
x_l = x_s2 - 575 / 2

nid_s1 = -1  # storage for node ID of support #1
nid_s2 = -1  # storage for node ID of support #2
nid_l = -1   # storage for node ID of load application

tol = lambda xa, xb: abs(xa - xb) < 1e-3
inside = lambda x, xa, xb: (x - xa) * (x - xb) < 0

for x1, x2, d1, d2 in zip(
    x_points, x_points[1:], d_points, d_points[1:]):
    x_nodes.append(x1)
    d_nodes.append(d1)

    if tol(x_s1, x1):
        nid_s1 = len(x_nodes) - 1
    elif inside(x_s1, x1, x2):
        x_nodes.append(x_s1)
        d_nodes.append(d1 + (x_s1 - x1) / (x2 - x1) * (d2 - d1))
        nid_s1 = len(x_nodes) - 1

    if tol(x_s2, x1):
        nid_s2 = len(x_nodes) - 1
    elif inside(x_s2, x1, x2):
        x_nodes.append(x_s2)
        d_nodes.append(d1 + (x_s2 - x1) / (x2 - x1) * (d2 - d1))
        nid_s2 = len(x_nodes) - 1

    if tol(x_l, x1):
        nid_l = len(x_nodes) - 1
    elif inside(x_l, x1, x2):
        x_nodes.append(x_l)
        d_nodes.append(d1 + (x_l - x1) / (x2 - x1) * (d2 - d1))
        nid_l = len(x_nodes) - 1

x_nodes.append(x_points[-1])
d_nodes.append(d_points[-1])

We can now build a stiffness matrix for the model. There are now 15 nodes and each node has 2 DOF, so the matrix will be 30 x 30. We’ll use a sparse matrix. We’ll assume that all elements are round and the material has a modulus of 30 GPa.

k_model= np.zeros((2 * len(x_nodes), 2 * len(x_nodes)))

for i, (x1, x2, d1, d2) in enumerate(
    zip(x_nodes, x_nodes[1:], d_nodes, d_nodes[1:])):
    # Each element connects the two adjacent nodes
    k_elm = elm_k(
        L = x2 - x1,
        d1 = d1,
        d2 = d2,
        alpha = 0.0490874 * 30e3
    )
    for ii in range(4):
        for jj in range(4):
            k_model[i * 2 + ii, i * 2 + jj] += k_elm[ii,jj]

We can visualize the stiffness matrix. As expected, all of the elements are near the diagonal.

plt.matshow(k_model)
plt.title("Visualization of Stiffness Matrix")

Text(0.5, 1.0, 'Visualization of Stiffness Matrix')

Next, we will create the load vector. This vector will have all elements set to zero except for the entry corresponding to the first DOF of the loading node.

p_model = np.zeros(2 * len(x_nodes))
p_model[nid_l * 2] = -8.9075  # 2 lb in N

Next, we’ll take away the constrained DOFs from the stiffness matrix and the load vector. In our case, those DOFs are the transverse displacement of the constrained nodes.

mask = [i for i, _ in enumerate(p_model)
        if i != nid_s1 * 2 and i != nid_s2 * 2]

p_const = p_model[mask]
k_const = k_model[mask, :]
k_const = k_const[:, mask]

Now, we can solve for the deflections:

import scipy.linalg
d_const = scipy.linalg.solve(k_const, p_const)

Now, we can add back in the constrained DOFs into the displacement solution. These will be zero because these DOFs were constrained.

d_model = np.zeros(2 * len(x_nodes))
d_model[mask] = d_const

Now, we can plot the results:

plt.plot(x_nodes, d_model[0::2])
plt.grid()
plt.title("Deflection")
plt.xlabel("x")
plt.ylabel("Vertical Deflection")
plt.show()

Stroup values are normally given in thousandths of an inch, which we can calculate as follows:

-d_model[nid_l * 2] / 25.4 * 1000

301.20020904559294

Conclusion

This blog post describes a way of numerically finding the relationship between the stiffness of a violin bow and its taper. We used the finite element method to do so. I’m planning on developing an online calculator for performing this computation. I plan to use an early version of py-script to do so, but since I’ve never used py-script, it’s possible that it will take a while to figure it out.

Speeding up Quadrature

2021-09-18T00:00:00-04:00

Up until recently, I hadn’t really thought about the way that numerical integration was performed. Sure, I knew about some techniques like using the trapezoid rule to perform numerical integration, and without thinking about it too much, I had just assumed that the integration routines like R’s integrate function used this technique too. But, I was wrong — most libraries that implement numerical integration use adaptive quadrature.

Adaptive quadrature is actually a rather interesting technique. I won’t go into too much detail here, but the function being integrated (the integrand) is evaluated at a number of points within the integration range, and the function values are multiplied by a set of weights. In mathematical terms:

$$ \int_a^b f\left(x\right) dx \approx \sum_i^n w_i f\left(x_i\right) $$

Where the weights, $w_i$ and the evaluation points, $x_i$ are tabulated values. These values can be taken from references such as Abramowitz (1972) .

The GNU Scientific Library uses two different sets of $w_i$ and $x_i$: the first set are 15-point Kronrod weights, and the second set are 7-point Gausian weights. The estimate of the integral is computed using these two sets of weight and the absolute value of the difference between the two results is an upper bound on the error.

If the error is too great, the range is sub-divided and the integral of each sub-divided range is summed to produce the complete integral — as are the error estimates. This sub-division procedure is the “adaptive” part of adaptive quadrature.

I’ve been working on a computational problem that involves the computation of an expression of the following form:

$$ \frac{ \int_{-\infty}^\lambda g(t)A(t)dt + \int_{\lambda}^\infty h(t)A(t)dt }{ \int_{-\infty}^{\infty}A(t)dt } $$

In my particular problem, $A(t)$ is expensive to compute, while $g(t)$ and $h(t)$ are relatively computationally cheap.

In my use case, I need to compute this integral many times with slightly different $g(t)$ and $h(t)$ functions, but with the $A(t)$ function identical each time.

For now, let’s ignore the integration bounds for these four integrals. We’ll revisit the bounds shortly. The quadrature estimate of first integral (containing $g(t)$) will be:

$$ \int g(t) A(t) dt \approx \sum_i^n w_i f(x_i) = \sum_i^n w_i g(x_i) A(x_i) $$

Thus, we can pre-compute the values of $A(x_i)$ once and avoid computing them again. A similar procedure can be used for the other three integrals in the original expression.

I’ve implemented this approach of pre-computing the values of $A(x_i)$ in C++. I’ve run this several times with different repetitions and compared the speed to a “naive” approach where the complete integration is performed each time. The results are as follows:

Repetitions	Naive Approach	Pre-Computing $A(x_i)$
1	0.485 ms	1.625 ms
10	3.65 ms	1.42 ms
100	57.4 ms	1.34 ms
1000	317 ms	1.24 ms
10000	2645 ms	1.11 ms

Using the naive approach, the time scales roughly linearly with the number of repetitions, while the approach where we pre-compute the value of $A(x_i)$ is roughly constant regardless of the number of repetitions. The specific values shown here are based on a single-run of the code, so the results will be affected the whatever else my PC was doing at the time, but we can still see general trends.

Returning to the discussion of the integration bounds: first, the bounds of the two integrals in the numerator and the bounds of the integral in the denominator are all different. To account for this, we compute the integral using the widest bounds, subdivide the range as required to achieve a suitable error estimate. Then for the smaller range, we choose the subdivisions that are within the new range, adding a smaller subdivision at one end if needed.

Second, you’ll notice that some of the integration bounds are infinite. This is handled by a clever trick that I would not have though of myself — a change of variables. In my code, I’ve used a $\tan$ transformation; in the GNU Scientific Library, they use a different transform that contains a singularity. This singularity is okay if you’re not altering the integration bounds after starting the computation (which GSL does not), but can lead to trouble otherwise. After this change of variables, the integration becomes:

$$ \int_{-\infty}^\infty f(x) dx = \int_{-\pi/2}^{\pi/2} f(\tan(t)) \cos^2(t) dt $$

With this transformation, the integration bounds become finite.

Using these few tricks, the quadrature for this particular problem can be sped up significantly. These tricks won’t work for all problems, though.

Long-Running Vignettes for R Packages

2021-06-21T00:00:00-04:00

I’m going to release a new version of cmstatr soon. This new version includes, amongst other things, a new vignette. In R packages, a vignette is a type of long-form documentation. This particular vignette includes a simulation study that helps to demonstrate the validity of a particular statistical method. This simulation study takes a long time to run, though. It takes long enough that I don’t want to sit and wait for it to run every time I check that package, and I don’t want to waste resources on the CRAN servers and force their servers to re-run my vignette every time they check that package.

Jeroen Oorms wrote a blog post at rOpenSci about this topic. I decided to follow the advice in that blog post and pre-compute the new vignette on my computer, and avoid having to re-run it every time the package is checked. The blog post doesn’t include all of the necessary information for vignettes that include graphs, though. This present blog post is intended to fill in that gap.

The basic idea is that you take your long-running vignette and rename it with the extension .Rmd.orig so that R (and CRAN) doesn’t try to build it, because it doesn’t recognize it as an RMarkdown file. Then you write a script that invokes knitr to to run the executable code in the vignette and write a .Rmd file where the code is no longer executable. With this approach, when R tried to re-build the vignette, none of the code is executable, and it runs almost instantly.

In the case of the new vignette being added to cmstatr, the filename of the vignette is hk_ext.Rmd.

The first step is easy. Just rename the vignette from hk_ext.Rmd to hk_ext.Rmd.orig.

If were were to run the function knitr::knit("hk_ext.Rmd.orig", output = "hk_ext.Rmd"), it would create the .Rmd file with the executable code turned into non-executable code, and with the results of the code included. The figures would be located in folder figures/ and referenced by the resulting markdown file. However, the path to figures/ will be relative to the current working directory. This is a problem, since the current working directory will (likely) be the root directory of the package, and the vignettes are stored in the vignettes/ sub-folder.

We can fix this problem by using the following script to re-build the vignette. I’ve saved this script with the very verbose filename rebuild-long-running-vignette.R.

old_wd <- getwd()

setwd("vignettes/")
knitr::knit("hk_ext.Rmd.orig", output = "hk_ext.Rmd")
knitr::purl("hk_ext.Rmd.orig", output = "hk_ext.R")

setwd(old_wd)

This sets the working directory to the vignettes/ sub-folder, rebuilds the vignette then sets the working directory back to what it originally was.

We also need to make a change to the setup chunk of our vignette (hk_ext.Rmd.orig). This will tell knitr to put the resulting figures in the same folder as the vignette, rather than a sub-folder.

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = ""  # Added this line to the standard setup chunk
)

Now to rebuild the vignette, you just run the script rebuild-long-running-vignette.R. This script should be added to .Rbuildignore so that it doesn’t get included in the built package. Similarly, the .Rmd.orig file needs to be added to the .Rbuildignore file.

The other issue is remembering to update the vignette, now that it’s not automatic. I personally use devtools to release packages to CRAN. When you run devtools::release() it asks you a bunch of standard questions. It’s possible to add extra questions according to the documentation. So, I’ve added the following un-exported function to the package:

release_questions <- function() {
  c(
    "Did you re-build the hk_ext.Rmd using `rebuild-long-running-vignette.R`?"
  )
}

Calculating Extended Hanson—Koopmans Tolerance Limits

2021-06-12T00:00:00-04:00

Calculating tolerance limits — such as A-Basis and B-Basis is an important part of developing and certifying composite structure for aircraft. When the data doesn’t fit a convenient parametric distribution like a Normal or Weibull distribution, one often resorts to non-parametric methods. Several non-parametric methods exist for determining tolerance limits.

Vangel’s 1994 paper Vangel (1994) discusses a non-parametric method for determining tolerance limits. This article provides a brief summary of that work and discusses the implementation of that method in the R language, as well as some choices that can be made in the implementation.

This method of calculating non-parametric tolerance limits is an extension of the Hanson—Koopmans method Hanson and Koopmans (1964) . The lower tolerance can be calculated using the following formula:

$$ T_L = x_{(j)}\left[\frac{x_{(i)}}{x_{(j)}}\right]^z $$

where $x_{(i)}$ and $x_{(j)}$ indicate the $i$th and $j$th order statistic of the sample (that is, the $i$th smallest and the $j$th smallest value).

The values of $j$ and $z$ need to be determined somehow.

There is a function $H(z)$ defined as follows:

$$ H(z) = Pr\left[T(z) \ge \log(1 - \beta)\right] $$

where $\beta$ is the content of the desired tolerance limit. The details are outside the scope of this article, but we can write a function that solves the following equation for $z$.

$$ H(z) = \gamma $$

where $\gamma$ is the confidence of the desired tolerance limit.

It turns out that we obtain different values of $z$ depending on which values of $i$ and $j$ we choose.

Vangel’s approach is to set $i=1$ in all cases, then to find the value of $j$ that would produce a tolerance limit that is nearest to the population quantile assuming that the data is distributed according to a standard normal distribution.

We’ll investigate this approach through simulation. First, we’ll load a few packages.

library(tidyverse)
library(cmstatr)

Next, we’ll set the value of $i=1$ and a value of the content and confidence of our tolerance limit. We’ll choose B-Basis tolerance limits as an example.

i <- 1
p <- 0.90
conf <- 0.95

The expected value of the $i$th order statistic for a normally distributed sample can be calculated using the following function (see Harter (1961) ). We’ll need this function soon.

expected_order_statistic <- function(i, n) {
  int <- function(x) {
    x * pnorm(-x) ^ (i - 1) * pnorm(x) ^ (n - i) * dnorm(x)
  }
  integral <- integrate(int, -Inf, Inf)
  stopifnot(integral$message == "OK")
  factorial(n) / (factorial(n - i) * factorial(i - 1)) * integral$value
}

When using Vangel’s approach, we need to minimize the value of the following function.

fcn <- function(j, n) {
  e1 <- expected_order_statistic(i, n)
  e2 <- expected_order_statistic(j, n)
  z <- hk_ext_z(n, i, j, p, conf)
  abs(z * e1 + (1 - z) * e2 - qnorm(p))
}

We can plot the above function versus $j$ for the value of $n=17$:

data.frame(
  j = seq(7, 11, by = 0.1)
) %>%
  mutate(fcn = Vectorize(fcn)(j, 17)) %>%
  ggplot(aes(x = j, y = fcn)) +
  geom_line() +
  geom_point(
    data = data.frame(j = 7:11) %>%
      mutate(fcn = Vectorize(fcn)(j, 17)),
    mapping = aes(x = j, y = fcn)
  )

In this particular case, we can see that $j=9$ produces the minimum value of this function (for integer values of $j$). But this function at $j=8$ is not much worse.

Of note, there is a table of optimum values of $j$ for various values of $n$ published in CMH-17-1G [@CMH-17-1G] ¹. For most values of $n$, the optimum value from the function above matches the published value. However, for samples of size 17, 20, 23, 24 and 28, the function above disagrees with the published values by one unit. We will focus the simulation effort on samples of these sizes. For sample sizes of interest, the following values of $j$ and $z$ are published in CMH-17-1G.

published_r_n <- tribble(
  ~n, ~j_pub, ~z_pub,
  17, 8, 1.434,
  20, 10, 1.253,
  23, 11, 1.143,
  24, 11, 1.114,
  28, 12, 1.010
)

We can create an R function that returns the “optimum” value of $j$ where all the integer values of $j$ are considered, then the integer with the lowest value of that function is returned. Such an R function is as follows:

optim_j <- function(n) {
  j <- 2:n
  f <- sapply(2:n, function(j) Vectorize(fcn)(j, n))
  j[f == min(f)]
}

For values of $n$ of interest, we’ll generate a large number of samples (10,000) drawn from a normal distribution. We can calculate the true population quantile, since we know the population parameters. We can use the two variations of the nonparametric tolerance limit approach to calculate tolerance limits. The proportion of those tolerance limits that are below the population quantile should equal the selected confidence level. We’ll restrict the simulation to values of $n$ where we find different values of $j$ compared with those publised in CMH-17-1G.

mu_normal <- 100
sd_normal <- 6

set.seed(1234567)  # make this reproducible

sim_normal <- pmap_dfr(published_r_n, function(n, j_pub, z_pub) {
  j_opt <- optim_j(n)
  z_opt <- hk_ext_z(n, i, j_opt, p, conf)

  map_dfr(1:10000, function(i_sim) {
    tibble(
      n = n,
      x = list(sort(rnorm(n, mu_normal, sd_normal))),
      j_pub = j_pub,
      j_opt = j_opt,
      z_pub = z_pub,
      z_opt = z_opt,
    )
  }
  )
}) %>%
  rowwise() %>%
  mutate(
    T_pub = x[j_pub] * (x[i] / x[j_pub]) ^ z_pub,
    T_opt = x[j_opt] * (x[i] / x[j_opt]) ^ z_opt
  )
sim_normal

## # A tibble: 50,000 × 8
## # Rowwise: 
##        n x          j_pub j_opt z_pub z_opt T_pub T_opt
##    <dbl> <list>     <dbl> <int> <dbl> <dbl> <dbl> <dbl>
##  1    17 <dbl [17]>     8     9  1.43  1.40  85.4  85.4
##  2    17 <dbl [17]>     8     9  1.43  1.40  89.5  89.3
##  3    17 <dbl [17]>     8     9  1.43  1.40  83.2  83.5
##  4    17 <dbl [17]>     8     9  1.43  1.40  83.6  83.8
##  5    17 <dbl [17]>     8     9  1.43  1.40  83.4  83.8
##  6    17 <dbl [17]>     8     9  1.43  1.40  84.1  84.4
##  7    17 <dbl [17]>     8     9  1.43  1.40  82.6  82.9
##  8    17 <dbl [17]>     8     9  1.43  1.40  87.5  87.6
##  9    17 <dbl [17]>     8     9  1.43  1.40  83.9  83.8
## 10    17 <dbl [17]>     8     9  1.43  1.40  86.9  87.2
## # … with 49,990 more rows

We can plot the distribution of the tolerance limits that result from our R code and from the values of $j$ and $z$ published in CMH-17-1G. We see that the distributions are very similar.

sim_normal %>%
  pivot_longer(cols = T_pub:T_opt, names_to = "Approach") %>%
  ggplot(aes(x = value, color = Approach)) +
  geom_density() +
  facet_wrap(n ~ .) +
  ggtitle("Distribution of Tolerance Limits for Various Values of n")

In this article, we’re calculating the B-Basis (lower 90/95 tolerance limit). So, the population quantile that we’re approximating is:

x_p_normal <- qnorm(1 - p, mu_normal, sd_normal)
x_p_normal

## [1] 92.31069

We can now determine what proportion of the calculated tolerance limits were below the population quantile.

sim_normal %>%
  mutate(below_pub = T_pub < x_p_normal,
         below_opt = T_opt < x_p_normal) %>%
  group_by(n) %>%
  summarise(
    prop_below_pub = sum(below_pub) / n(),
    prop_below_opt = sum(below_opt) / n()
  )

## # A tibble: 5 × 3
##       n prop_below_pub prop_below_opt
##   <dbl>          <dbl>          <dbl>
## 1    17          0.964          0.967
## 2    20          0.967          0.965
## 3    23          0.960          0.960
## 4    24          0.959          0.957
## 5    28          0.954          0.954

In all cases, the tolerance limits are conservative when the data are normally distributed. Remember that we expect that 95% of the tolerance limits should be below the population quantile: here we see a slightly higher proportion than 95%.

We can repeat this with a distribution that is far from normal. Let’s try it withe the $\chi^2$ distribution.

df_chisq <- 6

set.seed(2345678)  # make this reproducible

sim_chisq <- pmap_dfr(published_r_n, function(n, j_pub, z_pub) {
  j_opt <- optim_j(n)
  z_opt <- hk_ext_z(n, i, j_opt, p, conf)

  map_dfr(1:10000, function(i_sim) {
    tibble(
      n = n,
      x = list(sort(rchisq(n, df_chisq))),
      j_pub = j_pub,
      j_opt = j_opt,
      z_pub = z_pub,
      z_opt = z_opt,
    )
  }
  )
}) %>%
  rowwise() %>%
  mutate(
    T_pub = x[j_pub] * (x[i] / x[j_pub]) ^ z_pub,
    T_opt = x[j_opt] * (x[i] / x[j_opt]) ^ z_opt
  )
sim_chisq

## # A tibble: 50,000 × 8
## # Rowwise: 
##        n x          j_pub j_opt z_pub z_opt T_pub T_opt
##    <dbl> <list>     <dbl> <int> <dbl> <dbl> <dbl> <dbl>
##  1    17 <dbl [17]>     8     9  1.43  1.40 1.00  1.03 
##  2    17 <dbl [17]>     8     9  1.43  1.40 1.35  1.34 
##  3    17 <dbl [17]>     8     9  1.43  1.40 1.39  1.40 
##  4    17 <dbl [17]>     8     9  1.43  1.40 1.39  1.28 
##  5    17 <dbl [17]>     8     9  1.43  1.40 1.39  1.43 
##  6    17 <dbl [17]>     8     9  1.43  1.40 0.283 0.297
##  7    17 <dbl [17]>     8     9  1.43  1.40 0.514 0.497
##  8    17 <dbl [17]>     8     9  1.43  1.40 0.264 0.268
##  9    17 <dbl [17]>     8     9  1.43  1.40 1.68  1.61 
## 10    17 <dbl [17]>     8     9  1.43  1.40 0.661 0.692
## # … with 49,990 more rows

The population quantile is:

x_p_chisq <- qchisq(1 - p, df_chisq)
x_p_chisq

## [1] 2.204131

The distribution of the tolerance limits calculated using the values of $j$ and $z$ that we calculate and those published. Again, the distributions are very similar.

sim_chisq %>%
  pivot_longer(cols = T_pub:T_opt, names_to = "Approach") %>%
  ggplot(aes(x = value, color = Approach)) +
  geom_density() +
  facet_wrap(n ~ .) +
  ggtitle("Distribution of Tolerance Limits for Various Values of n")

We can now determine what proportion of the calculated tolerance limits were below the population quantile.

sim_chisq %>%
  mutate(below_pub = T_pub < x_p_chisq,
         below_opt = T_opt < x_p_chisq) %>%
  group_by(n) %>%
  summarise(
    prop_below_pub = sum(below_pub) / n(),
    prop_below_opt = sum(below_opt) / n()
  )

## # A tibble: 5 × 3
##       n prop_below_pub prop_below_opt
##   <dbl>          <dbl>          <dbl>
## 1    17          0.963          0.965
## 2    20          0.959          0.959
## 3    23          0.959          0.958
## 4    24          0.955          0.955
## 5    28          0.953          0.953

Again with this distribution, we see that the tolerance limits are conservative.

Finally, let’s try again using a t-Distribution.

df_t <- 3
offset_t <- 150

set.seed(4567)  # make this reproducible

sim_t <- pmap_dfr(published_r_n, function(n, j_pub, z_pub) {
  j_opt <- optim_j(n)
  z_opt <- hk_ext_z(n, i, j_opt, p, conf)

  map_dfr(1:10000, function(i_sim) {
    tibble(
      n = n,
      x = list(sort(rt(n, df_t) + offset_t)),
      j_pub = j_pub,
      j_opt = j_opt,
      z_pub = z_pub,
      z_opt = z_opt,
    )
  }
  )
}) %>%
  rowwise() %>%
  mutate(
    T_pub = x[j_pub] * (x[i] / x[j_pub]) ^ z_pub,
    T_opt = x[j_opt] * (x[i] / x[j_opt]) ^ z_opt
  )
sim_t

## # A tibble: 50,000 × 8
## # Rowwise: 
##        n x          j_pub j_opt z_pub z_opt T_pub T_opt
##    <dbl> <list>     <dbl> <int> <dbl> <dbl> <dbl> <dbl>
##  1    17 <dbl [17]>     8     9  1.43  1.40  140.  140.
##  2    17 <dbl [17]>     8     9  1.43  1.40  147.  147.
##  3    17 <dbl [17]>     8     9  1.43  1.40  144.  144.
##  4    17 <dbl [17]>     8     9  1.43  1.40  147.  147.
##  5    17 <dbl [17]>     8     9  1.43  1.40  147.  147.
##  6    17 <dbl [17]>     8     9  1.43  1.40  145.  145.
##  7    17 <dbl [17]>     8     9  1.43  1.40  146.  146.
##  8    17 <dbl [17]>     8     9  1.43  1.40  147.  147.
##  9    17 <dbl [17]>     8     9  1.43  1.40  146.  146.
## 10    17 <dbl [17]>     8     9  1.43  1.40  143.  143.
## # … with 49,990 more rows

The population quantile is:

x_p_t <- qt(1 - p, df_t) + offset_t
x_p_t

## [1] 148.3623

The distribution of the tolerance limits using the two approaches are as follows. Again, the distributions are very similar.

sim_t %>%
  pivot_longer(cols = T_pub:T_opt, names_to = "Approach") %>%
  ggplot(aes(x = value, color = Approach)) +
  geom_density() +
  facet_wrap(n ~ .) +
  ggtitle("Distribution of Tolerance Limits for Various Values of n")

We can now determine what proportion of the calculated tolerance limits were below the population quantile.

sim_t %>%
  mutate(below_pub = T_pub < x_p_t,
         below_opt = T_opt < x_p_t) %>%
  group_by(n) %>%
  summarise(
    prop_below_pub = sum(below_pub) / n(),
    prop_below_opt = sum(below_opt) / n()
  )

## # A tibble: 5 × 3
##       n prop_below_pub prop_below_opt
##   <dbl>          <dbl>          <dbl>
## 1    17          0.958          0.959
## 2    20          0.953          0.952
## 3    23          0.953          0.953
## 4    24          0.954          0.954
## 5    28          0.953          0.953

For this distribution, the tolerance limits are still conservative.

From this simulation work, it appears that both approaches to selecting the value of $j$ preform equally well. The tolerance limits produced using each approach for a particular sample will be different, but both approaches seem to be equally valid.

The R package cmstatr contains the function hk_ext_z_j_opt which returns $j$ and $z$ for calculating tolerance limits with the optimization method described here (after version 0.8.0²). While the tolerance limits found for some particular samples may differ slightly from that produced by the tables published in CMH-17-1G, both results appear are equally valid.

It should be noted that CMH-17-1G uses $r$ and $k$ instead of $j$ and $z$ as used in this article and in Vangel’s paper. ↩
cmstatr version 0.8.0 and earlier used a slightly different function that was optimized. That version of the code produces slightly different values of $j$ for certain values of $n$. ↩

Basis Values From Censored Data

2021-02-09T00:00:00-05:00

Earlier, I wrote a post about using a likelihood-based approach to calculating Basis values. In that post, I hinted that likelihood-based approaches can be useful when dealing with censored data.

First of all, what does censoring mean? It means that the value reported is either artificially high or artificially low. There are a few reasons that this could happen. It happens often with lifetime data: with fatigue tests, you set a number of cycles at which the specimen “runs out” and you stop the test; with studies of mortality, some of the subjects will still be alive when you do the analysis. In these cases, the true value is greater than the observed result, but you don’t know by how much. These are examples of right-censored data.

Data can also be left-censored, meaning that the true value is less than the observed value. This can happen if some of the values are too small to be measured. Perhaps the instrument that you’re using can’t detect values below a certain amount.

There is also interval-censored data. This often occurs in survey data. For example, you might have data for individuals aged 40-44, but you don’t know where they fall within that range.

In this post, we’re going to deal with right-censored data.

At my day job, I often deal with data from testing of metallic inserts installed in honeycomb sandwich panel. These metallic inserts have a hole in their centers that will accept a screw. Their purpose is to allow a screw to be fastened to the panel, and the strength of this connections is one of the important considerations.

We determine the strength of the insert through testing. The usual test coupon that we use has two of these inserts installed, and we pull them away from each other to measure the shear strength. This is a convenient way of applying the load, but I’ve long thought that it must give low results. The loading of the coupon looks like this:

The reason that I’ve thought that this test method will give artificially low results is the fact that there are two inserts. The test ends when either one of these two insert fails: the other insert must be stronger than the one that failed first.

To illustrate this, let’s do a slightly silly thought experiment. Let’s imagine that we’re making a set of these coupons. We decide that we’re going to install one insert in each coupon first, then come back tomorrow and install the other insert in each coupon. Tomorrow comes around, and we decide to let the brand new intern install the second insert. The intern hasn’t yet been fully trained, and they accidentally install the wrong type of insert in teh second hole, but unfortunately they look identical to the correct type. The correct type of insert has a strength that is always $1000 lbf$, but we don’t know that yet. The wrong type of insert always has a strength of exactly $500 lbf$. When we do our tests, all of the coupons fail on the side that the intern installed (the wrong insert) and the strength of each coupon is $500 lbf$. We conclude that the mean strength of these inserts is $500 lbf$ with a very low variance.

But, we’d be wrong.

In this thought experiment, the actual mean strength of the inserts (considering both the correct and incorrect types of inserts) is $750 lbf$ and there’s actually a pretty high variance. We were simply unable to observe the strength of the stronger screws because of censoring.

In a more realistic case, we’re actually going to be dealing with parts that have strengths drawn from the same continuous distribution. As we move on, we’re going to assume that the strength of each individual insert is a random variable drawn from the same continuous distribution (that is, they are IID).

Let’s create some simulated data. We’ll start by loading a few R packages that we’ll need.

library(tidyverse)
library(cmstatr)
library(stats4)

Next, we’ll create a sample of $40$ simulated insert strengths. These will be drawn from a normal distribution with a mean of $1000$ and a standard deviation of $100$.

pop_mean <- 1000
pop_sd <- 100

set.seed(123)  # make this example reproducible

strength <- rnorm(40, pop_mean, pop_sd)

Now let’s calculate the mean of this sample. We expect it to be fairly close to 1000, and indeed it is.

mean(strength)

## [1] 1004.518

And we can also calculate the standard deviation:

sd(strength)

## [1] 89.77847

For the strength of most aircraft structures, we are concerned with a lower tolerance bound of the strength. For multiple load-path structure, we need to calculate the B-Basis strength, which is the $95/%$ lower confidence bound on the 10-th percentile of the strength.

Since we know the actual strength of all 40 inserts, we can calculate the B-Basis based on these actual insert strengths. Ideally, the B-Basis value that we calculate later will be close to this value.

basis_normal(x = strength)

## `outliers_within_batch` not run because parameter `batch` not specified
## `between_batch_variability` not run because parameter `batch` not specified

## 
## Call:
## basis_normal(x = strength)
## 
## Distribution:  Normal    ( n = 40 )
## B-Basis:   ( p = 0.9 , conf = 0.95 )
## 852.1482

Now, we’ll take these $40$ insert strengths and put them into $20$ coupons: each with two inserts. The observed coupon strength will be set to the lower of the two inserts installed in that coupon, because the coupon will fail as soon as either one of the installed inserts fails.

dat <- data.frame(
  ID = 1:20,
  strength1 = strength[1:20],
  strength2 = strength[21:40]
) %>%
  rowwise() %>%
  mutate(strength_observed = min(strength1, strength2)) %>%
  ungroup()
dat

## # A tibble: 20 × 4
##       ID strength1 strength2 strength_observed
##    <int>     <dbl>     <dbl>             <dbl>
##  1     1      944.      893.              893.
##  2     2      977.      978.              977.
##  3     3     1156.      897.              897.
##  4     4     1007.      927.              927.
##  5     5     1013.      937.              937.
##  6     6     1172.      831.              831.
##  7     7     1046.     1084.             1046.
##  8     8      873.     1015.              873.
##  9     9      931.      886.              886.
## 10    10      955.     1125.              955.
## 11    11     1122.     1043.             1043.
## 12    12     1036.      970.              970.
## 13    13     1040.     1090.             1040.
## 14    14     1011.     1088.             1011.
## 15    15      944.     1082.              944.
## 16    16     1179.     1069.             1069.
## 17    17     1050.     1055.             1050.
## 18    18      803.      994.              803.
## 19    19     1070.      969.              969.
## 20    20      953.      962.              953.

Let’s look at the summary statistics for this data:

dat %>%
  summarise(
    mean = mean(strength_observed),
    sd = sd(strength_observed),
    cv = cv(strength_observed)
  )

## # A tibble: 1 × 3
##    mean    sd     cv
##   <dbl> <dbl>  <dbl>
## 1  954.  75.1 0.0788

Hmmm. We see the mean is much lower than the mean of the individual insert strength. Remember that the mean insert strength was $1005$, but the mean strength of the coupons is $954$.

Next, we’ll naively calculate a B-Basis value from the measured strength. We’ll assume a normal distribution.

dat %>%
  basis_normal(strength_observed)

## 
## Call:
## basis_normal(data = ., x = strength_observed)
## 
## Distribution:  Normal    ( n = 20 )
## B-Basis:   ( p = 0.9 , conf = 0.95 )
## 809.1911

We’ll just keep this number in mind for now and we’ll move on to the idea of using a likelihood-based approach to calculate a B-Basis value, considering the fact that this data is censored.

The way that this data is censored might not be immediately obvious. But, each time we test one of these coupons, which contain two inserts, we actually get two pieces of data. We get the strength of one of the inserts. This is an exact value. But we also get a second piece of data. We know that the strength of the other insert is at least as high as the one that failed first. This is a right censored value.

In the previous post, I gave an expression for the likelihood function. However, that function only considers exact observations. The expression for the likelihood, considering censored data as follows (see Meeker et al. (2017) ).

$$ \mathcal{L}\left(\theta\right) = \prod_{i=1}^{n} \begin{cases} f\left(X_i;\,\theta\right) & \mbox{if } X_i \mbox{ is exact} \ F\left(X_i;\,\theta\right) & \mbox{if } X_i \mbox{ is left censored} \ 1 - F\left(X_i;\,\theta\right) & \mbox{if } X_i \mbox{ is right censored} \end{cases} $$

Where $f()$ is the probability density function and $F()$ is the cumulative density function.

We can implement a log-likelihood function based on this in R as follows:

log_likelihood_normal <- function(mu, sig, x, censored) {
  suppressWarnings(
    sum(map2_dbl(x, censored, function(xi, ci) {
      if (ci == "exact") {
        dnorm(xi, mean = mu, sd = sig, log = TRUE)
      } else if (ci == "left") {
        pnorm(xi, mean = mu, sd = sig, log.p = TRUE)
      } else if (ci == "right") {
        pnorm(xi, mean = mu, sd = sig, log.p = TRUE,
               lower.tail = FALSE)
      } else {
        stop("Invalid value of `censored`")
      }
    }))
  )
}

We can use this log-likelihood function to find the maximum-likelihood estimates (MLE) of the population parameters using the stats4::mle() function. First, we’ll find the MLE based only on the observed strength of each coupon, taken as a single exact value.

mle(
  function(mu, sig) {
    -log_likelihood_normal(mu, sig, dat$strength_observed, "exact")
  },
  start = c(1000, 100)
)

## 
## Call:
## mle(minuslogl = function(mu, sig) {
##     -log_likelihood_normal(mu, sig, dat$strength_observed, "exact")
## }, start = c(1000, 100))
## 
## Coefficients:
##        mu       sig 
## 953.70230  73.27344

(Note that the value of start is just a starting point for the numeric root finding.)

Here, we get the same value of mean that we previously calculated.

Now, we’ll repeat the MLE procedure, but now give it two pieces of data for each coupon: one exact value, and one right-censored value.

mle(
  function(mu, sig) {
    -log_likelihood_normal(mu,
                           sig,
                           c(dat$strength_observed, dat$strength_observed),
                           c(rep("exact", 20), rep("right", 20)))
  },
  start = c(1000, 100)
)

## 
## Call:
## mle(minuslogl = function(mu, sig) {
##     -log_likelihood_normal(mu, sig, c(dat$strength_observed, 
##         dat$strength_observed), c(rep("exact", 20), rep("right", 
##         20)))
## }, start = c(1000, 100))
## 
## Coefficients:
##         mu        sig 
## 1003.90717   88.51774

The mean estimated this way is remarkably close to the true value.

As we did in the previous blog post, we’ll next create a function that returns the profile likelihood based on a value of $t_p$ (the value that the proportion $p$ of the population is below).

profile_likelihood_normal <- function(tp, p, x, censored) {
  m <- mle(
    function(mu, sig) {
      -log_likelihood_normal(mu, sig, x, censored)
    },
    start = c(1000, 100)  # A starting guess
  )
  mu_hat <- m@coef[1]
  sig_hat <- m@coef[2]
  ll_hat <- log_likelihood_normal(mu_hat, sig_hat, x, censored)

  optimise(
    function(sig) {
      exp(
        log_likelihood_normal(
          mu = tp - sig * qnorm(p),
          sig = sig,
          x = x,
          censored = censored
        ) - ll_hat
      )
    },
    interval = c(0, sig_hat * 5),
    maximum = TRUE
  )$objective
}

The shape of this curve is as follows:

data.frame(
  tp = seq(700, 1000, length.out = 200)
) %>%
  rowwise() %>%
  mutate(R = profile_likelihood_normal(
    tp, 
    0.1,
    c(dat$strength_observed, dat$strength_observed),
    c(rep("exact", 20), rep("right", 20))
  )) %>%
  ggplot(aes(x = tp, y = R)) +
  geom_line() +
  ggtitle("Profile Likelihood for the 10th Percentile")

Next, we’ll find the value of $u$ that satisfies this equation:

$$ 0.05 = \frac{ \int_{-\infty}^{u}R(t_p) d t_p }{ \int_{-\infty}^{\infty}R(t_p) d t_p } $$

fn <- Vectorize(function(tp) {
  profile_likelihood_normal(
    tp,
    0.1,
    c(dat$strength_observed, dat$strength_observed),
    c(rep("exact", 20), rep("right", 20)))
})

denominator <- integrate(
  f = fn,
  lower = 0,
  upper = 1000
)

uniroot(
  function(upper) {
    trial_area <- integrate(
      fn,
      lower = 0,
      upper = upper
    )
    return(trial_area$value / denominator$value - 0.05)
  },
  interval = c(700, 1000)
)

## $root
## [1] 845.7739
## 
## $f.root
## [1] -1.8327e-08
## 
## $iter
## [1] 10
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05

This value of $846$ is much higher than the value of $809$ that we found earlier based on the coupon strength. But, this value of $846$ is a little lower than the B-Basis of $852$ that was based on the actual strength of all of the inserts installed.

One way to view the differences between these three numbers is as follows. The B-Basis strength is related to the 10-th percentile of the strength. But it is actually a confidence bound on the 10-th percentile. If we have only a little bit of information about the strength, there is a lot of uncertainty about the actual 10-th percentile, so the lower confidence bound is quite low. If we have a lot of information about the strength, the uncertainty is small, so the lower confidence bound is close to the actual 10-th percentile.

When we calculated a B-Basis from the observed coupon strength, we had 20 pieces of information. When we calculated a B-Basis from the actual insert strength, we had 40 pieces of information. When we calculated the B-Basis value considering the censored data, we had 40 pieces of information, but half that information wasn’t as informative as the other half: the exact values provide more information than the censored values.

Basis Values Using a Likelihood Approach

2021-02-09T00:00:00-05:00

All materials have some variability in their strength: some pieces of a given material are stronger than others. The design standards for civil aircraft mandate that one must account for this material variability. This is done by setting appropriate material allowables such that either $90\%$ or $99\%$ of the material will have a strength greater than the allowable with $95\%$ confidence. These values are referred to as B-Basis and A-Basis values, respectively. In the language of statistics, they are lower tolerance bounds on the material strength.

When you’re designing an aircraft part, one of the first steps is to determine the allowables to which you’ll compare the stress when determining the margin of safety. For many metals, A- or B-Basis values are published, and the designer will use those published values as the allowable. However, when it comes to composite materials, it is often up to the designer to determine the A- or B-Basis value themselves.

The most common way of calculating Basis values is to use the statistical methods published in Volume 1 of CMH-17 and implemented in the R package cmstatr (among other implementations). These methods area based on frequentest inference.

For example, if the data is assumed to be normally distributed, with this frequentest approach, you would calculate the B-Basis value using the non-central t-distribution (see, for example, Krishnamoorthy and Mathew (2008) ).

However, the frequentest approach is not the only way to calculate Basis values: a likelihood-based approach can be used as well. The book Statistical Intervals by Meeker et al. (2017) discusses this approach, among other topics.

The basic idea of a likelihood-based inference is that you can observe some data (by doing mechanical tests, or whatever), but you don’t yet know the population parameters, such as the mean and the variance. But, you can say that some possible values of the population parameters are more likely than others. For example, if you perform 18 tension tests of a material and the results are all around 100, the likelihood that the population mean is 100 is pretty high, but the likelihood that the population mean is 50 is really low. You can define a mathematical function to quantify this likelihood: this is called the likelihood function.

If you just need a point-estimate of the population parameters, you can find the highest value of this likelihood function: this is called the maximum likelihood estimate. If you need to find an interval or a bound (for example, the B-Basis, which is a lower tolerance bound), you can plot this likelihood function versus the population parameters and use this distribution of likelihood to determine a range of population parameters that are “sufficiently likely” to be within the interval.

The likelihood-based approach to calculating Basis values is more computationally expensive, but it allows you to deal with data that is left- or right-censored, and you can use the same computational algorithm for a wide variety of location-scale distributions. I’m planning on writing about calculating Basis values for censored data soon.

Example Data

For the purpose of this blog post, we’ll look at some data that is included in the cmstatr package. We’ll use this data to calculate a B-Basis value using the more traditional frequentest approach, then using a likelihood-based approach.

We’ll start by loading several R packages that we’ll need:

library(cmstatr)
library(tidyverse)
library(stats4)

Next, we’ll get the data that we’re going to use. We’ll use the “warp tension” data from the carbon.fabric.2 data set that comes with cmstatr. We’ll consider only the RTD environmental condition.

carbon.fabric.2 %>%
  filter(test == "WT" & condition == "RTD")

##    test condition batch panel thickness nplies strength modulus failure_mode
## 1    WT       RTD     A     1     0.113     14  129.224   8.733          LAB
## 2    WT       RTD     A     1     0.112     14  144.702   8.934      LAT,LWB
## 3    WT       RTD     A     1     0.113     14  137.194   8.896          LAB
## 4    WT       RTD     A     1     0.113     14  139.728   8.835      LAT,LWB
## 5    WT       RTD     A     2     0.113     14  127.286   9.220          LAB
## 6    WT       RTD     A     2     0.111     14  129.261   9.463          LAT
## 7    WT       RTD     A     2     0.112     14  130.031   9.348          LAB
## 8    WT       RTD     B     1     0.111     14  140.038   9.244      LAT,LGM
## 9    WT       RTD     B     1     0.111     14  132.880   9.267          LWT
## 10   WT       RTD     B     1     0.113     14  132.104   9.198          LAT
## 11   WT       RTD     B     2     0.114     14  137.618   9.179      LAT,LAB
## 12   WT       RTD     B     2     0.113     14  139.217   9.123          LAB
## 13   WT       RTD     B     2     0.113     14  134.912   9.116          LAT
## 14   WT       RTD     B     2     0.111     14  141.558   9.434    LAB / LAT
## 15   WT       RTD     C     1     0.108     14  150.242   9.451          LAB
## 16   WT       RTD     C     1     0.109     14  147.053   9.391          LGM
## 17   WT       RTD     C     1     0.111     14  145.001   9.318      LAT,LWB
## 18   WT       RTD     C     1     0.113     14  135.686   8.991    LAT / LAB
## 19   WT       RTD     C     1     0.112     14  136.075   9.221          LAB
## 20   WT       RTD     C     2     0.114     14  143.738   8.803      LAT,LGM
## 21   WT       RTD     C     2     0.113     14  143.715   8.893      LAT,LAB
## 22   WT       RTD     C     2     0.113     14  147.981   8.974      LGM,LWB
## 23   WT       RTD     C     2     0.112     14  148.418   9.118      LAT,LWB
## 24   WT       RTD     C     2     0.113     14  135.435   9.217      LAT/LAB
## 25   WT       RTD     C     2     0.113     14  146.285   8.920      LWT/LWB
## 26   WT       RTD     C     2     0.111     14  139.078   9.015          LAT
## 27   WT       RTD     C     2     0.112     14  146.825   9.036      LAT/LWT
## 28   WT       RTD     C     2     0.110     14  148.235   9.336      LWB/LAB

We really care only about the strength vector from this data, so we’ll save that vectory by itself in a variable for easy access later.

dat <- (carbon.fabric.2 %>%
          filter(test == "WT" & condition == "RTD"))[["strength"]]
dat

##  [1] 129.224 144.702 137.194 139.728 127.286 129.261 130.031 140.038 132.880
## [10] 132.104 137.618 139.217 134.912 141.558 150.242 147.053 145.001 135.686
## [19] 136.075 143.738 143.715 147.981 148.418 135.435 146.285 139.078 146.825
## [28] 148.235

Frequentest B-Basis

We can use the cmstatr package to calculate the B-Basis value from this example data. We’re going to assume that the data follows a normal distribution throughout this blog post.

basis_normal(x = dat, p = 0.9, conf = 0.95)

## 
## Call:
## basis_normal(x = dat, p = 0.9, conf = 0.95)
## 
## Distribution:  Normal    ( n = 28 )
## B-Basis:   ( p = 0.9 , conf = 0.95 )
## 127.5415

So using this approach, we get a B-Basis value of $127.54$.

Likelihood-Based B-Basis

The first step in implementing a likelihood-based approach is to define a likelihood function. This function is the product of the probability density function (PDF) at each observation ($X_i$), given a set of population parameters ($\theta$) (see Wasserman (2004) ).

$$ \mathcal{L}\left(\theta\right) = \prod_{i=1}^{n} f\left(X_i;\,\theta\right) $$

We’ll actually implement a log-likelihood function in R because taking a log-transform avoids some numerical issues. This log-likelihood function will take three arguments: the two parameters of the distribution (mu and sigma) and a vector of the data.

log_likelihood_normal <- function(mu, sig, x) {
  suppressWarnings(
    sum(
      dnorm(x, mean = mu, sd = sig, log = TRUE)
    )
  )
}

We can use this log-likelihood function to find the maximum-likelihood estimates (MLE) of the population parameters using the stats4::mle() function. This function takes the negative log-likelihood function and a starting guess for the parameters.

mle(
  function(mu, sig) {
    -log_likelihood_normal(mu, sig, dat)
  },
  start = c(130, 6.5)
)

## 
## Call:
## mle(minuslogl = function(mu, sig) {
##     -log_likelihood_normal(mu, sig, dat)
## }, start = c(130, 6.5))
## 
## Coefficients:
##         mu        sig 
## 139.626036   6.594905

We will be denoting these maximum likelihood estimates as $\hat\mu$ and $\hat\sigma$. They match the sample mean and sample standard deviation within a reasonable tolerance, but are not exactly equal.

mean(dat)

## [1] 139.6257

sd(dat)

## [1] 6.716047

The relative likelihood is the ratio between the value of the likelihood function evaluated at a given set of parameters to the value of the likelihood function evaluated at the MLE of the parameters. The relative likelihood would then be a function with two arguments: one for each of the parameters $\mu$ and $\sigma$. To reduce the number of arguments, Meeker et al. (2017) use a profile likelihood function instead. This the same as the likelihood ratio, but it is maximized with respect to $\sigma$, as defined below:

$$ R\left(\mu\right) = \max_\sigma \left[\frac{\mathcal{L}\left(\mu, \sigma\right)}{\mathcal{L}\left(\hat\mu, \hat\sigma\right)}\right] $$

When we’re trying to calculate a Basis value, we don’t really care about the mean as a population parameter. Instead, we care about a particular proportion of the population. Since a normal distribution (or other location-scale distributions) are uniquely defined by two parameters, Meeker et al. (2017) note that you can use two alternate parameters instead. In our case, we’ll keep $\sigma$ as one of the parameters, but we’ll use $t_p$ as the other instead. Here, $t_p$ is the value that the proportion $p$ of the population falls below. For example, $t_{0.1}$ would represent the 10-th percentile of the population.

We can convert between $\mu$ and $t_p$ as follows:

$$ \mu = t_p - \sigma \Phi^{-1}\left(p\right) $$

Given this re-parameterization, we can implement the profile likelihood function as follows:

profile_likelihood_normal <- function(tp, p, x) {
  m <- mle(
    function(mu, sig) {
      -log_likelihood_normal(mu, sig, x)
    },
    start = c(130, 6.5)
  )
  mu_hat <- m@coef[1]
  sig_hat <- m@coef[2]
  ll_hat <- log_likelihood_normal(mu_hat, sig_hat, x)

  optimise(
    function(sig) {
      exp(
        log_likelihood_normal(
          mu = tp - sig * qnorm(p),
          sig = sig,
          x = x
        ) - ll_hat
      )
    },
    interval = c(0, sig_hat * 5),
    maximum = TRUE
  )$objective
}

We can visualize the profile likelihood function:

data.frame(
  tp = seq(120, 140, length.out = 200)
) %>%
  rowwise() %>%
  mutate(R = profile_likelihood_normal(tp, 0.1, dat)) %>%
  ggplot(aes(x = tp, y = R)) +
  geom_line() +
  ggtitle("Profile Likelihood for the 10th Percentile")

The way to interpret this plot is that it’s quite unlikely that the true value of $t_p$ is 120, and it’s unlikely that it’s 140, but it’s pretty likely that it’s around 131.

However, when we’re calculating Basis values, we aren’t trying to find the most likely value of $t_p$: we’re trying to find a lower bound of the value of $t_p$.

The asymptotic distribution of $R$ is the $\chi^2$ distribution. If you’re working with large samples, you can use this fact to determine the lower bound of $t_p$. However, for the sample sizes that are typically used for composite material testing, the actual distribution of $R$ is far enough from a $\chi^2$ distribution, that you can’t actually do this.

Instead, we can use numerical integration to find the lower tolerance bound. We can find a value of $t_p$, which we’ll call $u$, where $0.05\%$ of the area under the $R$ curve is to its left. This will give the $95\%$ lower confidence bound on the population parameter. This can be written as follows. We’ll use numerical root finding to solve this expression for $u$.

$$ 0.05 = \frac{ \int_{-\infty}^{u}R(t_p) d t_p }{ \int_{-\infty}^{\infty}R(t_p) d t_p } $$

Since the value of $R$ vanishes as we move far from about 130, we won’t actually integrate from $-\infty$ to $\infty$, but rather integrate between two values are are relatively far from the peak of the $R$ curve.

We can implement this in the R language as follows. First, we’ll find the value of the denominator.

fn <- Vectorize(function(tp) {
  profile_likelihood_normal(tp, 0.1, dat)
})

denominator <- integrate(
  f = fn,
  lower = 100,
  upper = 150
)
denominator

## 4.339919 with absolute error < 8.9e-07

uniroot(
  function(upper) {
    trial_area <- integrate(
      fn,
      lower = 0,
      upper = upper
    )
    return(trial_area$value / denominator$value - 0.05)
  },
  interval = c(100, 150)
)

## $root
## [1] 127.4914
## 
## $f.root
## [1] -3.810654e-08
## 
## $iter
## [1] 14
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05

The B-Basis value that we get using this approach is $127.49$. This is quite close to $127.54$, which was the value that we got using the frequentest approach.

In a simple case like this data set, it wouldn’t be worth the extra effort of using a likelihood-based approach to calculating the Basis value, but we have demonstrated that this approach does work.

In a later blog post, we’ll explore a case where it is worth the extra effort. (Edit: that post is here)

cmstatr: Composite Material Data Statistics in R

2020-07-22T00:00:00-04:00

From what I’ve seen, a lot of the statistical analysis of data from composite material data is done in MS Excel. There are a number of very good tools for doing this analysis in MS Excel: ASAP, HYTEQ, STAT-17, and more recently, CMH17-STATS. I expect that the reason for the popularity of MS Excel for this application is that everyone in the industry has MS Excel installed on their computer and MS Excel is easy to use.

If you’ve read my blog before, you’ll know that I think that reproducibility is important for engineering calculations. In my view, this includes statistical analysis. If the analysis isn’t reproducible, how does a reviewer — either now or in the future — know if it’s right?

The current MS Excel tools are typically password protected so that users can’t view the macros that perform the calculations. I suspect that this was done with the best of intentions in order to prevent users from changing the code. But it also means that users can’t verify that the code is correct, or check if there are any unstated assumptions made.

To allow statistical analysis of composite material data using open-source software, I’ve written an package for the R programming language that implements the statistical methods described in CMH-17-1G. This package, cmstatr has been released on CRAN. There is also a brief discussion of this package in a paper published in the Journal of Open Source Software.

This R package allows statistical analysis to be performed using open-source tools — which can be verified by the user — and facilitates statistical analysis reports to be written at the same time that the analysis is performed by using R-Notebooks (see my earlier post).

I’ve tried to write the functions in a consistent manner so that it’s easier to learn how to use the package. I’ve also written functions to work well with the tidyverse set of packages.

There are some examples of how to use the cmstatr package in this vignette.

I hope that people find this package useful. If you use this package and find a bug, have feedback or would like a feature added, please raise an issue on GitHub.

Tracking Issues using Jupyter Notebooks

2020-04-12T00:00:00-04:00

Edit (26-May-2022): This post is largely obsolete, now that GitHub is able to render math in Markdown documents, including issues. I’m keeping this post up for historical reasons, but I’d now recommend that you now use GitHub Issues directly and include mathematica notation as needed.

I’m currently collaborating on a paper. My collaborator and I are writing the paper using LaTeX and we’re using git to track and share changes to the manuscript. We currently have a shared repository on GitHub.

GitHub has a lot of great features for collaborating on software — after all that’s why it was developed. The “Issues” features in a repository is a particularly useful feature. This allows you to discuss problems, and track the resolution of those problems. Text formatting is supported in GitHub Issues using Markdown. In many flavors of markdown, you can also embed math using LaTeX syntax. Unfortunately, GitHub flavored markdown does not support math (Edit: Note that GitHub flavored markdown now does support math). This is probably fine for the vast majority of software projects. However, it is a problem when we’re trying to discuss a mathematical model.

Several people on the internet have suggested various solutions to this shortcoming. Some have suggested using an external engine to render your math as an image, then embed that image in markdown. This works, but I think it’s cumbersome.

Several others have suggested using a Jupyter Notebook, which GitHub does actually render. I think that this is a better solution, and this is the solution that I’m planning on using with my collaborator.

Implementation Summary

In our git repository, I’m creating a folder called issues-open. Inside this folder is a set of Jupyter Notebooks, one per issue. Each collaborator can review these Notebooks, which conveniently get rendered on the GitHub web interface. When a collaborator has something to add to the issue, they can fire up their Jupyter instance and make some changes — either by adding new cells to the bottom of the notebook, or making changes to the existing text — and committing and pushing the changes. We’ve adopted the practice of starting each cell with a heading with the name of the author of that cell. This way, the Notebook looks a bit like a conversation.

Launching Jupyter Notebooks

We’re using a conda environment for Python so that we’re synced up on the versions of each package we’re using. So, the first step will be creating the conda environment from the environment YAML file. In our case, this would look like this:

conda env create -f environment.yml

This only needs to be done once on each computer. Once that’s been done, you just need to activate the environment. This is basically just telling your terminal that you want to use that version of Python. This can be accomplished like the following (obviously, replace the name of the environment with the correct name):

conda activate my-environment

Now, you can launch the Jupyter Notebook session using the following. Your web browser should pop up and allow you to create new notebooks and edit existing notebooks in the browser once you run this command.

jupyter notebook

Collaborating on Issues

The Jupyter Notebook interface is relatively straight forward and doesn’t need much discussion here. Most of the important features are available through the menus. There are keyboard shortcuts that come in handy, which can be found here.

Jupyter notebooks comprise a set of cells. The basic types of cells are markdown, code and raw. We’ll ignore raw cells here. Markdown cells contain text styled using markdown syntax. Code cells contain executable code. In our case, this will all be Python code.

If there is any code in the notebook, it’s important to realize that it runs interactively. You execute one code cell at a time. You don’t have to execute them in order either. So, if the code has side effects — like changing a global variable — the order that you run the cells in makes a difference. I think it’s good practice to restart your Python interpreter and re-run all the cells before committing a notebook in git. To do this, just click Kernel / Restart & Run All. This guarantees that the cells were run in order and have repeatable output.

The other advantage to restarting the kernel and re-running all the cells before committing is to avoid extraneous changes being tracked by git. The notebook files include a counter indicating the order in which the cells were executed. The first cell to be executed will have a counter value of 1, the second will have a value of 2, etcetera. If you execute the first five cells, then execute the first one again, it will now have a counter value of 6. If you’ve been playing around with a notebook for a while, all those counters will be incremented even higher. Even if you make no real changes to the notebook, git will register these counter changes as changes that need to be committed and tracked. You really only want the real changes to be tracked, and the easiest way to do this is to ensure that the code cells are executed in order starting from an execution count of one.

Closing an Issue

When it’s time to close an issue, whomever closes the issue simply moves the Jupyter Notebook discussing the issue to a folder called issues-closed. This should be a git-mv so that the history is maintained.

As an example, to close the issue discussed in the Notebook reorder-model-development.ipynb, the command would be:

git mv issues-open/reorder-model-development.ipynb issues-closed/

Pandoc Report Templates

2019-10-29T00:00:00-04:00

The main benefit of using Notebooks (R Notebooks or Jupyter Notebooks) is that the document is reproducible: the reader knows exactly how the results of the analysis were obtained. I wrote about the use of Notebooks in an earlier post.

Most organizations have a certain report format: a certain cover sheet layout, a certain font, a log of revisions, etcetera. For the most part, organizations have an MS Word template for this report format. If you want to use a Notebook for you analysis and to write your report, you have a few options:

You could write front matter in MS Word using your company’s report template and then attach the Notebook as an appendix.
You could also use Pandoc (more about what this is later) to convert the Notebook into a .docx file and then merge it into the report template.
You could create your own Pandoc template to convert a Notebook directly into a PDF with the correct formatting.

The first option of attaching a Notebook as an appendix to a report otherwise created in MS Word is effective but is means that you need to maintain two different files: the MS Word report and the Notebook itself. The second option of exporting the Notebook to MS Word and merging it into the template is problematic when it comes to document revisions. If the part of the analysis is revised, there is a temptation to change the affected part by either only re-exporting that section from the Notebook into docx, or worse, making the change directly in MS Word. In both cases, there is the possibility of breaking the reproducibility. For example, let’s say that in your report you define some constants at the beginning and do some math using these constants:

P = 1000
A1 = 2
A2 = 4

sigma1 = P / A1
print(sigma1)
# 500

sigma2 = P / A2
print(sigma2)
# 250

Now let’s say that you ask your new intern to revise the document so that $P = 1200$. They just edit the MS Word version of the report thinking that they will save some time. They don’t notice that $P$ is used twice in the calculation and only update the result from the first time it’s used. Now the report reads:

P = 1200
A1 = 2
A2 = 4

sigma1 = P / A1
print(sigma1)
# 600

sigma2 = P / A2
print(sigma2)
# 250

The report is now wrong. In a simple case like this, you’ll probably notice the error when you review your intern’s work, but if the math was significantly more complex, there is probably a fairly good chance that you wouldn’t pick up on the newly introduced error.

For this reason, I think that the best option is to create a Pandoc template for your company’s report template. This means that you’ll be creating a PDF directly from the Notebook. In order to revise the report, you have to re-run the Notebook — the whole Notebook.

For those unfamiliar with Pandoc, it is a program for converting between various file formats. It’s also free and open-source software. Commonly, it’s used for converting from Markdown into HTML or PDF (actually, Pandoc converts to a LaTeX format and LaTeX converts to PDF, but this happens transparently). Pandoc can also convert into MS Word (.docx) and several other formats.

When I decided to create a corporate format for use with notebooks, I looked at the types of notebooks that we use. Generally, statistics are done in an R-Notebook and other analysis is done in a Jupyter notebook. Unfortunately, R-Notebooks and Jupyter Notebooks use different templates. R-Notebooks use pandoc templates, while Jupyter uses its own template. Fortunately, there is a workaround. Jupyter is able to export to markdown, which can be read by pandoc and translated to PDF using a pandoc template. Thus, I made the decision to write a pandoc template.

When pandoc converts a markdown file to PDF, it actually uses LaTeX. The pandoc template is actually a template for converting markdown into LaTeX. Pandoc then calls pdflatex to turn this .tex file into a PDF.

When I first started figuring out how to write a template for converting markdown to PDF, I thought I was going to have to write a LaTeX class or style. I got scared. LaTeX classes are not for the faint of heart. But, I soon realized that I didn’t actually have to do that. The pandoc template that I needed to write was just a regular LaTeX document that has some parameters that pandoc can fill in. I’m not sure that I could figure out how to write a LaTeX class in a reasonable amount of time, but I sure can write a document using LaTeX. This is something that I learned to do when I wrote my undergraduate thesis, and while I don’t write LaTeX often anymore, it’s really not that hard.

A very basic LaTeX file would look something like this:

\documentclass{article}
\begin{document}

\title{My Report Title}
\author{A. Student}

\maketitle

\section{Introduction}
Some text

\end{document}

A pandoc template is just a LaTeX file, but with placeholder for the content that pandoc will insert. These placeholders are just variables surrounded with dollar signs. For example, pandoc has a variable called body. This variable will contain the body of the report. We would simply put $body$ in the part of the template where we want pandoc to insert the body of the report.

Pandoc also supports for and if statements. A common pattern is to check for the existence of a variable and use it if it does exist and use a default value if it does not. The syntax for this would look something like:

$if(myvar)$
    $myvar$
$else$
    Default text
$endif$

I’ve written the above code on multiple lines for readability, but it could be written on a single line too.

Similarly, if a variable is a list, you’d use a for statement to iterate over the list. We’ll cover this later when we talk about adding logs of revisions.

Defining New Template Variables

Pandoc defines a number of variables by default. However, you’ll likely need to define some variables of your own. First of all, you’ll likely need to define a variable for the report number and the revision.

To create the variable, it’s just a matter of defining it in the YAML header of the markdown file. Variables can either have a single value or they can be lists. Elements of a list start with dash at the beginning of the line.

Once we add the report number (which we’ll call report-no) and the revision (which we’ll call rev) to the YAML header, the YAML header will look like the following:

title: "Report Title"
author: "A. Student"
report-no: "RPT-001"
rev: B

(Bonus points if you immediately though of William Sealy Gosset when you read that).

We’ll probably want to add a log of revisions to the report. The contents of this log of revisions will have to come from somewhere, and the YAML header is the most logical place. The log of revisions will be a list with one element of the list corresponding to each revision in the log. Lists can have nested members. In our case, an entry within the log of revisions will have a revision letter, a date and a description. Including the log of revisions, the YAML header will look like this:

title: "Report Title"
author: "A. Student"
report-no: "RPT-001"
rev: B
rev-log:
-   rev: A
    date: 1-Jun-2019
    desc: Initial release
-   rev: B
    date: 18-Jun-2019
    desc: Updated loads based on fligt test data

We can now use these variables in our pandoc template. Using the variables report-no and rev are straight forward and will be just the same as using the default variables (like title and author).

Using the list variables will require the use of a for statement. In the case of a log of revisions, each revision will get a row in a LaTeX table. Using the variable rev-log, this table will look like this:

\begin{tabular}{| m{0.25in} | m{0.95in} | m{4.0in} |}
    \hline
    Rev Ltr & Date & Description \\
    $for(rev-log)$
        \hline
        $rev-log.rev$ & $rev-log.date$ & $rev-log.desc$ \\
    $endfor$
    \hline
\end{tabular}

In the above LaTeX code, everything between $for(...)$ and $endfor$ gets repeated for each item in the list rev-log. We can access the nested members using dot notation.

Using the Pandoc Template from an R-Notebook

RStudio handles a lot of the interface with pandoc. Adding the following to the YAML header of the R-Notebook should cause RStudio to use your new template when it compiles the R-Notebook to PDF. This should be all you need to do.

output:
  pdf_document:
    template: my_template_file.tex
    toc_depth: 3
    fig_caption: true
    keep_tex: false
    df_print: kable

Using the Pandoc Template from a Jupyter Notebook

Using your new pandoc template from a Jupyter Notebook is a bit more complicated because Jupyter doesn’t work directly with pandoc. First of all, we need to tell nbconvert to convert to markdown. I think that it’s best to re-run the notebook at the same time (to make sure that it is, in fact, fully reproducible. You can do this using nbconvert as follows:

jupyter nbconvert --execute --to markdown my-notebook.ipynb

But, Jupyter notebooks don’t have YAML headers like R-Notebooks do, so we need a place to put all the variables that the template needs. The easiest way to do this is to create a cell at the beginning of the notebook with the cell type set as raw, then enter the YAML header into this cell, including the starting end ending fences (---). This cell would, then, have a content similar to the following. Cells of type raw simply get copied to the output, so this becomes the YAML header in the resulting markdown file.

---
title: "Report Title"
author: "A. Student"
report-no: "RPT-001"
rev: B
rev-log:
-   rev: A
    date: 1-Jun-2019
    desc: Initial release
-   rev: B
    date: 18-Jun-2019
    desc: Updated loads based on flight test data
---

Once you’ve used nbconvert to create the markdown file, you can call pandoc. You’ll have to provide the template as a command-line argument and also specify the output filename (so that pandoc knows you want a pdf) and also give the code highlighting style. The call to pandoc will look something like this.

`pandoc` my-notebook.md -N --template=my_template_file.tex -o my-notebook.pdf --highlight-style=tango

Documentation of Your Template

A “trick” that I’ve used is to add some documentation about how to use the template inside the template itself. It’s pretty unlikely that the user will actually open up the template, but it’s relatively likely that the user will forget one of the variables that the template expects. Since pandoc allows if/else statements, I’ve added the following to my template:

$if(abstract)$
    \abstract{$abstract$}
$else$
    \abstract{
        The documentation for using the template goes here
    }
$endif$

This means that if the user forgets to define the abstract variable, the cover page of the report (where the abstract normally goes in my case) will contain the documentation for the template.

Change Bars: Future Work

One of the things that I haven’t yet figured out are change bars. In my organization, we put vertical bars in the margin of reports to indicate what part of a report has been revised. There are LaTeX packages for (manually) inserting change bars into documents. However, I haven’t yet figured out how to automatically insert these into a report generated using pandoc. I’m sure there’s a way, though.

Conclusion

I hope that this demystifies the process of writing a pandoc template to allow you to create reports directly from Jupyter Notebooks or R-Notebooks in your company’s report format.

(Edited to fix a few typos)

Package Adequacy for Engineering Calculations

2019-06-29T00:00:00-04:00

If you do engineering calculations or analysis using a language like R or Python, chances are that you’re going to use some packages. Packages are collections of code that someone else has written that you can use in your code. For example, if you need to solve a system of linear equations by inverting a matrix and you’re using Python, you might use numpy. Or if you’re using R and you need to fit a linear model to some data, you would probably use the stats package.

If you’re involved in “engineering,” you need a high level of confidence that the results that you’re getting are correct. Note that in this post — and my blog in general — that when I say “engineering,” I don’t mean software engineering: I mean design and analysis of of structures or systems that have an effect on safety. I work in civil aeronautics, mainly dealing with composites, but also dealing with metallic structure regularly. Depending on the particular type of engineering that you’re engaged in and the particular problem at hand, the consequences of getting the wrong answer could be fairly severe. You better be sure that both the interpreter and the packages are correct. Probably the best way to do this is to validate the results using another method: are there published results for a similar problem that you can use as a benchmark? Perhaps you can do some physical testing? But even if you’re doing you due diligence and validating the results somehow, you will still waste a lot of your time if there were a problem with either the interpreter or one of the packages.

Compiled languages — like C or FORTRAN —are compiled into machine code that runs directly on the processor. Interpreted languages, like Python, R or JavaScript, are not compiled into machine code, but instead an interpreter (a piece of software) reads each line of code and figures out how to run it when you run the code (not ahead of time). As far as interpreters go, if you’re using CPython (the “standard” Python interpreter) or GNU-R (the “standard” R interpreter), I think there is a rather low risk that there are any errors in the interpreter. These interpreters are written by a bunch of smart people, and both are open source, so the code that makes up the interpreters themselves are read by a much larger group of smart people. Furthermore, both interpreters are widely used and have been around for a while, so it’s very likely that significant bugs that are likely to change the result of an engineering calculation would have been found by users by now and would have been fixed.

Packages are more of a risk than interpreters are. Again, if you’re using a very widely used package that has been around for a while, like numpy (in Python) or stats (in R), there’s a pretty good chance that any bugs that would affect your calculations would have been found by now — and packages like these are maintained by groups of dedicated people.

If you’re using R, chances are that you’re getting your packages from CRAN. You should be reading the CRAN page for the package that you’re using. You can find an example of such a page here. There are a few things that you should look for to help you evaluate the reliability of the package (in addition to reference manual and any vignettes that explain how to use the package). The first is the priority of the package. Not all packages have a priority, but if the priority is “base” or “recommended,” the package is maintained by the r-core team and is almost certainly used by a lot of people. You can be fairly comfortable with these packages.

The second thing that you should look at on the CRAN page for a package is the CRAN Checks. CRAN will test all the packages every time a new version of R is released and it tests all the packages routinely to determine if a change in one package caused errors in another packages. You can see an example CRAN Check for my package rde here.

This practice is called continuous integration. It does all of these checks on several different operating systems — Windows, OSX, and several Linux distributions. If you open the CRAN Checks results for a package, you’ll see a table of all the various combinations of R version and operating system that have been tested along with the amount of time that it took to run the test and a status for each. If the Status is “OK,” then there were no errors identified. If the Status is “NOTE,” “WARNING,” or “ERROR.” There might be something wrong and it may or may not be serious. If you click on the Status link, you’ll see details and can evaluate for yourself.

I think that these CRAN checks are actually a very strong point for the R ecosystem. It ensures that package maintainers know when something outside of their package breaks their code. And, it enforces a certain level of quality: package maintainers are given a certain amount of time to fix errors, and if they don’t the package gets removed from CRAN.

The CRAN checks do a few things. First, they check that the package can, in fact, be loaded (maybe there’s an error that prevents you from using it at all). There are a few other things that it does, but the most important in terms of reliability of the package is that the CRAN checks will run any test created by the package maintainer. These tests are called unit tests. They are test that determine if the code in the package actually has the expected behavior. Package maintainers don’t have to write unit tests, but the good ones do. You can look at what tests the package maintainer has written by downloading the code of the package (you can download it from CRAN). The test are in a folder called tests. Tests basically work by providing some input to the package’s functions, and checking that the result is correct. For R packages, the testthat framework is a popular testing framework. For packages that use the testthat framework, you’ll see a number of statements that use the expect_... family of functions. Some of these tests will likely ensure that the package works at all — checking things like the return type for functions, or that a function actually does raise certain errors when invalid arguments are passed to it. Some of the tests should also ensure that the package provides correct results. When I write tests for a package, I always write both types of tests. For the tests that ensure that the results are correct, I often either check cases that have closed-form solutions, or check that the code in the package produces results that are approximately equal to example results published in articles or books. You’ll need to read through the tests to decide if they provide enough assurance that the package is correct.

If you decide that the tests for a package are not sufficient, you have three options.

You could choose not the use that package: maybe there is another that does something similar.
You can write tests yourself and contribute those tests back to the package maintainer. After all, R packages are open-source and users are encouraged to contribute back to the community. Most package maintainers would be happy to receive a patch that adds more tests: writing tests is not fun, and most people would be grateful if someone else offers to do it.
You could also manually test the package. The difficulty here is ensuring that you re-test the package every time you update the version of this package on your system.

In the python world, continuous integration isn’t as well integrated into the ecosystem. Most packages that you install probably come from PyPI. As far as I know, PyPI doesn’t do any continuous integration: it’s up to the package maintainer to run their tests regularly. Package maintainers can do one of two things: they can run the tests on their own machine before releasing a new version to PyPI, or they can use a continuous integration service like Travis-CI or CircleCI. Many of the continuous integration services provide the service for free for open source projects, so many Python packages do use a continuous integration services. Packages that use a continuous integration service normally advertise it in their README file. You’ll still need to assess whether the tests are adequate, and if the package doesn’t use continuous integration, you’ll have to either run the test yourself, or trust that the package maintainer did.

If you have already written tests for your package, setting up continuous integration using Travis-CI is quite straight forward. I haven’t personally used CircleCI, but I would imagine that it’s similarly easy to use. You can see the continuous integration results from my pcakge rde on Travis-CI here.

Whether you’re using Python or R, there are ways of ensuring that the packages you use for engineering calculations are adequate for your needs. Some people seem to be a little bit scared of open source packages and software for engineering calculations, but in a lot of ways, open source software is actually better for this since you have the ability of verifying it yourself and making a decision about whether to use it.

Automating Software Validation Reports

2019-06-20T00:00:00-04:00

I’ve been working on a Python package to analyze adhesively bonded joints recently. This package will be used to analyze adhesive joints in certain aircraft structure and will be used to substantiate the design of structural repairs, amongst other uses. Because of this, the output of this package needs to be validated against test data. This validation also needs to be documented in an engineering report.

I’ve been thinking about how to do this. On one hand, I’ve been thinking about the types of (mechanical) tests that we’ll need to run to validate the model and the various test configurations that we’ll need to include in the validation test matrix. On the other hand, I’ve also been thinking about change change management of the package and ensuring that validation report stays up to date.

I’m imagining the scenario where we run the validation testing and find that the model and the test results agree within, say, 10%. Maybe that’s good enough for the purpose (depending on the direction of the disagreement). We can then write our validation report and type out the sentence “the test data and the model were found to agree within 10%.” Then, I’m imagining that we make a refinement to the model formulation and release a new version of the package that now agrees with the test data within 5%. Now, we have a validation report for the old version of the package, but no report describing the validation of the new version. We’d need to go back through the validation report, re-run the model for all the validation cases and update the report.

When we update the validation report manually, there’s probably a pretty good chance that some aspect of the update gets missed. Maybe it’s as simple as a one of the model outputs doesn’t get updated in the revised validation report. It’s also potentially rather time consuming to update this report. It would be faster to make this validation report a Jupyter Notebook (which I’ve previously talked about). I haven’t yet written about it here, but it is possible to have a Jupyter Notebook render to a PDF using a corporate report format, so it’s even possible to make this validation report look like it should (Edit: I’ve now written about this here). We could also set up a test in the package to re-run the Jupyter Notebook, and perhaps integrate it into a continuous integration system so that the Notebook gets re-run every time a change is made to the package. This would mean that the validation report is always up to date.

When you write a Jupyter Notebook, it usually has some code that produces a result — either a numeric result, or a graph — and then you have some text that you’ve written which explains the result. The problem is that this text that you’ve written doesn’t doesn’t respond to changes in the result. Sure, there are ways of automatically updating individual numbers inside the text that you’ve written, but sometimes the way that the result of the code changes warrants a change in the sentiment of the text. Maybe the text needs to change from “the model shows poor agreement with experimental results and shouldn’t be used in this case” to “the model shows excellent agreement with experimental results and has been validated.” There’s no practical way that this type of update to the text could be automated. But if the update to the result of the code in the Notebook has been automated, there’s a good chance that the text and the results from the code will end up disagreeing — especially if the report is more than a few pages.

The Solution

So, what can be done to rectify this? We want to have the ease of having the results of the code automatically update, but we want to make sure that those results and the text of the report match. One approach to this problem — and the approach that I intend to use for the adhesive joint analysis package — is to add assert statements to the Notebook. This way, if the assertion fails, the Notebook won’t automatically rebuild and our attention will be drawn to the issue.

As an example, if the text says that the model is conservative, meaning that the strain predicted by the model is higher than the strain measured by strain gauges installed on the test articles from the validation testing, we could write the following assert statement in the Jupyter Notebook:

assert(model_strain > experimental_strain)

Now, if we later make a change to the model that causes it to under-predict strain, we’ll be alerted to this and prompted to update the validation report.

Implementing the Solution

To run a Jupyter Notebook from code (for example in a test suite), I’ve use the following code in the past. This code was based on code found on The Data Incubator Blog

def _notebook_run(self, path):
    kernel_name = "python{}".format(sys.version_info[0])
    file_dir = os.path.dirname(__file__)
    errors = []

    with open(path) as f:
        nb = nbformat.read(f, as_version=4)
        nb.metadata.get("kernelspec", {})["name"] = kernel_name
        ep = ExecutePreprocessor(kernel_name=kernel_name, timeout=3600)

        try:
            ep.preprocess(nb, {"metadata": {"path": file_dir}})

        except CellExecutionError as e:
            if "SKIP" in e.traceback:
                errors.append(str(e.traceback))
            else:
                raise e

        return nb, errors


_notebook_run("file-name-of-my-notebook.ipynb")

This code will run the Notebook file-name-of-my-notebook.ipynb and will raise an error if an error is encountered. If this is inside a unittest2 or NoseTest test suite, this will cause a test failure.

Conclusion

Validating software used in a way that affects an aircraft design is very important in ensuring the safety of that design. Keeping the validation report up to date can be tedious, but can be automated using Jupyter Notebooks. The conclusions drawn in the validation report need to match the results of the software being validated. One approach to ensuring that this is always true is to add assert statements to the Jupyter Notebook that forms the validation report.

Reproducibility of Engineering Calculations

2019-06-20T00:00:00-04:00

Reproducibility in engineering work doesn’t seem to get the attention that it deserves. I can’t count the number of times that I’ve read an old engineering report in search of a particular result, only to find that the calculation that lead to that result is only barely described, or there is just a screenshot of an Excel workbook with a few input numbers and a final result. When I find things like this, it makes me a little nervous: did the original author use the correct formula when computing this result? What assumptions did the author make and neglect to document? What approximations were made? Was the original review of the report diligent enough to check this particular result?

Let’s take a hypothetical example. For simplicity, let’s assume that we’re analyzing some sort of bracket. It’s 2 inches wide, 0.125 inches thick and 5 inches long. It’s cantilievered with a load applied 2 inches from the free edge. We care about both the deflection and the maximum stress. The formulae for deflection and stress are given by Roark¹. We’ll adapt those equations slightly:

$$ \delta_a = \frac{-P}{6 E I} (2 L^3 - 3 L^2 a + a^3) $$

$$ \sigma = \frac{M_B \frac{t}{2}}{I} = \frac{P (L - a) \frac{t}{2}}{I} $$

Given these equations and the data above, we could quite easily do the calculation in an spreadsheet program like MS-Excel. But, if we want to include our calculation in a report (most likely as a screenshot of the spreadsheet), our report will probably just look like this:

This shows the “right” answer, but if you’re reviewing the report, how do you know that the answer is right? If you’re reviewing the report before it’s released, you can probably get a copy of the Excel file and check the formulae in the cells. You’ll spend a few minutes deciphering the formula to figure out if it’s correct. But, if you’re reading the report later, especially if you’re outside the company that wrote it, good luck. You’re going to have to get out a pen, paper and your calculator to repeat the calculation and figure out if it’s right. This problem is even worse if the author of the report hard coded in a few of the input values (i.e. length, width, elastic modulus, etc.) into the formulae.

There are a few ways to address this problem of reproducibility. We’ll explore two of these ways. The first is to use software like MathCAD, or it’s free alternative SMath-Studio. Both of these products are WYSIWYG math editors that are unit aware. With either of these, your could do your calculations in the MathCAD or SMath-Studio and paste a screenshot of this into your report.

Now, the input data and the formula would be shown directly in the report. The added benefit is that, since these pieces of software are unit aware, you can’t make simple unit errors —- if you forget an exponent, the units shown in the result won’t be what you expect, so you know that you’ve made a mistake.

The other way to approach this problem is to use something called a notebook. If you’re comfortable enough to write simple code in Python, you could use a jupyter notebook. If you’re doing some data analysis or statistics, you might prefer to write some code in R (though, you could use pandas if you prefer to use Python). While you use R with jupyter notebooks (as well as several other languages), in my opionion R Studio’s R Notebooks are a little bit better to work with. If you were to do the same calculation with a notebook (in this case, we’ll use a jupyter notebook and Python), it would look like this:

There are a few advantages of using a notebook. First, you can use a programming language with a little bit more power than MathCAD or SMath-Studio — if you need to do an iterative calculation or find the root of system of non-linear equations, you can do it with a language like Python or R — and do so in a way that’s not too difficult for the reader to understand. The other advantage of using a notebook is that notebooks are intended to mix code, results and text. You could actually write your whole report using a notebook! You could explain your approach to solving the problem, include the code used to solve the problem and then show the results all in the same document. No need to copy-and-paste anything and no need to store multiple files (like a word document and a SMath-Studio file).

Text written in a notebook (either a jupyter notebook or an R Notebook) is written using using something called markdown. This is a “lightweight” way of formatting text. If you want a bullet list, you just type an asterix at the beginning of each line; if you want a heading, you start the line with a hash symbol (or two for a sub-heading). And, most importantly for engineering reports, you can include formulae using LaTeX from within markdown just by enclosing the formula with two dollar signs before and after it — no need to suffer through using the MS-Word Equation Editor.

If you need a corporate format for your report, there are ways to create PDFs from either a jupyter notebook or an R Notebook using a custom format. I plan on writting about this in a later post. Stay tuned. (Edit: I’ve written about this here)

We’ve explored a few ways of making an engineering report more reproducible. Neither of the solutions explored are idea for every scenario — some scenarios are more suited to one of the solutions or the other — but both will improve many engineering reports.

W. Young and R. Budynas, Roark’s Formulas for Stress and Strain, Seventh Edition. New York: McGraw-Hill, 2002. ↩

rde: Now on CRAN

2018-07-09T00:00:00-04:00

For the last couple of years, we’ve been using the statistical programming language R when we do statistical analysis or data visualizations at work. We typically deal with small data — most of the time, our data sets are high-tens or low-hundreds of rows of data.

A lot of the time, we create R Notebooks with our analysis and visualizations. This works well for us: the R Notebook contains the code used to do the analysis, the results of the analysis and the visualizations, all in one place. This eliminates questions like: “did you remove outliers before making the graph?” Or, “did you check that the data are distributed normally before you did that test?” A reviewer of the R Notebook can see exactly what was done.

By default, the R Notebook produces an html file that you can open in your browser. You can email this html file to a colleague, and they can see your results and graphs, as well as exactly how you obtained them. If you made a logical mistake, or an inappropriate assumption, your colleague has the opportunity to find it.

There is also a button in the html file that the R Notebook gets exported to that says “Download Rmd.” This allows your colleague to open the notebook in R Studio and run your code. If you sent your data.

The one problem with just emailing R Notebooks to a colleague is that the R Notebook does not include the data. This might be okay if the data source is a file on a network, or a database that you both have access to, but in a lot of cases — at least in my work — the data is a CSV or Excel file. Now, if I want to send an R Notebook to a colleague to review, I need to remember to send the data file along with it.

Enter rde.

I wrote the package rde (which stands for Reproducible Data Embedding) to tackle this problem. This package allows you to embed data right in your R Notebook (or any other R code). It does so by compressing the data and then base-64 encoding it into an ASCII string. This string can be pasted into the R Notebook and converted back into the original data when someone re-runs the Notebook.

I won’t go into all the details of how to use the package. If you’d like to learn more, you can read the package vignette.

This isn’t the first R Package that I’ve written, but it is the first one that I’ve submitted to CRAN. When you install an R package using install.packages(), you’re installing it from CRAN. I think that CRAN is one of the best parts of the R ecosystem since it does continuous integration for all of the packages hosted there. This helps ensure that all the packages continue to work as R is updated and as other packages are updated. I’ll likely talk about this more in a future blog post.

If you’re an R user and you think that the package rde would help you in your workflow, check it out. You can install it by typing install.packages("rde") in R. If you find a bug, please file an issue on GitHub. And, if you would like to add functionality or improve it in some way, feel free to send me a pull request.

Welcome to Kloppenborg.ca

2018-06-27T22:00:00-04:00

Welcome to kloppenborg.ca

I plan to use this website as a blog where I discuss topics related to engineering, technology and whatever else I’m thinking about at the time.

If you find any of the posts here interesting, feel free to share them. If you don’t feel free to ignore them.