Let me integrate you out so that no one would notice you
An old fashioned mathematical technique might be to integrate out parameters from equations, and as more the equation is complicated as more it is difficult to get the analytical solution.
Basic idea
Consider a variable exchange $x \mapsto z$ that satisfies \begin{align} \int_{-\infty}^{\infty} f(x) dx \Longrightarrow \int_{-1}^{1} f\Bigl(x(z)\Bigr) x'(z) dz. \nonumber \end{align} An easy solution for that is \begin{align} z = 2 \left(\frac{1}{1+\exp(-x)} \right)-1 \Longleftrightarrow x = \log{\left ( \frac{1+z}{1-z} \right )} \nonumber \end{align} and its first order derivative is \begin{align} x'(z) = \frac{2}{1 - z^{2}} \nonumber \end{align} It is well known that \begin{align} \int_{-\infty}^{\infty} \frac{1}{\sqrt{2 \pi}} \exp\left(-\frac{x^2}{2} \right) dx = 1, \nonumber \end{align} and apply the variable exchange and the number we get by approximating 7-point Gauss rule with a 15-point Kronrod rule is: \begin{align} \int_{-1}^{1} \frac{1}{\sqrt{2 \pi}} \exp\left(-\frac{\left(\log{\left ( \frac{1+z}{1-z} \right )}\right)^2}{2} \right) \left|\frac{2}{1 - z^{2}}\right| dz = 0.99999692068548 \approx 1 \nonumber \end{align} Note that $1-0.99999692068548=3.07931451959398 \times 10^{-6}$.Test model
\begin{align} y = a_i x + \epsilon_i \\ a_i \sim {\cal N}(a,\omega) \\ \epsilon_i \sim {\cal N}(0,\sigma) \nonumber \end{align} Likelihood is \begin{align} \prod_{i=1}^n{\cal N}_d(y \mid a_i x, \sigma){\cal N}_d(a_i \mid a x,\omega) \nonumber \end{align}Analytical solution
Integrate out $a_i$ from the model \begin{align} p_i &= \int_{-\infty}^{\infty}{\cal N}_d(y \mid a_i x, \sigma){\cal N}_d(a_i \mid a x,\omega) d a_i \\ &= \int_{-\infty}^{\infty} \left[ \frac{1}{\sqrt{2 \pi} \sigma}\exp \left(-\frac{y - a_i x}{2 \sigma^2} \right) \frac{1}{\sqrt{2 \pi} \omega}\exp \left(-\frac{a_i - a}{2 \omega^2} \right) \right]d a_i \\ &=\frac{1}{\sqrt{2 \pi (\omega^2 x^2 + \sigma^2)}} \exp \left( -\frac{(y-ax)^2}{2(\omega^2 x^2 + \sigma^2)} \right) \nonumber \end{align}
> summary(fit1)$summary
mean se_mean sd 2.5% 25% 50%
a 0.9888168 0.0006589843 0.03325222 0.9221687 0.9669033 0.9888208
omega 0.4188159 0.0021747104 0.07486459 0.2419583 0.3795231 0.4274141
sigma 0.5302711 0.0005498567 0.01996546 0.4930311 0.5166785 0.5294182
lp__ -530.7041447 0.0428310403 1.40084589 -534.4240874 -531.3030933 -530.2699704
75% 97.5% n_eff Rhat
a 1.0112437 1.0524001 2546.192 1.002017
omega 0.4695590 0.5430313 1185.086 1.000445
sigma 0.5427861 0.5725930 1318.438 1.002972
lp__ -529.6919616 -529.2135198 1069.704 1.002579
Numerical solution
Compare the results with one using the following numerical integral. \begin{align} \int_{-\infty}^{\infty} f(x) dx \Longrightarrow \int_{-1}^{1} f\left(\log{\left ( \frac{1+z}{1-z} \right )}\right) \left|\frac{2}{1 - z^{2}}\right| dz. \nonumber \end{align}
> summary(fit2)$summary
mean se_mean sd 2.5% 25% 50%
a 0.9895701 0.0007101089 0.03497257 0.9192264 0.9664499 0.9899563
omega 0.4267222 0.0013115691 0.06287198 0.3079255 0.3800736 0.4278155
sigma 0.5293555 0.0004180750 0.01877470 0.4938143 0.5165164 0.5290596
lp__ -877.2053140 0.0278415504 1.21027072 -880.2556886 -877.7230284 -876.8940771
75% 97.5% n_eff Rhat
a 1.0127247 1.0576070 2425.521 1.000030
omega 0.4714495 0.5469434 2297.904 1.000609
sigma 0.5421843 0.5668702 2016.683 1.000074
lp__ -876.3328600 -875.8352797 1889.636 1.000688
Script
Data generation
set.seed(1234)
n = 1000
a = 1
omega = 0.5
sigma = 0.5
x = runif(n,-1,1)
a_i = rnorm(n,a,omega)
e_i = rnorm(n,0,sigma)
y = x*a_i + e_i
## 7-point Gauss rule with a 15-point Kronrod rule
X = c(
-0.991455371120813,
-0.949107912342759,
-0.864864423359769,
-0.741531185599394,
-0.586087235467691,
-0.405845151377397,
-0.207784955007898,
0,
0.207784955007898,
0.405845151377397,
0.586087235467691,
0.741531185599394,
0.864864423359769,
0.949107912342759,
0.991455371120813
)
C = c(
0.022935322010529,
0.063092092629979,
0.104790010322250,
0.140653259715525,
0.169004726639267,
0.190350578064785,
0.204432940075298,
0.209482141084728,
0.204432940075298,
0.190350578064785,
0.169004726639267,
0.140653259715525,
0.104790010322250,
0.063092092629979,
0.022935322010529
)
d = list(
n = n,
x = x,
y = y,
X = X,
C = C
)
plot(x,y)
Model for analytical solution
model1 = '
functions{
real ll(real a, real sigma, real omega, real x, real y){
real res;
res = -0.5*log(pi())-0.5*(log(omega^2*x^2+sigma^2)) - 0.5*(y - a*x)^2/(omega^2*x^2+sigma^2);
return res;
}
}
data{
int n;
real x[n];
real y[n];
}
parameters{
real a;
real omega;
real sigma;
}
model{
for (i in 1:n)
target += ll(a,sigma,omega,x[i],y[i]);
}
'
StanModel1 <- stan_model(model_code = model1)
fit1 <- sampling(StanModel1
,data = d
)
Model for numerical integral
model2 = '
functions{
real pdf(real exz, real z, real a, real sigma, real omega, real x, real y){
real res;
res = 1/(sqrt(2*pi())*sigma)*exp(-(y-exz*x)^2/(2*sigma^2))*1/(sqrt(2*pi())*omega)*exp(-(exz-a)^2/(2*omega^2))*(2/(1-z^2));
return res;
}
real int_pdf(vector X, vector C, real a, real sigma, real omega, real x, real y){
real res;
vector[15] Z;
for (i in 1:15)
Z[i] = log((1+X[i])/(1-X[i]));
res = 0;
for (i in 1:15)
res += C[i]*pdf(Z[i],X[i],a,sigma,omega,x,y);
return log(res);
}
}
data{
int n;
real x[n];
real y[n];
vector[15] X;
vector[15] C;
}
parameters{
real a;
real omega;
real sigma;
}
model{
for (i in 1:n)
target += int_pdf(X,C,a,sigma,omega,x[i],y[i]);
}
'
StanModel2 <- stan_model(model_code = model2)
fit2 <- sampling(StanModel2
,data = d
)
Comments
Post a Comment