So the unknown parameters A_1, A_2, A_3 are sampled from a prior and then the equation

1=F\left(t,A_{1},A_{2},A_{3}\right)=A_{3}t\mathrm{plogis}\left(A_{2}t\mathrm{plogis}\left(A_{1}t^{2}\right)\right)

is solved for t and the solution is observed with negligible error.

That’s correct in theory but super narrow distribution are likely to sample very slow.

For a very narrow distribution you can approximate integrating out one parameter. Let’s say, A_3 because it’s easy to solve A_3 from the equation given t, A_1 and A_2. The observational model is

t^{\star}\sim N\left(t,\sigma\right)

and in the limit of \sigma \to 0 this looks like

A_{3}^{\star}\sim N\left(A_{3},\left|\frac{\partial A_{3}}{\partial t}\right|\sigma\right)

where A_{3}^{\star} is the value corresponding to the observed root (keeping A_1, A_2 constant).

Now integrate and the probability is

\int\mathbb{P}\left(A_{1},A_{2},A_{3}\right)\mathrm{d}A_{3}\propto\mathbb{P}\left(A_{1},A_{2},A_{3}^{\star}\right)\times\left|\frac{\partial A_{3}}{\partial t}\right|\times\sigma

The distribution is modified by a multiplier that looks like a Jacobian adjustment.

Implicit function differentiation says that

\frac{\partial A_{3}}{\partial t}=-\left(\frac{\partial F}{\partial A_{3}}\right)^{-1}\left(\frac{\partial F}{\partial t}\right)

The derivatives are (I think)

\frac{\partial F}{\partial A_{3}}=\frac{F}{A_{3}}

\frac{\partial F}{\partial t}=\frac{F}{t}+A_{2}t\mathrm{plogis}^{\prime}\left(A_{2}t\mathrm{plogis}\left(A_{1}t^{2}\right)\right)\left(\mathrm{plogis}\left(A_{1}t^{2}\right)+2A_{1}t^{2}\mathrm{plogis}^{\prime}\left(A_{1}t^{2}\right)\right)

and of course at the solution F=1.

So, something like this

```
data {
int N;
vector[N] rt;
real median_A_3;
real<lower=0> sigma;
real alpha;
real<lower=0> beta;
}
parameters {
real<lower=0> A_2;
vector<lower=0,upper=1>[N] x;
}
transformed parameters {
vector[N] A_1 = alpha + x*beta;
vector[N] A_3;
for (i in 1:N) {
A_3[i] = 1.0/(rt[i]*plogis(A_2*rt[i]*plogis(A_1[i]*rt[i]*rt[i])));
}
}
model {
x ~ uniform(0,1);
//A_2 ~ prior() (if it's a parameter)
A_3 ~ lognormal(median_A_3, sigma);
// integrated likelihood, up to a constant
for (i in 1:N) {
real rt_ = rt[i];
real A1rt2 = A_1[i]*rt_*rt_;
real pA1rt2 = plogis(A1rt2);
real d = A_2*rt_*plogis_der(pA1rt2)*(pA1rt2 + 2*A1rt2*plogis_der(A1rt2));
target += log(A_3[i]) + log(1/rt_ + d);
}
}
```