Estimate parameter with result and another parameter(s)

[JulianChambrier]

Hello everybody,

I wanted to know if it was possible with Aesara to estimate a parameter (dataset or associated distribution) by knowing a parameter and the result.
I illustrate with a drawing. For example, in the diagram below, I would like to know if with the result box (in red) and the box on the bottom right (all the points of a parameter), it was possible to estimate the other parameter (box on the bottom left). The estimation of this box could possibly be points in space or a probability distribution. The idea is that the parameters (the two black boxes) give the result (the red box) when combined together

[rlouf]

Are you asking if you can do optimisation to find a point estimate, or Bayesian inference to find a distribution? (Your diagram appears truncated by the way)

It would help me to understand better if you could give a mathematical description of your problem with a toy example, or link to an example implementation with another library.

[JulianChambrier]

Hello,

Ideally, I would like to find the best distribution (with the best parameters) of a parameter.

I give you the example below to try to illustrate my problem.

I am trying to illustrate for example the cost of work on a given task, represented by C. We know a data set for C of some samples (c1, c2, ..., cn). For example, 150 000€, 12 000€, ... (without any other information).
This cost C should be formalized as the product of the daily cost of people (C_P), the number of days worked on the task (N_T) and the number of people allocated to the task (N_P). This formula depends on a business expertise. Other parameters probably influence C but for the moment we try to solve with these 3 parameters.

C_P and N_T both follow probability distributions that we have estimated with the data sets we have. Here let us suppose that C_P follows a distribution of PERT with min = 900, mode = 1500, max=2500 and N_T follows a distribution of PERT with min = 4, mode =5, max=6.

We have points for C.

For N_P, we do not know much.

The problem is therefore C = C_P * N_T * N_P.

We would try to estimate a distribution for N_P. Let be :
N_P = C / (C_P * N_T). N_P can follow any distribution (in this case we would lean towards Poisson).

Would Aesara be able to do something like that. Maybe with a Bayesian inference or something else.

Thanks for your help,

Sincerely

[rlouf]

Ok, let's assume N_P is Poisson distributed, and for now that C_P and N_T follow a normal distribution with known parameters mu_c, sigma2_c and mu_n, sigma2_n respectively (you can obviously adapt this to your problem). Then you can write your probabilistic model:

import aesara.tensor as at

srng = at.random.RandomStream(0)
n_t = snrg.normal(mu_n, sigma2_n)
c_p = srng.normal(mu_c, sigma2_c)

# We need an extra assumption about the distribution of `C`, say normal:
mu_rv = srng.normal(0, 1)  # adjust this prior to your data
sigma_rv = srng.halfnomal(1.) # same
c = srng.normal(mu_rv, sigma_rv)

# We assume that the parameter of the Poisson distribution
# follows a Gamma distribution (conjugate prior)
n_p = c / (c_p * n_p)

To perform bayesian inference you will need to compute the joint probability of this model conditioned on your data, and the test value of n_p. You can use AePPL for that:

import aeppl

# `n_t` and `c_p` are stochastic variables since you use a distribution fitted to the data, 
# so are not conditioned on.

# Let us define a variable that will hold the values taken by `c`
mu_vv = mu_rv.clone()
sigma_vv = sigma_rv.clone()

logprob = aeppl.joint_logprob({c: c_data, sigma_rv: sigma_vv, mu_rv: mu_vv})

logprob is a computation graph, that you can compile to a function that can be called:

import aesara

# We want to compute the value of `logprob` given  the values`n_p_vv`
logprob_fn = aesara.function((sigma_vv, mu_vv), logprob)
print(logprob_fn(1., 1.))
# scalar value

You can use any optimization algorithm on this function in python to find the most likely value for mu_vv and sigma_vv. Alternatively you can perform full bayesian inference using AeMCMC.

Assume you went the optimization route, so you have two values mu_val and sigma_val for mu_vv and sigma_vv. You can now compute values for n_p by compiling the original model function specifying that n_p is the output, mu_rv and sigma_rv are the inputs (you will provide values for them):

sample_n_p_fn = aesara.function((mu_rv, sigma_rv), n_p)
sample_n_p_fn(mu_val, sigma_val)

[rlouf]

Sorry I misunderstood the problem, I corrected my response.

[JulianChambrier]

Hello, sorry for my late reply. Some parts of the code do not work : - AttributeError: module 'aesara.tensor' has no attribute 'div'.

Moreover, what type is waiting for c_data ?

I try to understand how to solve all this problem, but I didn't understand everything in this approach

[rlouf]

Hello, sorry for my late reply. Some parts of the code do not work : - AttributeError: module 'aesara.tensor' has no attribute 'div'.

Yeah, I always make the mistake because of #1212. Corrected.

Moreover, what type is waiting for c_data ?

That's homework :)

I try to understand how to solve all this problem, but I didn't understand everything in this approach

Is this a question about Aesara or about solving your problem in general?

[JulianChambrier]

I think I'm not quite comfortable with the Bayesian approach and how it can help us with my problem. Could you explain it to me?

[rlouf]

It's hard to summarize in a github discussion, but the gist is that you would get a distribution for the priors of $c$ and then would get a distribution for $n_p$ by sampling forward. If you're interested McElreath's Statistical Rethinking is a good introduction.

Anyway, you can also use non-bayesian methods here: you can run logprob_fn through any optimizer implemented in Scipy to get point estimates for mu and sigma and get a point estimate for n_p by passing these values to sample_n_fn.