In a micro step-ahead, the Metropolis-within-Gibbs code can now start from a previous run. This is advantageous for doing the computation in batches and testing whether the chain is already converged or we have enough samples for a nice posterior distribution plot.
I have also tested the more elaborate hierarchical Bayesian analysis of coronal loops oscillations. In one of our previous paper (Arregui & Asensio Ramos 2011), we applied a simple Bayesian analysis of single coronal loops oscillations. The measured period ($P$) and damping time ($\tau_d$) were related, in our model, with the Alfvén time $\tau_A$, density contrast between the interior and exterior of the coronal loop $\xi$, and the ratio between the length and radius of the tube $l/R$ $$P=\tau_A \sqrt{2} \sqrt{\frac{\xi+1}{\xi}}$$ and $$\frac{\tau_d}{P} = \frac{2}{\pi} \frac{\xi+1}{\xi-1} \frac{1}{l/R}$$ Under the assumption that $P$ and $\tau_d$ are measured quantities, we applied a MCMC code to sample the posterior distributions. We observed a quite strong prior sensitivity for some of the parameters, but reliable inferences for others.
The idea now is to go much deeper in the analysis. On the one hand, we plan to get closer to the "measurement" phase. The period and the damping time are inferred quantities from the analysis of the fluctuations of coronal loops obtained from satellite images of the corona. They are obtained after assuming that the oscillation (measured in pixels vs. time from the initial time) is given by the following complex approximate formula:
$$x(t)=x_\mathrm{osc}(t)+x_\mathrm{trend}(t)$$
The trend absorbs anything that is not purely oscillatory and is given by the following polynomial fit:
$$x_\mathrm{trend}(t)=a_0+a_1 t + a_2 t^2+a_3 t^3+ \cdots$$
The oscillatory part is modeled as a damped oscillation:
$$x_\mathrm{osc}(t)=A_0 \sin \left( \frac{2 \pi}{P} t + \phi \right) \exp \left( - \frac{t}{\tau_d} \right)$$
Consequently, assuming the previous expressions for the period and damping time, the i-th loop is described by the set of parameters:
$$\{(A_0)_i, (\tau_A)_i, \xi_i, (l/R)_i, \phi_i, \mathbf{a}_i\}$$
Depending on the complexity of the trend, we can have between 7 and 11 parameters to describe not more than several dozen observed points.
Our idea is then to consider that these parameters are realizations of some parametric distributions (thus the hierarchical character) and try to learn the parameters of these priors. We have proposed the following priors that we are now trying to learn from the data:
$$(\tau_A)_i \sim \mathrm{IG}(\gamma,\delta)$$
$$(\xi)_i \sim \mathrm{IG}(\epsilon,\eta,1)$$
$$(l/R)_i \sim \mathrm{Beta}(\alpha,\beta,0,2)$$
$$(A_0)_i \sim \mathrm{Gamma}(0.01,0.01)$$
$$(a_j)_i \sim \mathrm{N}(0,\sigma^2)$$
with the data following a normal distribution around the previous oscillatory+trend model. We consider that the previous distributions are sufficiently general and we have to play now with the hyperpriors for the parameters $\alpha$, $\beta$, $\gamma$, $\delta$, $\epsilon$ and $\eta$.
I have also tested the more elaborate hierarchical Bayesian analysis of coronal loops oscillations. In one of our previous paper (Arregui & Asensio Ramos 2011), we applied a simple Bayesian analysis of single coronal loops oscillations. The measured period ($P$) and damping time ($\tau_d$) were related, in our model, with the Alfvén time $\tau_A$, density contrast between the interior and exterior of the coronal loop $\xi$, and the ratio between the length and radius of the tube $l/R$ $$P=\tau_A \sqrt{2} \sqrt{\frac{\xi+1}{\xi}}$$ and $$\frac{\tau_d}{P} = \frac{2}{\pi} \frac{\xi+1}{\xi-1} \frac{1}{l/R}$$ Under the assumption that $P$ and $\tau_d$ are measured quantities, we applied a MCMC code to sample the posterior distributions. We observed a quite strong prior sensitivity for some of the parameters, but reliable inferences for others.
The idea now is to go much deeper in the analysis. On the one hand, we plan to get closer to the "measurement" phase. The period and the damping time are inferred quantities from the analysis of the fluctuations of coronal loops obtained from satellite images of the corona. They are obtained after assuming that the oscillation (measured in pixels vs. time from the initial time) is given by the following complex approximate formula:
$$x(t)=x_\mathrm{osc}(t)+x_\mathrm{trend}(t)$$
The trend absorbs anything that is not purely oscillatory and is given by the following polynomial fit:
$$x_\mathrm{trend}(t)=a_0+a_1 t + a_2 t^2+a_3 t^3+ \cdots$$
The oscillatory part is modeled as a damped oscillation:
$$x_\mathrm{osc}(t)=A_0 \sin \left( \frac{2 \pi}{P} t + \phi \right) \exp \left( - \frac{t}{\tau_d} \right)$$
Consequently, assuming the previous expressions for the period and damping time, the i-th loop is described by the set of parameters:
$$\{(A_0)_i, (\tau_A)_i, \xi_i, (l/R)_i, \phi_i, \mathbf{a}_i\}$$
Depending on the complexity of the trend, we can have between 7 and 11 parameters to describe not more than several dozen observed points.
Our idea is then to consider that these parameters are realizations of some parametric distributions (thus the hierarchical character) and try to learn the parameters of these priors. We have proposed the following priors that we are now trying to learn from the data:
$$(\tau_A)_i \sim \mathrm{IG}(\gamma,\delta)$$
$$(\xi)_i \sim \mathrm{IG}(\epsilon,\eta,1)$$
$$(l/R)_i \sim \mathrm{Beta}(\alpha,\beta,0,2)$$
$$(A_0)_i \sim \mathrm{Gamma}(0.01,0.01)$$
$$(a_j)_i \sim \mathrm{N}(0,\sigma^2)$$
with the data following a normal distribution around the previous oscillatory+trend model. We consider that the previous distributions are sufficiently general and we have to play now with the hyperpriors for the parameters $\alpha$, $\beta$, $\gamma$, $\delta$, $\epsilon$ and $\eta$.
No comments:
Post a Comment