Bayesian thermometer

Salekhard (66.53 °N, 66.67 °E) is located in Yamalia at the Northern end of the Ural mountain range. 95% of the daily average temperatures observed at Salekhard weather station since 1882 fall within the shaded area of the above graph.  The red line at 5°C is an indicative lower threshold for plant growth. Clearly, the growing season in Yamalia is temperature limited, variable and short.

Tree ring data have been collected from specimens of living and sub-fossil larch trees (Larix Sibirica) found in the vicinity. This is known as the Polar Ural dataset. Specialists have created growth index chronologies extending back over more than a thousand years, allowing inferences to be made about local climatic history. The data confirm that annual tree growth is correlated with “summer temperature” (usually taken as the average of June-July-August).  The Polar Ural dataset includes X-ray measurements of  annual maximum latewood density (MXD). It turns out that MXD is more strongly correlated with temperature than are noisier tree ring widths, making the Polar Ural data an especially useful tree ring thermometer.

Plants have a more complex response to temperature than does, say, a mercury thermometer. In horticulture daily growth is usually taken to be zero below a base temperature threshold T_B and linear above. However this response is unlikely to be uniform over the growing season. Indeed for most plant species, warmth is more beneficial in the early “vegetative” part of the growing season than in late summer or autumn. Putting these facts together, I assume that annual tree growth is sensitive to an effective growing season temperature S defined by

(1)   \begin{equation*}S(\{T_1,...T_{D}\}) = \frac{1}{D}\sum_{i=1}^{D} w_i \begin{cases} 0 & T_i < T_B \\ T_i-T_B & T_i > T_B \end{cases}\end{equation*}

where T_i is temperature on day i and D is an interval of days which spans the growing season (e.g. for Polar Urals I took D=200 from 11 Apr to 27 Oct).  Daily weights w_i (\sum_i w_i = 1) can be modelled flexibly using a two parameter \mathrm{Beta} function (defined on the interval [0,1]):

(2)   \begin{equation*}w_i  = \mathrm{B}(\frac{i}{D}, \alpha,\alpha + \epsilon) \end{equation*}

\alpha > 1 and \epsilon = 0 produces a centrally peaked distribution, while \epsilon controls the degree of asymmetry. Positive values of \epsilon give a higher weight to earlier parts of the growing season and visa versa.

Although wood density is expected to increase with S, it cannot exceed some physical limit, implying an additional non-linearity in the relationship between S and MXD growth index. A simple model of expected MXD is:

(3)   \begin{equation*} \overline{\mathrm{MXD}} = A(1-e^{-\lambda S})\end{equation*}

For the same reason, errors about expected MXD are skewed to the downside. I assume a skew-Normal distribution:

(4)   \begin{equation*} \mathrm{MXD} \sim \mathcal{N}_{skew} \left(\overline{\mathrm{MXD}},\omega,\mu \right)\end{equation*}

For zero skew (\mu=0) this reduces to a standard normal distribution with mean \overline{\mathrm{MXD}} and standard deviation \omega.

With suitable priors, joint posterior probability distributions for the 7 parameters of Equation 1-4 were obtained using the remarkable Markov Chain Hamiltonian Monte Carlo package stan.

Rendered by

For stan aficionados, the calculation used 4 chains each with 1000 iterations of which 500 were warmup. All \hat{\mathrm{R}} values of <1.01 suggesting good convergence. The fit to MXD data is illustrated below using mean values of the parameters (showing 5% and 50% confidence intervals of the skew Normal error distribution).

A few more things.

  1. The posterior distribution of \alpha and \epsilon can be used to estimate the point in the growing season where growth is most sensitive to temperature. This turns out to be 27 June \pm 17 days. As expected, this is in the earlier part of the growing season.
  2. Prior and posterior distributions of T_B are shown in the graph below. T_B is quite well constrained by the data at biologically reasonable values i.e. within the range expected for a cold weather plant species.
  3. The posterior “summer temperature” S determined by Bayesian analysis is similar but not identical to the standard JJA summer temperature often used by specialists to reconstruct of past climates close to the tree-line in the Northern Hemisphere. The relationship is illustrated below using mean posterior parameters values.


In a follow up post I will look at the calibration problem of reconstructing S from observed \mathrm{MXD}.

Briffa et al 2013


weather GHCN


Posted in Climate, ecosystem model, Forests, Statistics | Leave a comment

On Colour

  • Chromaticity diagrams represent all colours that can be perceived by the human eye. The Luv colour space (“CIE 1976”, illustrated above at a particular luminosity L) is believed to be approximately “perceptually uniform” which means that the distance between any two colour points is proportional to the colour difference perceived by humans.
  • The colours of the rainbow (spectrum of visible light) correspond to saturated colours on the curved boundary of the uv-diagram. These monochromatic colours (indicated above at 20nm intervals) are not distributed uniformly in uv space. In other words the colour acuity of the human eye depends on wavelength.  There is a strong peak at 484nm (cyan) and another broader peak at 592nm (orange).
  • The colours of the rainbow are shown below, firstly as a linear function of wavelength and secondly in perceptually uniform space. Notice that the cyan, blue-green and orange areas expand, while green, red and violet shrink.
  • Digital photographers are familiar with tristimulus sRGB, Adobe RGB and proPhoto colour triangles. In uv space, these cover 33%, 39% and 77% of perceivable colours respectively. Of course, colours outside the sRGB gamut cannot be displayed correctly on any sRGB device. In particular, strictly speaking none of the colours of the rainbow are displayed correctly in sRGB (or ARGB). 
  • Assuming that it is perceptually uniform, the uv-diagram can be used to compute the colour error for the visible light spectrumARGB reduces the colour error at all wavelengths relative to sRGB, because the sRGB gamut is a subset of the ARGB gamut. However there is still a significant error at the 484nm colour acuity peak. proPhoto has the nice property that the colour error is zero in the vicinity of the cyan and orange colour acuity peaks. Unfortunately no proPhoto displays exist at the present time.
  • There is an excellent discussion of the difficulties in rendering the visible light spectrum here.

R Code

Computations were done using R’s colorscience and rgeos packages.


Posted in Uncategorized | Leave a comment