Engineer needs help modeling uncertainty of physical system measurement ( TypeError: Variables that depend on other nodes cannot be used for observed data. )

1

Apologies for being vague, and for my total lack of stats knowledge. :slight_smile: I will try to avoid accidental misuse of stats terms by avoiding them as much as possible. :slight_smile: And if anything is unclear I will gladly clarify as much as I am able.

I have a somewhat complicated physical system I need to know the internal state of. I can not make direct measurements of this state, but I can make indirect, approximate measurements of some function of that state. And of course the initial calibration of my system is approximate.

So far I have used non-probabilistic methods (Levenberg Marquardt) to fit a model of the system to my measurements. This works well, but I need more than a best fit. I would like a distribution of my belief of what the internal state is, given my measurements.

(My actual problem is more complicated than what I’m describing here, but I think if I can solve this I will be able to figure the rest out on my own.)

I model my system as a nonlinear deterministic function of both its internal state and my initial calibration. For now I can assume my model is correct and exactly captures the physical system. The internal state has about 10 dimensions, and the initial calibration has about 20 dimensions. I can provide reasonable estimates of my calibration error, and I know a feasible region for the internal state. I think I should be able to get by with Gaussians or other well behaved distributions.

I can take upwards of a few hundred different measurements, but there is a catch. Each measurement is of a slightly different function of my model/system. I can not repeat a measurement (except in simulation). All measurements are corrupted by noise, and I can provide a prior of this noise.

Essentially, in pseudocode I have:

# internal_state is a 10 dimensional random variable
# calibration is a 20 dimensional random variable
# functions is a list of known nonlinear deterministic functions
# sigma is some known constant

measurements = []
for i in range(len(functions)):
	func = functions[i]
	measurements.append( func(internal_state, calibration) + N(0, sigma) )

# given our measurements, what is the distribution of our belief of what internal_state is?

I have tried to look at an extremely simplified toy variant of my problem (and I’m probably going about it wrong):

mu = [2., 3.]
with model:
    internal_state = pm.Normal("internal_state", mu=mu[0])
    calibration = pm.Normal("calibration", mu=mu[1])

    def toy_function(a, b):
    	# In reality I would have many different but similar functions
	    return (1.1 + a) * (0.9 + b) + a

    measurement_error = pm.Normal("measurement_error", sigma=0.01)
    measurement = pm.Normal("measurement", 
                            mu = toy_system(mu[0], mu[1]), 
                            observed = toy_system(internal_state, calibration) + measurement_error)
    idata = pm.sample()
    # Given measurement, what do we know about internal_state ?

However this gives me an error (PyMC v5.21.1):

> TypeError: Variables that depend on other nodes cannot be used for observed data.

If I take this error message at face value it feels like a show stopper, but the restriction itself seems odd to me. I looked online for help with this error, but the discussions I found were outside my depth.

I feel like this should be a solvable problem, as I have a reasonable guess of what all the sources of error are, but I do not know how to proceed. The resources online solve a good many problems that sure seem much more complicated than this toy task, but they don’t seem engineering/physics focused, except for some regression examples that did not leave me with a good idea of how to tackle this.

I assume I am simply struggling with a knowledge gap, but if there is another tool or approach that is better suited for this please let me know.

Thank you!

2

Your thinking about the problem seems to be somewhat backwards. You never observe random variables, you observe draws from the distributions those RVs represent, then try to learn something about the latent random variables you posit in your model based on that. It fundamentally doesn’t make sense to write something like observed = f(RV), because it reverses what is “data” and what is “model”.

A PyMC model describes a generative process. What you want to write down is the way in which your measurements come about. The measurements, on the other hand, are totally fixed, and everything else in the model is going to bend to better explain how they could have arisen. This is a useful organizing principle for me when tackling a new problem, because I sometimes have to rearrange how I think about a problem.

In your case, I would start by trying to recast the whole problem in terms like measurement = f(???) + error. If error ~ N(0, sigma), then this is equivalent to measurement ~ N(f(???), sigma) due to how normal distributions behave under affine transformations. Now you have observation distribution for your model, measurement_hat = pm.Normal('measurement_hat', mu=mu, sigma=0.01, observed=measurements), you just have to figure out what goes in the mean function. I guess this is toy_function, but I don’t really understand how it is being used.

I also don’t really understand your setup with mu. Is that what is actually observed? Are internal_state and calibration unknown, or just themselves noisy? Where is the actual measured quantity in the model?