The PEcAn RTM package
Introduction
The PEcAnRTM package provides tools for analyses involving common radiative transfer models. The highlights of this package are its ability to efficiently run a suite of related leaf and canopy radiative transfer models, as well as to perform maximum likelihood and, particularly, Bayesian inversions of these models.
Installation
The easiest way to install this package is via
install_github from the devtools package.
install.packages("devtools")
library(devtools)
install_github("ashiklom/pecan", subdir = "modules/rtm")
# Defaults to branch 'master'.
# For custom branches, add `ref = "branchname"`This package relies on a modern (>= 2003) Fortran compiler, as
determined by your R installation. On most Unix systems, this is
standard, but R may specify a particular compiler version that you don’t
have, resulting in an installation error. To fix this, simply add the
following to your system ~/R/Makevars file.
FC = gfortranOverview of features
Simulating reflectance
Available radiative transfer models are called by passing a vector of their parameters. Similar models with different versions (e.g. PROSPECT 4, 5, 5B) also have a “version” argument. These models return a matrix of reflectance, transmittance, and/or absorption spectra (depending on the model) at 1 nm resolution over the wavelength range 400 to 2500 nm.
The PROSPECT family of models returns the reflectance (column 1) and transmittance (column 2) for an individual leaf as a function of 4 to 6 parameters (depending on the version).
##
## Attaching package: 'PEcAnRTM'## The following object is masked from 'package:graphics':
##
## matplotwl <- 400:2500
params <- c(
"N" = 1.4, "Cab" = 40, "Car" = 15,
"Cbrown" = 0.5, "Cw" = 0.002, "Cm" = 0.004
)
p4 <- prospect(params[c(-3, -4)], version = 4)
p4## reflectance transmittance
## 400 0.04099898 0.0001814716
## 401 0.04103723 0.0001699630
## 402 0.04106248 0.0001575602
## 403 0.04111370 0.0001452140
## 404 0.04120401 0.0001348522
## ...
## reflectance transmittance
## 2496 0.1422667 0.2801837
## 2497 0.1417965 0.2792691
## 2498 0.1420623 0.2795135
## 2499 0.1413868 0.2782765
## 2500 0.1418849 0.2788870Notice only the first and last few lines of the spectra are printed,
and the wavelength is annotated on the side. This is is because all RTM
simulations in this package are of a special matrix class
spectra.
## [1] "spectra" "matrix"This class allows traditional subsetting via index with single
brackets ([]), as well as subsetting by wavelength with
double brackets.
## reflectance
## 400 0.04099898
## 401 0.04103723
## 402 0.04106248
## 403 0.04111370
## 404 0.04120401
## ...
## reflectance
## 445 0.04545556
## 446 0.04542921
## 447 0.04540287
## 448 0.04537665
## 449 0.04533726## transmittance
## 500 0.008937186
## 501 0.010279518
## 502 0.011772740
## 503 0.013438540
## 504 0.015276985
## ...
## transmittance
## 516 0.05336301
## 517 0.05776873
## 518 0.06224375
## 519 0.06679217
## 520 0.07136411The package also provides special plotting (plot,
matplot) and combining (cbind) methods. Note
that the cbind method automatically matches wavelengths,
which facilitates working with spectra with different wavelengths.
p5 <- prospect(params[-4], version = 5)
p5b <- prospect(params, version = "5B")
p_multi <- cbind(p4, p5, p5b)
matplot(p_multi[, c(1, 3, 5)], lty = 1:3, col = 2, ylim = c(0, 1))
matplot(1 - p_multi[, c(2, 4, 6)], lty = 1:3, col = 3, add = TRUE)
legend("topright", c("Reflectance", "Transmittance"), col=c(2, 3), lty = 1)
legend("top", c("4", "5", "5B"), lty = 1:3)
The SAIL family of models returns the bidirectional (1),
hemispherical directional (2), directional hemispherical (3), and
bidirectional hemispherical (4) reflectance factors for a canopy with a
given set of approximately 20 parameters. It is often coupled to the
PROSPECT model as PRO4SAIL. (Again, note that the return type is a
spectra, which leads matplot to automatically
put wavelengths on the x axis.)
## N Cab Car Cbrown Cw Cm LIDFa LIDFb
## 1.500 40.000 8.000 0.000 0.010 0.009 -0.350 -0.150
## TypeLIDF LAI q tts tto psi psoil
## 1.000 3.000 0.010 30.000 10.000 0.000 1.000p4s <- pro4sail(sail.params)
matplot(p4s, xlab="Wavelength (nm)", ylab="Reflectance")
legend("topright", as.character(1:4), col=1:4, lty=1:4)
The above example illustrates the use of defparam to get
the default parameters for a particular model. Similarly,
model.list is a data.frame containing all
currently available models.
## modname modcode fullname
## 1 prospect_d 1153 PROSPECT D
## 2 prospect_5b 1152 PROSPECT 5B
## 3 prospect_5 1151 PROSPECT 5
## 4 prospect_4 1141 PROSPECT 4
## 5 pro4sail 2111 PRO4SAIL
## 6 pro4saild 2113 PRO4SAILD
## par.names
## 1 N Cab Car Canth Cbrown Cw Cm
## 2 N Cab Car Cbrown Cw Cm
## 3 N Cab Car Cw Cm
## 4 N Cab Cw Cm
## 5 N Cab Car Cbrown Cw Cm LIDFa LIDFb TypeLIDF_r LAI q tts tto psi psoil
## 6 N Cab Car Canth Cbrown Cw Cm LIDFa LIDFb TypeLIDF_r LAI q tts tto psi psoil
## par.default
## 1 N=1.5 Cab=40 Car=8 Canth=8 Cbrown=0 Cw=0.01 Cm=0.009
## 2 N=1.5 Cab=40 Car=8 Cbrown=0 Cw=0.01 Cm=0.009
## 3 N=1.5 Cab=40 Car=8 Cw=0.01 Cm=0.009
## 4 N=1.5 Cab=40 Cw=0.01 Cm=0.009
## 5 N=1.5 Cab=40 Car=8 Cbrown=0 Cw=0.01 Cm=0.009 LIDFa=-0.35 LIDFb=-0.15 TypeLIDF=1 LAI=3 q=0.01 tts=30 tto=10 psi=0 psoil=1
## 6 N=1.5 Cab=40 Car=8 Canth=1 Cbrown=0 Cw=0.01 Cm=0.009 LIDFa=-0.35 LIDFb=-0.15 TypeLIDF=1 LAI=3 q=0.01 tts=30 tto=10 psi=0 psoil=1Inversion
A novel feature of this package is the ability to perform a Bayesian inversion of a Radiative Transfer Model. Here are several advantages of the Bayesian approach:
- Parameter uncertainty: The output of a Bayesian analysis is a full joint probability distribution of the model parameters, which includes a robust estimate of their uncertainy and covariance between parameters.
- Prior knowledge: If previous, independent estimates of parameters are available, these can be used to inform the model.
- Partitioning variability: Random effects models provide a powerful framework for understanding the sources of variability and uncertainty in a data set.
An inversion can be performed using the invert.auto
function, which uses the Metropolis Hastings MCMC algorithm to invert an
arbitrary model. There are a lot of configuration options to
invert.auto, so the recommended way to perform an inversion
is to start with a default settings list, provided by the package
itself. In the following sample, we demonstrate the default inversion of
the PROSPECT model.
The model to be inverted is, in this case, just a one-line call to the PROSPECT 5 model with the params vector. It returns a vector of reflectance values. The requirement is that the “model” output be as long as the length (or number of rows) of the observation vector (or matrix). In this case, the PROSPECT model returns 2101 reflectance values, so our observation also has to have this many points.
## function (params, seed = NULL)
## prospect(params, 5)[, 1]
## <bytecode: 0x559981c9e8d0>
## <environment: namespace:PEcAnRTM>The recommended way to set settings is to start with a default object
and modify it. For example, in the following block, we reduce the length
of the run to make it go a little faster. threshold = ...
sets the maximum value of the multivariate Gelman Diagnostic used to
assess convergence – it defaults to 1.1, but we set it higher here for
demonstrative purposes).
invert.options$n.tries <- 1 # Number of attempts
invert.options$nchains <- 2 # Number of MCMC chains
invert.options$ngibbs <- 5000 # Number of iterations per chain
invert.options$burnin <- 1000 # Length of burnin period
invert.options$do.lsq.first <- TRUE # Initialize with results from a fast least-squares optimization algorithm
invert.options$threshold <- 1.3 # Maximum value for Gelman diagnosticThe full list of inversion options is below. See the documentation
for invert.auto for a full description of each of
these.
## [1] "model" "inits.function" "prior.function" "param.mins"
## [5] "param.maxs" "ngibbs" "nchains" "burnin"
## [9] "ngibbs.max" "ngibbs.min" "ngibbs.step" "return.samples"
## [13] "target" "do.lsq" "save.samples" "quiet"
## [17] "adapt" "adj_min" "n.tries" "do.lsq.first"
## [21] "threshold"Now, we load some test data (for Acer rubrum leaves
(testspec_ACRU, accessible through
data(testspec)).
data(testspec)
observed <- testspec_ACRU[,1]
plot(wl, observed, xlab="Wavelength", ylab="Reflectance", type='l')
To perform an inversion, just pass the observation matrix and the
inversion settings into invert.auto.
quiet = TRUE suppresses the progress bar, which is done
here to clean up the knitted document. We also set
parallel = FALSE here because this vignette is rebuilt on
shared hardware where only one processor might be available; for runs on
your own machine, you will probably want to set it TRUE to use all the
cores on your machine.
if(file.exists("inversion.output.rds")){
inversion.output <- readRDS("inversion.output.rds")
} else {
inversion.output <- invert.auto(observed = observed,
invert.options = invert.options,
quiet = TRUE,
parallel = FALSE)
saveRDS(inversion.output, "inversion.output.rds")
}## [1] "Running chain 1 of 2"## Loading required namespace: MASS## [1] "Running chain 2 of 2"
## [1] "The following parameters did not converge: N (2.319), Cab (16.325), Cw (2.014), residual (11.297), mpsrf (5.597)"
## [1] "The following parameters did not converge: N (2.496), Cab (13.627), Car (2.526), Cw (1.972), Cm (1.305), residual (11.558), mpsrf (5.561)"
## [1] "The following parameters did not converge: N (2.354), Cab (11.718), Car (4.226), Cw (2.194), residual (12.562), mpsrf (5.786)"
## [1] "The following parameters did not converge: N (2.087), Cab (11.004), Car (5.419), Cw (1.909), residual (11.331), mpsrf (5.054)"
## [1] "The following parameters did not converge: N (1.536), Cab (6.358), Car (3.944), Cw (1.902), residual (9.023), mpsrf (2.221)"
## [1] "The following parameters did not converge: Cab (5.468), Car (3.827), Cw (1.675), residual (7.280), mpsrf (1.861)"
## [1] "The following parameters did not converge: Cab (4.201), Car (3.679), Cw (1.435), residual (5.595), mpsrf (1.643)"
## [1] "The following parameters did not converge: Cab (2.984), Car (3.316), residual (4.211), mpsrf (1.449)"
## [1] "The following parameters did not converge: Cab (2.009), Car (2.882), residual (3.249), mpsrf (1.338)"
## [1] "The following parameters did not converge: Cab (1.632), Car (1.971), residual (1.924)"
## [1] "The following parameters did not converge: Cab (1.356), Car (1.443)"
## [1] "Converged with all Gelman diag <= 1.040"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.010275 1.500855 1.742733 1.078381 1.091234 1.511129 FALSE
## 47 1.014927 1.398001 1.568968 1.060891 1.088559 1.277160 FALSE
## 48 1.018952 1.348195 1.436491 1.070806 1.079634 1.105966 FALSE
## 49 1.025754 1.265303 1.290486 1.069895 1.076361 1.035175 FALSE
## 50 1.033380 1.177513 1.193167 1.060902 1.069248 1.050862 FALSE
## 51 1.039638 1.113622 1.156769 1.055004 1.067374 1.112852 FALSE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 TRUE TRUE FALSE FALSE
## 47 TRUE TRUE FALSE FALSE
## 48 TRUE TRUE FALSE FALSE
## 49 TRUE TRUE FALSE FALSE
## 50 TRUE TRUE FALSE FALSE
## 51 TRUE TRUE FALSE FALSE
## PSRF residual > 1.10
## 46 TRUE
## 47 TRUE
## 48 TRUE
## 49 FALSE
## 50 FALSE
## 51 TRUE
## [1] "Converged with all Gelman diag <= 1.024"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.025759 1.277181 1.301350 1.067686 1.076087 1.040817 FALSE
## 47 1.034443 1.188738 1.200573 1.061510 1.075279 1.043894 FALSE
## 48 1.040124 1.118901 1.155472 1.055312 1.068866 1.114430 FALSE
## 49 1.046232 1.097783 1.103697 1.033141 1.059952 1.138026 FALSE
## 50 1.055114 1.106982 1.087086 1.028798 1.069712 1.214408 FALSE
## 51 1.061628 1.154562 1.074914 1.024133 1.076420 1.247946 FALSE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 TRUE TRUE FALSE FALSE
## 47 TRUE TRUE FALSE FALSE
## 48 TRUE TRUE FALSE FALSE
## 49 FALSE TRUE FALSE FALSE
## 50 TRUE FALSE FALSE FALSE
## 51 TRUE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 FALSE
## 47 FALSE
## 48 TRUE
## 49 TRUE
## 50 TRUE
## 51 TRUE
## [1] "Converged with all Gelman diag <= 1.012"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.045105 1.112414 1.118471 1.043627 1.068482 1.130412 FALSE
## 47 1.049161 1.103049 1.090549 1.030378 1.064622 1.163718 FALSE
## 48 1.057752 1.150982 1.073077 1.026262 1.073024 1.238870 FALSE
## 49 1.062234 1.164246 1.070555 1.017689 1.063556 1.274032 FALSE
## 50 1.058353 1.159564 1.065277 1.016529 1.056889 1.300055 FALSE
## 51 1.060801 1.152422 1.065388 1.012291 1.036022 1.257024 FALSE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 TRUE TRUE FALSE FALSE
## 47 TRUE FALSE FALSE FALSE
## 48 TRUE FALSE FALSE FALSE
## 49 TRUE FALSE FALSE FALSE
## 50 TRUE FALSE FALSE FALSE
## 51 TRUE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 TRUE
## 47 TRUE
## 48 TRUE
## 49 TRUE
## 50 TRUE
## 51 TRUE
## [1] "Converged with all Gelman diag <= 1.004"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.062764 1.154530 1.065666 1.019989 1.070744 1.264882 FALSE
## 47 1.061610 1.160923 1.070236 1.017157 1.063496 1.308672 FALSE
## 48 1.056078 1.186242 1.054588 1.012277 1.040144 1.260518 FALSE
## 49 1.075979 1.129199 1.042767 1.011445 1.032864 1.249020 FALSE
## 50 1.088591 1.104867 1.016547 1.009186 1.032687 1.287773 FALSE
## 51 1.115382 1.136678 1.004468 1.009068 1.034571 1.294312 TRUE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 TRUE FALSE FALSE FALSE
## 47 TRUE FALSE FALSE FALSE
## 48 TRUE FALSE FALSE FALSE
## 49 TRUE FALSE FALSE FALSE
## 50 TRUE FALSE FALSE FALSE
## 51 TRUE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 TRUE
## 47 TRUE
## 48 TRUE
## 49 TRUE
## 50 TRUE
## 51 TRUE
## [1] "Converged with all Gelman diag <= 1.003"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.061225 1.150667 1.066094 1.012302 1.035586 1.256250 FALSE
## 47 1.083239 1.106268 1.027819 1.011887 1.032831 1.261330 FALSE
## 48 1.103394 1.128585 1.009354 1.008734 1.036283 1.304683 TRUE
## 49 1.112459 1.142599 1.001704 1.007133 1.034044 1.289415 TRUE
## 50 1.117454 1.141364 1.001764 1.006769 1.034963 1.281200 TRUE
## 51 1.130763 1.113914 1.003272 1.011068 1.035217 1.288482 TRUE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 TRUE FALSE FALSE FALSE
## 47 TRUE FALSE FALSE FALSE
## 48 TRUE FALSE FALSE FALSE
## 49 TRUE FALSE FALSE FALSE
## 50 TRUE FALSE FALSE FALSE
## 51 TRUE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 TRUE
## 47 TRUE
## 48 TRUE
## 49 TRUE
## 50 TRUE
## 51 TRUE
## [1] "Converged with all Gelman diag <= 1.006"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.115914 1.129813 1.008111 1.008225 1.035894 1.297018 TRUE
## 47 1.114047 1.143266 1.001352 1.005927 1.036907 1.283718 TRUE
## 48 1.121717 1.138197 1.004315 1.007543 1.036380 1.294847 TRUE
## 49 1.127311 1.108302 1.002617 1.011455 1.034630 1.285551 TRUE
## 50 1.137443 1.078515 1.005791 1.013889 1.037377 1.295195 TRUE
## 51 1.133841 1.062286 1.005633 1.015951 1.042321 1.298115 TRUE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 TRUE FALSE FALSE FALSE
## 47 TRUE FALSE FALSE FALSE
## 48 TRUE FALSE FALSE FALSE
## 49 TRUE FALSE FALSE FALSE
## 50 FALSE FALSE FALSE FALSE
## 51 FALSE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 TRUE
## 47 TRUE
## 48 TRUE
## 49 TRUE
## 50 TRUE
## 51 TRUE
## [1] "Converged with all Gelman diag <= 1.007"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.125280 1.135346 1.004693 1.009110 1.035286 1.296000 TRUE
## 47 1.129095 1.102366 1.003009 1.012653 1.035130 1.307844 TRUE
## 48 1.132833 1.077568 1.005720 1.014701 1.036227 1.302234 TRUE
## 49 1.129021 1.059673 1.005579 1.016999 1.041444 1.292592 TRUE
## 50 1.114802 1.045912 1.005202 1.018359 1.040384 1.277404 TRUE
## 51 1.113703 1.059761 1.006738 1.015060 1.035772 1.253700 TRUE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 TRUE FALSE FALSE FALSE
## 47 TRUE FALSE FALSE FALSE
## 48 FALSE FALSE FALSE FALSE
## 49 FALSE FALSE FALSE FALSE
## 50 FALSE FALSE FALSE FALSE
## 51 FALSE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 TRUE
## 47 TRUE
## 48 TRUE
## 49 TRUE
## 50 TRUE
## 51 TRUE
## [1] "Converged with all Gelman diag <= 1.008"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.132712 1.076216 1.006025 1.014751 1.036528 1.303374 TRUE
## 47 1.127668 1.058215 1.005808 1.017233 1.040866 1.287174 TRUE
## 48 1.116565 1.044158 1.006215 1.018161 1.040278 1.279028 TRUE
## 49 1.112307 1.055867 1.006918 1.014579 1.034386 1.244630 TRUE
## 50 1.110387 1.048077 1.008861 1.012838 1.029689 1.216196 TRUE
## 51 1.113096 1.053154 1.008077 1.011560 1.037453 1.242241 TRUE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 FALSE FALSE FALSE FALSE
## 47 FALSE FALSE FALSE FALSE
## 48 FALSE FALSE FALSE FALSE
## 49 FALSE FALSE FALSE FALSE
## 50 FALSE FALSE FALSE FALSE
## 51 FALSE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 TRUE
## 47 TRUE
## 48 TRUE
## 49 TRUE
## 50 TRUE
## 51 TRUE
## [1] "Converged with all Gelman diag <= 1.008"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.115707 1.043455 1.005330 1.017971 1.040397 1.279511 TRUE
## 47 1.112092 1.055958 1.006945 1.014511 1.034118 1.244235 TRUE
## 48 1.108442 1.049768 1.008611 1.012551 1.029771 1.212259 TRUE
## 49 1.110573 1.053321 1.008097 1.011974 1.037151 1.241525 TRUE
## 50 1.122022 1.051827 1.007975 1.013487 1.054528 1.229665 TRUE
## 51 1.109865 1.038031 1.008239 1.014229 1.054555 1.214701 TRUE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 FALSE FALSE FALSE FALSE
## 47 FALSE FALSE FALSE FALSE
## 48 FALSE FALSE FALSE FALSE
## 49 FALSE FALSE FALSE FALSE
## 50 FALSE FALSE FALSE FALSE
## 51 FALSE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 TRUE
## 47 TRUE
## 48 TRUE
## 49 TRUE
## 50 TRUE
## 51 TRUE
## [1] "Converged with all Gelman diag <= 1.010"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.111901 1.049105 1.008733 1.012413 1.029573 1.217240 TRUE
## 47 1.112767 1.053246 1.008035 1.011535 1.037285 1.242570 TRUE
## 48 1.121485 1.051919 1.007903 1.013474 1.054459 1.229249 TRUE
## 49 1.109840 1.038093 1.008251 1.014190 1.054640 1.214147 TRUE
## 50 1.113090 1.040201 1.010475 1.015610 1.054696 1.220498 TRUE
## 51 1.100429 1.057444 1.010390 1.015603 1.046901 1.162533 TRUE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 FALSE FALSE FALSE FALSE
## 47 FALSE FALSE FALSE FALSE
## 48 FALSE FALSE FALSE FALSE
## 49 FALSE FALSE FALSE FALSE
## 50 FALSE FALSE FALSE FALSE
## 51 FALSE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 TRUE
## 47 TRUE
## 48 TRUE
## 49 TRUE
## 50 TRUE
## 51 TRUE
## [1] "Converged with all Gelman diag <= 1.003"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.119547 1.058820 1.007808 1.011928 1.049633 1.238282 TRUE
## 47 1.110856 1.038769 1.008329 1.015064 1.054775 1.223306 TRUE
## 48 1.116532 1.037261 1.010106 1.015966 1.057340 1.229082 TRUE
## 49 1.100955 1.059301 1.010798 1.014756 1.045597 1.167768 TRUE
## 50 1.090820 1.059423 1.006381 1.020940 1.048724 1.150023 FALSE
## 51 1.079316 1.064668 1.003312 1.025311 1.043861 1.148721 FALSE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 FALSE FALSE FALSE FALSE
## 47 FALSE FALSE FALSE FALSE
## 48 FALSE FALSE FALSE FALSE
## 49 FALSE FALSE FALSE FALSE
## 50 FALSE FALSE FALSE FALSE
## 51 FALSE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 TRUE
## 47 TRUE
## 48 TRUE
## 49 TRUE
## 50 TRUE
## 51 TRUE
## [1] "Converged with all Gelman diag <= 1.006"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.120242 1.041585 1.010321 1.014147 1.058171 1.236081 TRUE
## 47 1.097318 1.052809 1.011760 1.014469 1.042272 1.174968 FALSE
## 48 1.097997 1.058188 1.006206 1.020250 1.051193 1.155236 FALSE
## 49 1.082401 1.065512 1.002939 1.023541 1.045353 1.137504 FALSE
## 50 1.074574 1.047845 1.002009 1.017255 1.030985 1.142216 FALSE
## 51 1.065480 1.047268 1.006134 1.009992 1.021344 1.101988 FALSE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 FALSE FALSE FALSE FALSE
## 47 FALSE FALSE FALSE FALSE
## 48 FALSE FALSE FALSE FALSE
## 49 FALSE FALSE FALSE FALSE
## 50 FALSE FALSE FALSE FALSE
## 51 FALSE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 TRUE
## 47 TRUE
## 48 TRUE
## 49 TRUE
## 50 TRUE
## 51 TRUE
## [1] "Converged with all Gelman diag <= 1.004"The output of invert.auto is a list of two objects: A
list of summary statistics for each parameter and an
mcmc.list object of the samples for diagnostic purposes or
calculation of other summary statistics.
par(mfrow=c(1,1))
samples.mat <- as.matrix(inversion.output$samples)[-(2000:0),1:5]
colnames(samples.mat) <- params.prospect5
pairs(samples.mat, pch=".")
means <- unlist(inversion.output$results[grep("mu", names(inversion.output$results))])[1:5]
prospect.sim <- prospect(means, 5)[,1] # reflectance
plot(wl, observed, type='l', col=1, xlab="wavelength (nm)", ylab="reflectance")
lines(wl, prospect.sim, type='l', col=2)
legend("topright", c("observed", "predicted"), lty=1, col=1:2)