Last updated: 2018-05-12

workflowr checks: (Click a bullet for more information)
  • R Markdown file: up-to-date

    Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

  • Repository version: 5b37742

    Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.

    Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
    
    Ignored files:
        Ignored:    .DS_Store
        Ignored:    .Rhistory
        Ignored:    .Rproj.user/
        Ignored:    analysis/.DS_Store
        Ignored:    analysis/BH_robustness_cache/
        Ignored:    analysis/FDR_Null_cache/
        Ignored:    analysis/FDR_null_betahat_cache/
        Ignored:    analysis/Rmosek_cache/
        Ignored:    analysis/StepDown_cache/
        Ignored:    analysis/alternative2_cache/
        Ignored:    analysis/alternative_cache/
        Ignored:    analysis/ash_gd_cache/
        Ignored:    analysis/average_cor_gtex_2_cache/
        Ignored:    analysis/average_cor_gtex_cache/
        Ignored:    analysis/brca_cache/
        Ignored:    analysis/cash_deconv_cache/
        Ignored:    analysis/cash_fdr_1_cache/
        Ignored:    analysis/cash_fdr_2_cache/
        Ignored:    analysis/cash_fdr_3_cache/
        Ignored:    analysis/cash_fdr_4_cache/
        Ignored:    analysis/cash_fdr_5_cache/
        Ignored:    analysis/cash_fdr_6_cache/
        Ignored:    analysis/cash_plots_cache/
        Ignored:    analysis/cash_plots_fdp_cache/
        Ignored:    analysis/cash_sim_1_cache/
        Ignored:    analysis/cash_sim_2_cache/
        Ignored:    analysis/cash_sim_3_cache/
        Ignored:    analysis/cash_sim_4_cache/
        Ignored:    analysis/cash_sim_5_cache/
        Ignored:    analysis/cash_sim_6_cache/
        Ignored:    analysis/cash_sim_7_cache/
        Ignored:    analysis/correlated_z_2_cache/
        Ignored:    analysis/correlated_z_3_cache/
        Ignored:    analysis/correlated_z_cache/
        Ignored:    analysis/create_null_cache/
        Ignored:    analysis/cutoff_null_cache/
        Ignored:    analysis/design_matrix_2_cache/
        Ignored:    analysis/design_matrix_cache/
        Ignored:    analysis/diagnostic_ash_cache/
        Ignored:    analysis/diagnostic_correlated_z_2_cache/
        Ignored:    analysis/diagnostic_correlated_z_3_cache/
        Ignored:    analysis/diagnostic_correlated_z_cache/
        Ignored:    analysis/diagnostic_plot_2_cache/
        Ignored:    analysis/diagnostic_plot_cache/
        Ignored:    analysis/efron_leukemia_cache/
        Ignored:    analysis/fitting_normal_cache/
        Ignored:    analysis/gaussian_derivatives_2_cache/
        Ignored:    analysis/gaussian_derivatives_3_cache/
        Ignored:    analysis/gaussian_derivatives_4_cache/
        Ignored:    analysis/gaussian_derivatives_5_cache/
        Ignored:    analysis/gaussian_derivatives_cache/
        Ignored:    analysis/gd-ash_cache/
        Ignored:    analysis/gd_delta_cache/
        Ignored:    analysis/gd_lik_2_cache/
        Ignored:    analysis/gd_lik_cache/
        Ignored:    analysis/gd_w_cache/
        Ignored:    analysis/knockoff_10_cache/
        Ignored:    analysis/knockoff_2_cache/
        Ignored:    analysis/knockoff_3_cache/
        Ignored:    analysis/knockoff_4_cache/
        Ignored:    analysis/knockoff_5_cache/
        Ignored:    analysis/knockoff_6_cache/
        Ignored:    analysis/knockoff_7_cache/
        Ignored:    analysis/knockoff_8_cache/
        Ignored:    analysis/knockoff_9_cache/
        Ignored:    analysis/knockoff_cache/
        Ignored:    analysis/knockoff_var_cache/
        Ignored:    analysis/marginal_z_alternative_cache/
        Ignored:    analysis/marginal_z_cache/
        Ignored:    analysis/mosek_reg_2_cache/
        Ignored:    analysis/mosek_reg_4_cache/
        Ignored:    analysis/mosek_reg_5_cache/
        Ignored:    analysis/mosek_reg_6_cache/
        Ignored:    analysis/mosek_reg_cache/
        Ignored:    analysis/pihat0_null_cache/
        Ignored:    analysis/plot_diagnostic_cache/
        Ignored:    analysis/poster_obayes17_cache/
        Ignored:    analysis/real_data_simulation_2_cache/
        Ignored:    analysis/real_data_simulation_3_cache/
        Ignored:    analysis/real_data_simulation_4_cache/
        Ignored:    analysis/real_data_simulation_5_cache/
        Ignored:    analysis/real_data_simulation_cache/
        Ignored:    analysis/rmosek_primal_dual_2_cache/
        Ignored:    analysis/rmosek_primal_dual_cache/
        Ignored:    analysis/seqgendiff_cache/
        Ignored:    analysis/simulated_correlated_null_2_cache/
        Ignored:    analysis/simulated_correlated_null_3_cache/
        Ignored:    analysis/simulated_correlated_null_cache/
        Ignored:    analysis/simulation_real_se_2_cache/
        Ignored:    analysis/simulation_real_se_cache/
        Ignored:    analysis/smemo_2_cache/
        Ignored:    data/LSI/
        Ignored:    docs/.DS_Store
        Ignored:    docs/figure/.DS_Store
        Ignored:    output/fig/
    
    Untracked files:
        Untracked:  analysis/cash_plots_fdp_files/
        Untracked:  docs/cash_plots_fdp_files/
    
    
    Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
Expand here to see past versions:
    File Version Author Date Message
    rmd cc0ab83 Lei Sun 2018-05-11 update
    html 0f36d99 LSun 2017-12-21 Build site.
    html 853a484 LSun 2017-11-07 Build site.
    html 043bf89 LSun 2017-11-05 transfer
    rmd 3c00d15 LSun 2017-04-21 de bogged
    html 3c00d15 LSun 2017-04-21 de bogged
    rmd 239f195 LSun 2017-04-21 GD-ASH
    html 239f195 LSun 2017-04-21 GD-ASH

Problem setting

We are considering the normal means problem with heteroskedastic noise. Furthermore, The noise are not independent, with unknown pairwise correlation \(\rho_{ij}\).

\[ \begin{array}{rcl} \hat\beta_j |\beta_j, \hat s_j &=& \beta_j + \hat s_j z_j \ ;\\ z_j &\sim& N(0, 1) \text{ marginally} \ ;\\ z_j &:& \text{correlated} \ . \end{array} \]

Under ash framework, a prior is put on exchangeable \(\beta_j\):

\[ g(\beta_j) = \sum_k \pi_k g_k(\beta_j)\ . \]

According to previous exploration, correlated standard normal \(z\) scores can be seen as iid from a density composed of Gaussian derivatives.

\[ f_0(z_j) = \sum_{l=0}^\infty w_l \varphi^{(l)}(z_j) \ . \] By change of variables, the likelihood

\[ p(\hat\beta_j | \beta_j, \hat s_j) = \frac{1}{\hat s_j}f_0\left(\frac{\hat\beta_j - \beta_j}{\hat s_j}\right) = \sum_l w_l \frac{1}{\hat s_j}\varphi^{(l)}\left(\frac{\hat\beta_j - \beta_j}{\hat s_j}\right) \ . \]

Thus for each pair of observation \((\hat\beta_j, \hat s_j)\), the likelihood becomes

\[ \begin{array}{rcl} f(\hat\beta_j|\hat s_j) &=& \displaystyle\int p(\hat\beta_j | \beta_j, \hat s_j)g(\beta_j)d\beta_j\\ &=&\displaystyle\int \sum_l w_l \frac{1}{\hat s_j}\varphi^{(l)}\left(\frac{\hat\beta_j - \beta_j}{\hat s_j}\right) \sum_k \pi_k g_k(\beta_j)d\beta_j\\ &=& \displaystyle \sum_k \sum_l \pi_k w_l \int\frac{1}{\hat s_j} \varphi^{(l)}\left(\frac{\hat\beta_j - \beta_j}{\hat s_j}\right) g_k(\beta_j)d\beta_j\\ &:=& \displaystyle \sum_k \sum_l \pi_k w_l f_{jkl} \ . \end{array} \] Hence, it’s all boiled down to calculate \(f_{jkl}\), which is the convolution of a Gaussian derivative \(\varphi^{(l)}\) and a prior component \(g_k\), with some change of variables manipulation. Note that in usual ASH settings, \(g_k\) is either uniform or normal, both of which can be handled without conceptual difficulty.

If \(g_k\) is a uniform, the convolution of a Gaussian derivative and a uniform is just another Gaussian derivative in a lower order, such as

\[ \varphi^{(l)} * \text{Unif} \leadsto \varphi^{(l-1)} \ . \] On the other hand, if \(g_k\) is a normal, we can use a fact about convolution that

\[ \begin{array}{rcl} && \frac{d}{dx}(f*g) = \frac{df}{dx} * g = f * \frac{dg}{dx}\\ &\Rightarrow& \varphi^{(l)} * N(\mu, \sigma^2) = \left(\varphi * N(\mu, \sigma^2)\right)^{(l)} \leadsto \tilde\varphi^{(l)} \ . \end{array} \] Because the convolution of \(\varphi\), a Gaussian, and another Gaussian is still a Gaussian, the convolution of a Gaussian derivative and a normal gives another Gaussian derivative.

With \(f_{jkl}\) computed, the goal is then to maximize the joint likeliood of the observation \(\left\{\left(\hat\beta_1, \hat s_1\right), \cdots,\left(\hat\beta_n, \hat s_n\right) \right\}\) which is \[ \max\limits_{\pi, w}\prod_j f(\hat\beta_j|\hat s_j) = \prod_j \left(\displaystyle \sum_k \sum_l \pi_k w_l f_{jkl}\right) \ , \] subject to reasonable, well designed constraints on \(\pi_k\) and especially \(w_l\).

Characteristic function / Fourier transform

The exact form of \(f_{jkl}\) can be derived analytically. Below is a method using characteristic functions or Fourier transforms.

Let \(\zeta_X(t) = \zeta_F(t) =\zeta_f(t) := E[e^{itX}]\) be the characteristic function of the random variable \(X\sim dF = f\). Either the random variable, its distribution, or its density function can be put in the subscript, depending on the circumstances.

On the other hand, the characteristic function of a random variable is also closely related to the Fourier transform of its density. In particular,

\[ \zeta_f(t) = E[e^{itX}] = \int e^{itx}dF(x) = \int e^{itx}f(x)dx := \mathscr{F}_f(t) \ . \] Note that in usual definition, the Fourier transform is given by \[ \hat f(\xi) := \int f(x)e^{-2\pi ix\xi}dx \ , \] with a normalizing factor \(2\pi\) and a negative sign. However, even defined in different ways, \(\mathscr{F}_f(t)\) carries over many nice properties of \(\hat f(\xi)\), for example, as we’ll show,

\[ \mathscr{F}_{f^{(m)}}(t) = (-it)^m\mathscr{F}_f(t) \ , \]

where \(f^{(m)}\) is the \(m^\text{th}\) derivative of \(f\).

Also, \[ \mathscr{F}_{f*g}(t) = \mathscr{F}_f(t)\mathscr{F}_g(t) \ , \] where \(*\) stands for convolution.

With these properties, the Fourier transform tool could be very useful when dealing with Gaussian derivatives and their convolution with other (density) functions.

Under this definition, the inversion formula for the characteristic functions, also known as the inverse Fourier transform, is

\[ f(x) = \frac1{2\pi} \int e^{-itx}\zeta_f(t)dt = \frac{1}{2\pi}\int e^{-itx}\mathscr{F}_f(t)dt \ . \]

Then the characteristic function of \(\hat\beta_j | \hat s_j\)

\[ \zeta_{\hat\beta_j|\hat s_j}(t) = E\left[e^{it\hat\beta_j}\mid\hat s_j\right] = E\left[e^{it(\beta_j + \hat s_j z_j)}\mid\hat s_j\right] = E\left[e^{it\beta_j}\right]E\left[e^{it\hat s_jz_j}\mid\hat s_j\right] = \zeta_{\beta}(t)\zeta_{f_0}\left(\hat s_jt\right) \ . \]

Let’s take care of \(\zeta_{\beta}(t)\) and \(\zeta_{f_0}\left(\hat s_jt\right)\) one by one. Note that

\[ \begin{array}{rrcl} &\beta_j &\sim& \sum_k\pi_kg_k \\ \Rightarrow & \zeta_{\beta}(t) &=& \displaystyle\int e^{it\beta_j}p(\beta_j)d\beta_j = \int e^{it\beta_j}\sum_k\pi_kg_k(\beta_j)d\beta_j =\sum_k\pi_k\int e^{it\beta_j}g_k(\beta_j)d\beta_j \\ &&=& \sum_k\pi_k\zeta_{g_k}(t) \ . \end{array} \]

Meanwhile, using the fact that Gaussian derivatives should be absolutely integrable,

\[ \begin{array}{rrcl} &f_0(z_j) &=& \sum_l w_l \varphi^{(l)}(z_j) \\ \Rightarrow & \zeta_{f_0}(t) &=& \displaystyle \int e^{itz_j}f_0(z_j)dz_j = \int e^{itz_j}\sum_lw_l\varphi^{(l)}(z_j)dz_j = \sum_lw_l\int e^{itz_j}\varphi^{(l)}(z_j)dz_j\\ &&=& \sum_l w_l \mathscr{F}_{\varphi^{(l)}}(t) \ . \end{array} \] Each Fourier transform of gaussian derivatives

\[ \begin{array}{rcl} \mathscr{F}_{\varphi^{(l)}}(t) &=& \displaystyle\int e^{itz_j}\varphi^{(l)}(z_j)dz_j\\ &=& \displaystyle e^{itz_j}\varphi^{(l-1)}(z_j)|_{-\infty}^\infty - \int it e^{itz_j}\varphi^{(l-1)}(z_j)dz_j\\ &=&\displaystyle 0-it\int e^{itz_j}\varphi^{(l-1)}(z_j)dz_j\\ &=&\displaystyle -it\int e^{itz_j}\varphi^{(l-1)}(z_j)dz_j\\ &=&\displaystyle \cdots\\ &=&\displaystyle (-it)^l \int e^{itz_j}\varphi(z_j)dz_j\\ &=&(-it)^l \mathscr{F}_{\varphi}(t)\\ &=&(-it)^l \zeta_\varphi(t)\\ &=&\displaystyle (-it)^l e^{-\frac12t^2} \ . \end{array} \] Thus,

\[ \begin{array}{rrcl} &\zeta_{\hat\beta_j|\hat s_j}(t) &=&\zeta_{\beta}(t)\zeta_{f_0}(\hat s_jt)\\ &&=&\left(\sum_k\pi_k\zeta_{g_k}(t)\right) \left(\sum_lw_l(-i\hat s_jt)^l \zeta_\varphi(\hat s_jt)\right)\\ &&=&\sum_k\sum_l\pi_kw_l(-i\hat s_jt)^l \zeta_{g_k}(t) \zeta_\varphi(\hat s_jt)\\ &&=&\sum_k\sum_l\pi_kw_l(-i\hat s_jt)^l \zeta_{g_k}(t) e^{-\frac12\hat s_j^2t^2}\\ \Rightarrow & f(\hat\beta_j|\hat s_j) &=&\displaystyle\frac1{2\pi}\int e^{-it\hat\beta_j}\zeta_{\hat\beta_j|\hat s_j}(t)dt\\ &&=&\displaystyle\frac1{2\pi}\int e^{-it\hat\beta_j}\left(\sum_k\sum_l\pi_kw_l(-i\hat s_jt)^l \zeta_{g_k}(t) \zeta_\varphi(\hat s_jt)\right)dt\\ &&=&\displaystyle\sum_k\sum_l\pi_kw_l \frac1{2\pi}\int e^{-it\hat\beta_j}(-i\hat s_jt)^l \zeta_{g_k}(t) \zeta_\varphi(\hat s_jt)dt\\ &&:=& \sum_k\sum_l\pi_kw_lf_{jkl} \ . \end{array} \] It is essentially the equivalent expression of \(f_{jkl}\) as

\[ \begin{array}{rcl} f_{jkl} &=& \displaystyle \int\frac{1}{\hat s_j} \varphi^{(l)}\left(\frac{\hat\beta_j - \beta_j}{\hat s_j}\right) g_k(\beta_j)d\beta_j\\ &=& \displaystyle\frac1{2\pi}\int e^{-it\hat\beta_j}(-i\hat s_jt)^l \zeta_{g_k}(t) \zeta_\varphi(\hat s_jt)dt \ . \end{array} \]

Uniform mixture prior

If the prior of \(\beta\) is a mixture of uniforms, \(\beta \sim g = \sum_k\pi_kg_k\) where \(g_k = \text{Unif}\left[a_k, b_k\right]\),

\[ \begin{array}{rrcl} &\zeta_{g_k}(t) &=& \displaystyle\frac{e^{itb_k} - e^{ita_k}}{it(b_k - a_k)}\\ \Rightarrow &\zeta_{g_k}(t)\zeta_\varphi(\hat s_jt) &=& \displaystyle\frac{e^{itb_k} - e^{ita_k}}{it(b_k - a_k)} e^{-\frac12\hat s_j^2t^2}\\ &&=& \displaystyle\frac{e^{itb_k-\frac12\hat s_j^2t^2}}{it(b_k - a_k)} -\displaystyle\frac{e^{ita_k-\frac12\hat s_j^2t^2}}{it(b_k - a_k)}\\ \Rightarrow & f_{jkl} &=& \displaystyle\frac1{2\pi}\int e^{-it\hat\beta_j}(-i\hat s_jt)^l \zeta_{g_k}(t) \zeta_\varphi(\hat s_jt)dt\\ &&=& \displaystyle\frac1{2\pi}\int e^{-it\hat\beta_j}(-i\hat s_jt)^l\displaystyle\frac{e^{itb_k-\frac12\hat s_j^2t^2}}{it(b_k - a_k)}dt - \displaystyle\frac1{2\pi}\int e^{-it\hat\beta_j}(-i\hat s_jt)^l\displaystyle\frac{e^{ita_k-\frac12\hat s_j^2t^2}}{it(b_k - a_k)}dt \ . \end{array} \]

It looks pretty complicated, but if we take advantage of usual Fourier transform tricks things get clearer.

Let \(\varphi_{\mu, \sigma^2}\) denote the density function of \(N(\mu, \sigma^2)\),

\[ \begin{array}{rl} &\varphi_{\mu, \sigma^2}(z) = \frac1\sigma\varphi\left(\frac{z - \mu}{\sigma}\right)\\ \Rightarrow & \varphi_{\mu, \sigma^2}^{(m)}(z) = \frac1{\sigma^{m+1}}\varphi^{(m)}\left(\frac{z - \mu}{\sigma}\right) \\ \Rightarrow & \zeta_{\varphi_{\mu, \sigma^2}}(t) = e^{it\mu - \frac12\sigma^2t^2} \ . \end{array} \]

Take the first part of \(f_{jkl}\), \[ \begin{array}{rrcl} & (-i\hat s_jt)^l\displaystyle\frac{e^{itb_k-\frac12\hat s_j^2t^2}}{it(b_k - a_k)} &=&\displaystyle -\frac{\hat s_j^l}{b_k - a_k}(-it)^{l-1}\zeta_{\varphi_{b_k, \hat s_j^2}}(t)\\ &&=&\displaystyle -\frac{\hat s_j^l}{b_k - a_k}(-it)^{l-1}\mathscr{F}_{\varphi_{b_k, \hat s_j^2}}(t)\\ &&=&\displaystyle -\frac{\hat s_j^l}{b_k - a_k}\mathscr{F}_{\varphi_{b_k, \hat s_j^2}^{(l-1)}}(t)\\ \Rightarrow & \displaystyle\frac1{2\pi}\int e^{-it\hat\beta_j}(-i\hat s_jt)^l\displaystyle\frac{e^{itb_k-\frac12\hat s_j^2t^2}}{it(b_k - a_k)}dt &=& -\displaystyle\frac1{2\pi}\int e^{-it\hat\beta_j} \frac{\hat s_j^l}{b_k - a_k}\mathscr{F}_{\varphi_{b_k, \hat s_j^2}^{(l-1)}}(t)dt \\ &&=&\displaystyle -\frac{\hat s_j^l}{b_k - a_k} \frac1{2\pi} \int e^{-it\hat\beta_j} \mathscr{F}_{\varphi_{b_k, \hat s_j^2}^{(l-1)}}(t) dt\\ &&=&\displaystyle -\frac{\hat s_j^l}{b_k - a_k} \varphi_{b_k, \hat s_j^2}^{(l-1)}\left(\hat\beta_j\right)\\ &&=&\displaystyle -\frac{\hat s_j^l}{b_k - a_k} \frac1{\hat s_j^{l}}\varphi^{(l-1)}\left(\frac{\hat\beta_j- b_k}{\hat s_j}\right)\\ &&=&\displaystyle - \frac{\varphi^{(l-1)}\left(\frac{\hat\beta_j- b_k}{\hat s_j}\right)}{b_k-a_k} \ . \end{array} \] Therefore,

\[ \begin{array}{rcl} f_{jkl} &=& \displaystyle \frac1{2\pi}\int e^{-it\hat\beta_j}(-i\hat s_jt)^l\displaystyle\frac{e^{itb_k-\frac12\hat s_j^2t^2}}{it(b_k - a_k)}dt - \frac1{2\pi}\int e^{-it\hat\beta_j}(-i\hat s_jt)^l\displaystyle\frac{e^{ita_k-\frac12\hat s_j^2t^2}}{it(b_k - a_k)}dt\\ &=& \displaystyle \left(- \frac{\varphi^{(l-1)}\left(\frac{\hat\beta_j- b_k}{\hat s_j}\right)}{b_k-a_k}\right) - \left(- \frac{\varphi^{(l-1)}\left(\frac{\hat\beta_j- a_k}{\hat s_j}\right)}{b_k-a_k}\right)\\ &=&\displaystyle \frac{\varphi^{(l-1)}\left(\frac{\hat\beta_j- a_k}{\hat s_j}\right) - \varphi^{(l-1)}\left(\frac{\hat\beta_j- b_k}{\hat s_j}\right)}{b_k-a_k} \end{array} \]

Normal mixture prior

Likewise, if the prior of \(\beta\) is a mixture of normals, \(\beta \sim g = \sum_k\pi_kg_k\) where \(g_k = N(\mu_k, \sigma_k^2)\),

\[ \begin{array}{rrcl} &\zeta_{g_k}(t) &=& \displaystyle e^{it\mu_k - \frac12\sigma_k^2t^2}\\ \Rightarrow &\zeta_{g_k}(t)\zeta_\varphi(\hat s_jt) &=& e^{it\mu_k - \frac12\sigma_k^2t^2} e^{-\frac12\hat s_j^2t^2}\\ &&=& e^{it\mu_k - \frac12\left(\sigma_k^2+\hat s_j^2\right)t^2} \\ &&=& \zeta_{\varphi_{\mu_k, \sigma_k^2 + \hat s_j^2}}(t)\\ &&=& \mathscr{F}_{\varphi_{\mu_k, \sigma_k^2 + \hat s_j^2}}(t) \\ \Rightarrow & (-i\hat s_jt)^l\zeta_{g_k}(t)\zeta_\varphi(\hat s_jt) &=& (-i\hat s_jt)^l\mathscr{F}_{\varphi_{\mu_k, \sigma_k^2 + \hat s_j^2}}(t)\\ &&=& \hat s_j^l(-it)^l\mathscr{F}_{\varphi_{\mu_k, \sigma_k^2 + \hat s_j^2}}(t) \\ &&=& \hat s_j^l \mathscr{F}_{\varphi_{\mu_k, \sigma_k^2 + \hat s_j^2}^{(l)}}(t) \\ \Rightarrow & f_{jkl} &=& \displaystyle\frac1{2\pi}\int e^{-it\hat\beta_j}(-i\hat s_jt)^l \zeta_{g_k}(t) \zeta_\varphi(\hat s_jt)dt\\ &&=& \displaystyle\frac1{2\pi}\int e^{-it\hat\beta_j} \hat s_j^l \mathscr{F}_{\varphi_{\mu_k, \sigma_k^2 + \hat s_j^2}^{(l)}}(t)dt\\ &&=& \hat s_j^l \displaystyle\frac1{2\pi}\int e^{-it\hat\beta_j} \mathscr{F}_{\varphi_{\mu_k, \sigma_k^2 + \hat s_j^2}^{(l)}}(t)dt\\ &&=& \hat s_j^l \varphi_{\mu_k, \sigma_k^2 + \hat s_j^2}^{(l)}\left(\hat\beta_j\right)\\ &&=& \displaystyle \hat s_j^l \frac{1}{\left(\sqrt{\sigma_k^2 + \hat s_j^2}\right)^{l+1}} \varphi^{(l)}\left(\frac{\hat\beta_j - \mu_k}{\sqrt{\sigma_k^2 + \hat s_j^2}}\right)\\ &&=& \displaystyle \frac{\hat s_j^l}{\left(\sqrt{\sigma_k^2 + \hat s_j^2}\right)^{l+1}} \varphi^{(l)}\left(\frac{\hat\beta_j - \mu_k}{\sqrt{\sigma_k^2 + \hat s_j^2}}\right) \ . \end{array} \] In common applications, \(\mu_k\equiv\mu\equiv0\).

Optimization problem

The problem is now boiled down to maximize the joint likelihood

\[ \begin{array}{rcl} \max\limits_{\pi, w}\prod\limits_j f(\hat\beta_j|\hat s_j) &=& \max\limits_{\pi, w}\prod\limits_j \left(\sum_k\sum_l\pi_k w_l f_{jkl}\right)\\ &\Leftrightarrow& \max\limits_{\pi, w}\sum_j\log\left(\sum_k\sum_l\pi_k w_l f_{jkl}\right) \ , \end{array} \]

subject to appropriate constraints on \(\pi_k\) and especially \(w_l\), where the specific form of \(f_{jkl}\) depends on the mixture component of the prior \(g\) of \(\beta_j\). Here we consider two cases.

Uniform mixture prior

\[ \begin{array}{rrcl} &\beta_j &\sim & \sum_k \pi_k \text{ Unif }[a_k, b_k]\\ \Rightarrow & f_{jkl} &= & \displaystyle\frac{\varphi^{(l-1)}\left(\frac{\hat\beta_j-a_k}{\hat s_j}\right) - \varphi^{(l-1)}\left(\frac{\hat\beta_j-b_k}{\hat s_j}\right)}{b_k - a_k} \ . \end{array} \]

Normal mixture prior

\[ \begin{array}{rrcl} &\beta_j &\sim & \sum_k \pi_k N\left(\mu_k, \sigma_k^2\right) \\ \Rightarrow & f_{jkl} &= & \displaystyle\frac{\hat s_j^l}{\left(\sqrt{\sigma_k^2 + \hat s_j^2}\right)^{l+1}} \varphi^{(l)}\left(\frac{ \hat\beta_j - \mu_k }{ \sqrt{\sigma_k^2 + \hat s_j^2} }\right) \ . \end{array} \]

Biconvex optimization

Our goal is then to estimate \(\hat\pi\) by solving the following constrained biconvex optimization problem

\[ \begin{array}{rl} \max\limits_{\pi,w} & \sum_j\log\left(\sum_k\sum_l\pi_k w_l f_{jkl}\right)\\ \text{subject to} & \sum_k\pi_k = 1\\ & w_0 = 1\\ & \sum_l w_l \varphi^{l}(z) \geq 0, \forall z\in \mathbb{R}\\ & w_l \text{ decay reasonably fast.} \end{array} \]




This reproducible R Markdown analysis was created with workflowr 1.0.1