Bayesian inversion to produce estimates using a state-space time-series model

Questions v5
1

Hi,

I’m having difficulty getting pymc3 to work when inverting the outputs of an RTM. Please could someone offer me some advice?

Ideally I’d like to take a time series of reflectance values and get a time series of associated crop traits - by using the prosail model in combination with a state-space time series model and bayesian model.

import numpy as np
import pandas as pd
import pymc3 as pm
import theano.tensor as tt
import arviz as az
import prosail
from theano.compile.ops import as_op

# Preprocessing: normalize the reflectance data reflectance has shape (,12)
reflectance_data_normalized = (reflectance - reflectance.mean()) / reflectance.std()

@as_op(itypes=[tt.dscalar, tt.dscalar, tt.dscalar, tt.dscalar, tt.dscalar, tt.dscalar], otypes=[tt.dvector])
def prosail_model(n, cab, car, cbrown, cw, cm):
    # ProSail model with parameters: n, cab, car, cbrown, cw, cm
    # Returns simulated reflectances
    reflectances = []
    reflectance = prosail.run_prosail(n, cab, car, cbrown, cw, cm, w,
                                          lidfa=-1.0, hspot=0.01, tts=30., tto=10., psi=0., factor='SDR', rsoil=0.5, psoil=0.5, prospect_version='D')
    return reflectances

# Model block
with pm.Model() as model:
    # Priors for ProSail parameters
    n_prior = pm.Uniform("n", lower=0.0, upper=1.0)
    cab_prior = pm.Uniform("cab", lower=0.0, upper=100.0)
    car_prior = pm.Uniform("car", lower=0.0, upper=100.0)
    cbrown_prior = pm.Uniform("cbrown", lower=0.0, upper=100.0)
    cw_prior = pm.Uniform("cw", lower=0.0, upper=0.5)
    cm_prior = pm.Uniform("cm", lower=0.0, upper=0.5)

    # Combine ProSail parameters
    prosail_prior = tt.stack([n_prior, cab_prior, car_prior, cbrown_prior, cw_prior, cm_prior], axis=-1)

    # Calculate ProSail model output for each time step
    prosail_output = pm.Deterministic("prosail_output", calculate_prosail(prosail_prior))

    # Likelihood - assuming normally distributed noise
    sigma_obs = pm.HalfNormal("sigma_obs", sigma=0.1)
    obs = pm.Normal("obs", mu=prosail_output, sigma=sigma_obs, observed=reflectance_data_normalized)

    # Sample
    trace = pm.sample(2000, tune=1000, target_accept=0.9)

# Plot traces
az.plot_trace(trace)

My error is “TypeError: Expected an array-like object, but found a Variable: maybe you are trying to call a function on a (possibly shared) variable instead of a numeric array?” because I am struggling to get the prosail scipy functions to cooperate with pymc3.

2

You can’t just use generic numpy/scipy code with PyMC. As you discovered, these functions expect numeric values to compute on, but a PyMC graph is comprised of symbolic pytensor Variables. Thus the error. In some cases it works, but it’s a happy accident. The computational backend for PyMC is pytensor – see here for some discussion on how it all works.

You were on the right track when using as_op – you will need to wrap your prosail model in a Pytensor Op. Well first you need to upgrade your setup if possible – pymc3 and theano are several years out of date at this point. After that, you can check out this tutorial about wrapping arbitrary numpy code in a custom Ops for guidance on how to proceed.