Add inversion exampes to repo

[ashiklom]

Example script:

library(rrtm)
library(BayesianTools)
​
data("testspec", package = "PEcAnRTM")
​
obs <- testspec_ACRU[,1]
​
likelihood <- function(params) {
  mod <- prospect4(params[1], params[2], params[3], params[4])
  sum(dnorm(obs, mod[["reflectance"]], params[5], log = TRUE))
}
​
prior <- createPrior(
  density = function(x) {
    # Note: Truncated at zero
    N <- dnorm(x[1], 1, 5, log = TRUE)
    Cab <- dnorm(x[2], 40, 15, log = TRUE)
    Cw <- dlnorm(x[3], -4.456, 1.216, log = TRUE)
    Cm <- dlnorm(x[4], -5.15, 1.328, log = TRUE)
    rsd <- dexp(x[5], 10, log = TRUE)
    N + Cab + Cw + Cm + rsd
  },
  sampler = function(n = 1) {
    N <- rnorm(n, 1, 5)
    Cab <- rnorm(n, 40, 15)
    Cw <- rlnorm(n, -4.456, 1.216)
    Cm <- rlnorm(n, -5.15, 1.328)
    rsd <- rexp(n, 10)
    cbind(N, Cab, Cw, Cm, rsd)
  },
  lower = c(1, 0, 0, 0, 0)
)
setup <- createBayesianSetup(
  likelihood = likelihood,
  prior = prior
)
​
samps <- runMCMC(setup)

This should exist as a script / vignette in the repo somewhere.

[serbinsh]

In this example, what are these priors based on again? Are these general? Specific to temperate broadleaf? Could I use these for Arctic or should I consider broadening?

Just some general questions I have - I will explore this on my data and share

[ashiklom]

I think they're designed around my general guesses about the typical ranges of all these parameters. You can view the exact densities with code like:

# Cab density, from 0 to 120
curve(dnorm(x, 40, 15), 0, 120)
# Cw density from 0 to 0.2
curve(dlnorm(x, -4.456, 1.216), 0, 0.2)

You can use those kinds of plots as a basis to determine whether they are reasonable for your data or not. These priors should be OK as a first pass, but worth thinking carefully about what ranges you expect for your data.

[ashiklom]

Found a bug in my example — the residual standard deviation wasn’t actually being passed to the likelihood, so it was just sampling from the prior. It was:

sum(dnorm(obs, mod[["reflectance"]], log = TRUE))

Should be:

sum(dnorm(obs, mod[["reflectance"]], params[5], log=TRUE))

Fixed in the edited version above.

[serbinsh]

Thanks! adding the fix and retrying the inversion.

[ashiklom]

Two quick additions here:

(1) The distributions3 package is super useful for elegantly creating and sampling from the same distributions. So you might consider something like this for your priors (this should make it a bit easier to test different prior configurations):

library(distributions3)

# Define prior distribution objects
Nd <- Normal(1, 5)
Cabd <- Normal(40, 15)
Cwd <- LogNormal(-4.456, 1.216)
Cmd <- LogNormal(-5.15, 1.328)
rsdd <- Exponential(10)

prior <- createPrior(
  density = function(x) {
    # Note: Truncated at zero
    log_pdf(Nd, x[1]) +
      log_pdf(Cabd, x[2]) +
      log_pdf(Cwd, x[3]) +
      log_pdf(Cmd, x[4]) +
      log_pdf(rsdd, x[5])
  },
  sampler = function(n = 1) {
    N <- random(Nd, n)
    Cab <- random(Cabd, n)
    Cw <- random(Cwd, n)
    Cm <- random(Cmd, n)
    rsd <- random(rsdd, n)
    cbind(N, Cab, Cw, Cm, rsd)
  },
  lower = c(1, 0, 0, 0, 0)
)

(2) If you want to try a simple heteroskedastic variance model with no autocorrelation, you can try something like this:

# Priors
rsd_sloped <- Exponential(10)  # param 5
rsd_intd <- Exponential(10)  # param 6

# Likelihood
refl <- mod[["reflectance"]]
rsd <- refl * params[5] + params[6]
sum(dnorm(obs, refl, rsd, log = TRUE))

[serbinsh]

I will give these additions a try. Interestingly, thus far it seems I am getting fairly different pigment retrievals between PROSPECT-D and 5b. not sure if this is real yet as this also corresponds to moving from PEcAn RTM to rrtm etc but, for example, Cab is cut in half with 5b rrtm.

[serbinsh]

@ashiklom FYI - here is one of my current working examples, but I will clean this up into an example that uses a spec dataset saved withing rrtm.

https://github.com/TESTgroup-BNL/canopyRT/blob/main/Rscripts/PROSPECT/Invert_leaf_refl_spectra_rrtm_PROSPECT-5b.R

One thing I have done to remove depends on PEcAN is to copy the autoburnin code into the repo this script lives in. Should we just go ahead and steal that code for rrtm since it works the same for these inversions?

[ashiklom]

Sure, you're welcome to steal any of my code you'd like, especially in the name of reducing dependencies!

It may be worth taking a closer look at that autoburnin code at some point -- I think it might have a bug that might make it overly conservative (so we could reach actual convergence faster).