<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta charset="utf-8" />
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="pandoc" />
<meta name="author" content="Lei Sun" />
<meta name="date" content="2017-04-12" />
<title>Normal Means with Heteroskedastic and Correlated Noise</title>
<script src="site_libs/jquery-1.11.3/jquery.min.js"></script>
<meta name="viewport" content="width=device-width, initial-scale=1" />
<link href="site_libs/bootstrap-3.3.5/css/cosmo.min.css" rel="stylesheet" />
<script src="site_libs/bootstrap-3.3.5/js/bootstrap.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/html5shiv.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/respond.min.js"></script>
<script src="site_libs/jqueryui-1.11.4/jquery-ui.min.js"></script>
<link href="site_libs/tocify-1.9.1/jquery.tocify.css" rel="stylesheet" />
<script src="site_libs/tocify-1.9.1/jquery.tocify.js"></script>
<script src="site_libs/navigation-1.1/tabsets.js"></script>
<link href="site_libs/highlightjs-9.12.0/textmate.css" rel="stylesheet" />
<script src="site_libs/highlightjs-9.12.0/highlight.js"></script>
<link href="site_libs/font-awesome-4.5.0/css/font-awesome.min.css" rel="stylesheet" />
<style type="text/css">code{white-space: pre;}</style>
<style type="text/css">
pre:not([class]) {
background-color: white;
}
</style>
<script type="text/javascript">
if (window.hljs) {
hljs.configure({languages: []});
hljs.initHighlightingOnLoad();
if (document.readyState && document.readyState === "complete") {
window.setTimeout(function() { hljs.initHighlighting(); }, 0);
}
}
</script>
<style type="text/css">
h1 {
font-size: 34px;
}
h1.title {
font-size: 38px;
}
h2 {
font-size: 30px;
}
h3 {
font-size: 24px;
}
h4 {
font-size: 18px;
}
h5 {
font-size: 16px;
}
h6 {
font-size: 12px;
}
.table th:not([align]) {
text-align: left;
}
</style>
</head>
<body>
<style type = "text/css">
.main-container {
max-width: 940px;
margin-left: auto;
margin-right: auto;
}
code {
color: inherit;
background-color: rgba(0, 0, 0, 0.04);
}
img {
max-width:100%;
height: auto;
}
.tabbed-pane {
padding-top: 12px;
}
button.code-folding-btn:focus {
outline: none;
}
</style>
<style type="text/css">
/* padding for bootstrap navbar */
body {
padding-top: 51px;
padding-bottom: 40px;
}
/* offset scroll position for anchor links (for fixed navbar) */
.section h1 {
padding-top: 56px;
margin-top: -56px;
}
.section h2 {
padding-top: 56px;
margin-top: -56px;
}
.section h3 {
padding-top: 56px;
margin-top: -56px;
}
.section h4 {
padding-top: 56px;
margin-top: -56px;
}
.section h5 {
padding-top: 56px;
margin-top: -56px;
}
.section h6 {
padding-top: 56px;
margin-top: -56px;
}
</style>
<script>
// manage active state of menu based on current page
$(document).ready(function () {
// active menu anchor
href = window.location.pathname
href = href.substr(href.lastIndexOf('/') + 1)
if (href === "")
href = "index.html";
var menuAnchor = $('a[href="' + href + '"]');
// mark it active
menuAnchor.parent().addClass('active');
// if it's got a parent navbar menu mark it active as well
menuAnchor.closest('li.dropdown').addClass('active');
});
</script>
<div class="container-fluid main-container">
<!-- tabsets -->
<script>
$(document).ready(function () {
window.buildTabsets("TOC");
});
</script>
<!-- code folding -->
<script>
$(document).ready(function () {
// move toc-ignore selectors from section div to header
$('div.section.toc-ignore')
.removeClass('toc-ignore')
.children('h1,h2,h3,h4,h5').addClass('toc-ignore');
// establish options
var options = {
selectors: "h1,h2,h3",
theme: "bootstrap3",
context: '.toc-content',
hashGenerator: function (text) {
return text.replace(/[.\\/?&!#<>]/g, '').replace(/\s/g, '_').toLowerCase();
},
ignoreSelector: ".toc-ignore",
scrollTo: 0
};
options.showAndHide = true;
options.smoothScroll = true;
// tocify
var toc = $("#TOC").tocify(options).data("toc-tocify");
});
</script>
<style type="text/css">
#TOC {
margin: 25px 0px 20px 0px;
}
@media (max-width: 768px) {
#TOC {
position: relative;
width: 100%;
}
}
.toc-content {
padding-left: 30px;
padding-right: 40px;
}
div.main-container {
max-width: 1200px;
}
div.tocify {
width: 20%;
max-width: 260px;
max-height: 85%;
}
@media (min-width: 768px) and (max-width: 991px) {
div.tocify {
width: 25%;
}
}
@media (max-width: 767px) {
div.tocify {
width: 100%;
max-width: none;
}
}
.tocify ul, .tocify li {
line-height: 20px;
}
.tocify-subheader .tocify-item {
font-size: 0.90em;
padding-left: 25px;
text-indent: 0;
}
.tocify .list-group-item {
border-radius: 0px;
}
</style>
<!-- setup 3col/9col grid for toc_float and main content -->
<div class="row-fluid">
<div class="col-xs-12 col-sm-4 col-md-3">
<div id="TOC" class="tocify">
</div>
</div>
<div class="toc-content col-xs-12 col-sm-8 col-md-9">
<div class="navbar navbar-default navbar-fixed-top" role="navigation">
<div class="container">
<div class="navbar-header">
<button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar">
<span class="icon-bar"></span>
<span class="icon-bar"></span>
<span class="icon-bar"></span>
</button>
<a class="navbar-brand" href="index.html">truncash</a>
</div>
<div id="navbar" class="navbar-collapse collapse">
<ul class="nav navbar-nav">
<li>
<a href="index.html">Home</a>
</li>
<li>
<a href="about.html">About</a>
</li>
<li>
<a href="license.html">License</a>
</li>
</ul>
<ul class="nav navbar-nav navbar-right">
<li>
<a href="https://github.com/LSun/truncash">
<span class="fa fa-github"></span>
</a>
</li>
</ul>
</div><!--/.nav-collapse -->
</div><!--/.container -->
</div><!--/.navbar -->
<!-- Add a small amount of space between sections. -->
<style type="text/css">
div.section {
padding-top: 12px;
}
</style>
<div class="fluid-row" id="header">
<h1 class="title toc-ignore">Normal Means with Heteroskedastic and Correlated Noise</h1>
<h4 class="author"><em>Lei Sun</em></h4>
<h4 class="date"><em>2017-04-12</em></h4>
</div>
<p><strong>Last updated:</strong> 2018-05-12</p>
<strong>workflowr checks:</strong> <small>(Click a bullet for more information)</small>
<ul>
<li>
<details>
<p><summary> <strong style="color:blue;">✔</strong> <strong>R Markdown file:</strong> up-to-date </summary></p>
<p>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.</p>
</details>
</li>
<li>
<details>
<p><summary> <strong style="color:blue;">✔</strong> <strong>Repository version:</strong> <a href="https://github.com/LSun/truncash/tree/5b37742d7cdd7ba50dd641579d2705df67b8b441" target="_blank">5b37742</a> </summary></p>
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. <br><br> 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 <code>wflow_publish</code> or <code>wflow_git_commit</code>). 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:
<pre><code>
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/
</code></pre>
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.
</details>
</li>
</ul>
<details>
<summary> <small><strong>Expand here to see past versions:</strong></small> </summary>
<ul>
<table style="border-collapse:separate; border-spacing:5px;">
<thead>
<tr>
<th style="text-align:left;">
File
</th>
<th style="text-align:left;">
Version
</th>
<th style="text-align:left;">
Author
</th>
<th style="text-align:left;">
Date
</th>
<th style="text-align:left;">
Message
</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:left;">
rmd
</td>
<td style="text-align:left;">
<a href="https://github.com/LSun/truncash/blob/cc0ab8379469bc3726f1508cd81e4ecd6fef1b1a/analysis/ash_gd.rmd" target="_blank">cc0ab83</a>
</td>
<td style="text-align:left;">
Lei Sun
</td>
<td style="text-align:left;">
2018-05-11
</td>
<td style="text-align:left;">
update
</td>
</tr>
<tr>
<td style="text-align:left;">
html
</td>
<td style="text-align:left;">
<a href="https://cdn.rawgit.com/LSun/truncash/0f36d998db26444c5dd01502ea1af7fbd1129b22/docs/ash_gd.html" target="_blank">0f36d99</a>
</td>
<td style="text-align:left;">
LSun
</td>
<td style="text-align:left;">
2017-12-21
</td>
<td style="text-align:left;">
Build site.
</td>
</tr>
<tr>
<td style="text-align:left;">
html
</td>
<td style="text-align:left;">
<a href="https://cdn.rawgit.com/LSun/truncash/853a484bfacf347e109f6c8fb3ffaab5f4d6cc02/docs/ash_gd.html" target="_blank">853a484</a>
</td>
<td style="text-align:left;">
LSun
</td>
<td style="text-align:left;">
2017-11-07
</td>
<td style="text-align:left;">
Build site.
</td>
</tr>
<tr>
<td style="text-align:left;">
html
</td>
<td style="text-align:left;">
<a href="https://cdn.rawgit.com/LSun/truncash/043bf895b4a909cf842ddfa19a7187908f608fa8/docs/ash_gd.html" target="_blank">043bf89</a>
</td>
<td style="text-align:left;">
LSun
</td>
<td style="text-align:left;">
2017-11-05
</td>
<td style="text-align:left;">
transfer
</td>
</tr>
<tr>
<td style="text-align:left;">
rmd
</td>
<td style="text-align:left;">
<a href="https://github.com/LSun/truncash/blob/3c00d157ca8476c29ff324b50e05f1ff3809025f/analysis/ash_gd.rmd" target="_blank">3c00d15</a>
</td>
<td style="text-align:left;">
LSun
</td>
<td style="text-align:left;">
2017-04-21
</td>
<td style="text-align:left;">
de bogged
</td>
</tr>
<tr>
<td style="text-align:left;">
html
</td>
<td style="text-align:left;">
<a href="https://cdn.rawgit.com/LSun/truncash/3c00d157ca8476c29ff324b50e05f1ff3809025f/docs/ash_gd.html" target="_blank">3c00d15</a>
</td>
<td style="text-align:left;">
LSun
</td>
<td style="text-align:left;">
2017-04-21
</td>
<td style="text-align:left;">
de bogged
</td>
</tr>
<tr>
<td style="text-align:left;">
rmd
</td>
<td style="text-align:left;">
<a href="https://github.com/LSun/truncash/blob/239f1952ba3f76b6894649ae7cb530b2f63be4e4/analysis/ash_gd.rmd" target="_blank">239f195</a>
</td>
<td style="text-align:left;">
LSun
</td>
<td style="text-align:left;">
2017-04-21
</td>
<td style="text-align:left;">
GD-ASH
</td>
</tr>
<tr>
<td style="text-align:left;">
html
</td>
<td style="text-align:left;">
<a href="https://cdn.rawgit.com/LSun/truncash/239f1952ba3f76b6894649ae7cb530b2f63be4e4/docs/ash_gd.html" target="_blank">239f195</a>
</td>
<td style="text-align:left;">
LSun
</td>
<td style="text-align:left;">
2017-04-21
</td>
<td style="text-align:left;">
GD-ASH
</td>
</tr>
</tbody>
</table>
</ul>
</details>
<hr />
<div id="problem-setting" class="section level2">
<h2>Problem setting</h2>
<p>We are considering the normal means problem with heteroskedastic noise. Furthermore, The noise are not independent, with unknown pairwise correlation <span class="math inline">\(\rho_{ij}\)</span>.</p>
<p><span class="math display">\[
\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}
\]</span></p>
<p>Under <a href="https://academic.oup.com/biostatistics/article/18/2/275/2557030/False-discovery-rates-a-new-deal"><code>ash</code> framework</a>, a prior is put on exchangeable <span class="math inline">\(\beta_j\)</span>:</p>
<p><span class="math display">\[
g(\beta_j) = \sum_k \pi_k g_k(\beta_j)\ .
\]</span></p>
<p>According to <a href="gaussian_derivatives.html">previous exploration</a>, correlated standard normal <span class="math inline">\(z\)</span> scores can be <em>seen as iid from a density composed of Gaussian derivatives</em>.</p>
<p><span class="math display">\[
f_0(z_j) = \sum_{l=0}^\infty w_l \varphi^{(l)}(z_j) \ .
\]</span> By change of variables, the likelihood</p>
<p><span class="math display">\[
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) \ .
\]</span></p>
<p>Thus for each pair of observation <span class="math inline">\((\hat\beta_j, \hat s_j)\)</span>, the likelihood becomes</p>
<p><span class="math display">\[
\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}
\]</span> Hence, it’s all boiled down to calculate <span class="math inline">\(f_{jkl}\)</span>, which is the convolution of a Gaussian derivative <span class="math inline">\(\varphi^{(l)}\)</span> and a prior component <span class="math inline">\(g_k\)</span>, with some change of variables manipulation. Note that in usual <code>ASH</code> settings, <span class="math inline">\(g_k\)</span> is either uniform or normal, both of which can be handled without conceptual difficulty.</p>
<p>If <span class="math inline">\(g_k\)</span> is a uniform, the convolution of a Gaussian derivative and a uniform is just another Gaussian derivative in a lower order, such as</p>
<p><span class="math display">\[
\varphi^{(l)} * \text{Unif} \leadsto \varphi^{(l-1)} \ .
\]</span> On the other hand, if <span class="math inline">\(g_k\)</span> is a normal, we can use <a href="https://en.wikipedia.org/wiki/Convolution#Differentiation">a fact about convolution</a> that</p>
<p><span class="math display">\[
\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}
\]</span> Because the convolution of <span class="math inline">\(\varphi\)</span>, a Gaussian, and another Gaussian is still a Gaussian, the convolution of a Gaussian derivative and a normal gives another Gaussian derivative.</p>
<p>With <span class="math inline">\(f_{jkl}\)</span> computed, the goal is then to maximize the joint likeliood of the observation <span class="math inline">\(\left\{\left(\hat\beta_1, \hat s_1\right), \cdots,\left(\hat\beta_n, \hat s_n\right) \right\}\)</span> which is <span class="math display">\[
\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) \ ,
\]</span> subject to reasonable, well designed constraints on <span class="math inline">\(\pi_k\)</span> and <a href="gaussian_derivatives_5.html">especially <span class="math inline">\(w_l\)</span></a>.</p>
</div>
<div id="characteristic-function-fourier-transform" class="section level2">
<h2>Characteristic function / Fourier transform</h2>
<p>The exact form of <span class="math inline">\(f_{jkl}\)</span> can be derived analytically. Below is a method using characteristic functions or Fourier transforms.</p>
<p>Let <span class="math inline">\(\zeta_X(t) = \zeta_F(t) =\zeta_f(t) := E[e^{itX}]\)</span> be the characteristic function of the random variable <span class="math inline">\(X\sim dF = f\)</span>. Either the random variable, its distribution, or its density function can be put in the subscript, depending on the circumstances.</p>
<p>On the other hand, the characteristic function of a random variable is also closely related to the Fourier transform of its density. In particular,</p>
<p><span class="math display">\[
\zeta_f(t) = E[e^{itX}] = \int e^{itx}dF(x) = \int e^{itx}f(x)dx := \mathscr{F}_f(t) \ .
\]</span> Note that in usual definition, the Fourier transform is given by <span class="math display">\[
\hat f(\xi) := \int f(x)e^{-2\pi ix\xi}dx \ ,
\]</span> with a normalizing factor <span class="math inline">\(2\pi\)</span> and a negative sign. However, even defined in different ways, <span class="math inline">\(\mathscr{F}_f(t)\)</span> carries over many nice properties of <span class="math inline">\(\hat f(\xi)\)</span>, for example, as we’ll show,</p>
<p><span class="math display">\[
\mathscr{F}_{f^{(m)}}(t) = (-it)^m\mathscr{F}_f(t) \ ,
\]</span></p>
<p>where <span class="math inline">\(f^{(m)}\)</span> is the <span class="math inline">\(m^\text{th}\)</span> derivative of <span class="math inline">\(f\)</span>.</p>
<p>Also, <span class="math display">\[
\mathscr{F}_{f*g}(t) = \mathscr{F}_f(t)\mathscr{F}_g(t) \ ,
\]</span> where <span class="math inline">\(*\)</span> stands for convolution.</p>
<p>With these properties, the Fourier transform tool could be very useful when dealing with Gaussian derivatives and their convolution with other (density) functions.</p>
<p>Under this definition, the inversion formula for the characteristic functions, also known as the inverse Fourier transform, is</p>
<p><span class="math display">\[
f(x) = \frac1{2\pi} \int e^{-itx}\zeta_f(t)dt = \frac{1}{2\pi}\int e^{-itx}\mathscr{F}_f(t)dt \ .
\]</span></p>
<p>Then the characteristic function of <span class="math inline">\(\hat\beta_j | \hat s_j\)</span></p>
<p><span class="math display">\[
\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) \ .
\]</span></p>
<p>Let’s take care of <span class="math inline">\(\zeta_{\beta}(t)\)</span> and <span class="math inline">\(\zeta_{f_0}\left(\hat s_jt\right)\)</span> one by one. Note that</p>
<p><span class="math display">\[
\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}
\]</span></p>
<p>Meanwhile, using the fact that <em>Gaussian derivatives should be absolutely integrable</em>,</p>
<p><span class="math display">\[
\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}
\]</span> Each Fourier transform of gaussian derivatives</p>
<p><span class="math display">\[
\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}
\]</span> Thus,</p>
<p><span class="math display">\[
\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}
\]</span> It is essentially the equivalent expression of <span class="math inline">\(f_{jkl}\)</span> as</p>
<p><span class="math display">\[
\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}
\]</span></p>
<div id="uniform-mixture-prior" class="section level3">
<h3>Uniform mixture prior</h3>
<p>If the prior of <span class="math inline">\(\beta\)</span> is a mixture of uniforms, <span class="math inline">\(\beta \sim g = \sum_k\pi_kg_k\)</span> where <span class="math inline">\(g_k = \text{Unif}\left[a_k, b_k\right]\)</span>,</p>
<p><span class="math display">\[
\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}
\]</span></p>
<p>It looks pretty complicated, but if we take advantage of usual Fourier transform tricks things get clearer.</p>
<p>Let <span class="math inline">\(\varphi_{\mu, \sigma^2}\)</span> denote the density function of <span class="math inline">\(N(\mu, \sigma^2)\)</span>,</p>
<p><span class="math display">\[
\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}
\]</span></p>
<p>Take the first part of <span class="math inline">\(f_{jkl}\)</span>, <span class="math display">\[
\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}
\]</span> Therefore,</p>
<p><span class="math display">\[
\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}
\]</span></p>
</div>
<div id="normal-mixture-prior" class="section level3">
<h3>Normal mixture prior</h3>
<p>Likewise, if the prior of <span class="math inline">\(\beta\)</span> is a mixture of normals, <span class="math inline">\(\beta \sim g = \sum_k\pi_kg_k\)</span> where <span class="math inline">\(g_k = N(\mu_k, \sigma_k^2)\)</span>,</p>
<p><span class="math display">\[
\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}
\]</span> In common applications, <span class="math inline">\(\mu_k\equiv\mu\equiv0\)</span>.</p>
</div>
</div>
<div id="optimization-problem" class="section level2">
<h2>Optimization problem</h2>
<p>The problem is now boiled down to maximize the joint likelihood</p>
<p><span class="math display">\[
\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}
\]</span></p>
<p>subject to appropriate constraints on <span class="math inline">\(\pi_k\)</span> and <a href="gaussian_derivatives_5.html">especially <span class="math inline">\(w_l\)</span></a>, where the specific form of <span class="math inline">\(f_{jkl}\)</span> depends on the mixture component of the prior <span class="math inline">\(g\)</span> of <span class="math inline">\(\beta_j\)</span>. Here we consider two cases.</p>
<div id="uniform-mixture-prior-1" class="section level3">
<h3>Uniform mixture prior</h3>
<p><span class="math display">\[
\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}
\]</span></p>
</div>
<div id="normal-mixture-prior-1" class="section level3">
<h3>Normal mixture prior</h3>
<p><span class="math display">\[
\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}
\]</span></p>
</div>
<div id="biconvex-optimization" class="section level3">
<h3>Biconvex optimization</h3>
<p>Our goal is then to estimate <span class="math inline">\(\hat\pi\)</span> by solving the following constrained biconvex optimization problem</p>
<p><span class="math display">\[
\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}
\]</span></p>
</div>
</div>
<hr>
<p>
</p>
<hr>
<!-- To enable disqus, uncomment the section below and provide your disqus_shortname -->
<!-- disqus
<div id="disqus_thread"></div>
<script type="text/javascript">
/* * * CONFIGURATION VARIABLES: EDIT BEFORE PASTING INTO YOUR WEBPAGE * * */
var disqus_shortname = 'rmarkdown'; // required: replace example with your forum shortname
/* * * DON'T EDIT BELOW THIS LINE * * */
(function() {
var dsq = document.createElement('script'); dsq.type = 'text/javascript'; dsq.async = true;
dsq.src = '//' + disqus_shortname + '.disqus.com/embed.js';
(document.getElementsByTagName('head')[0] || document.getElementsByTagName('body')[0]).appendChild(dsq);
})();
</script>
<noscript>Please enable JavaScript to view the <a href="http://disqus.com/?ref_noscript">comments powered by Disqus.</a></noscript>
<a href="http://disqus.com" class="dsq-brlink">comments powered by <span class="logo-disqus">Disqus</span></a>
-->
<!-- Adjust MathJax settings so that all math formulae are shown using
TeX fonts only; see
http://docs.mathjax.org/en/latest/configuration.html. This will make
the presentation more consistent at the cost of the webpage sometimes
taking slightly longer to load. Note that this only works because the
footer is added to webpages before the MathJax javascript. -->
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
"HTML-CSS": { availableFonts: ["TeX"] }
});
</script>
<hr>
<p>
This reproducible <a href="http://rmarkdown.rstudio.com">R Markdown</a>
analysis was created with
<a href="https://github.com/jdblischak/workflowr">workflowr</a> 1.0.1
</p>
<hr>
</div>
</div>
</div>
<script>
// add bootstrap table styles to pandoc tables
function bootstrapStylePandocTables() {
$('tr.header').parent('thead').parent('table').addClass('table table-condensed');
}
$(document).ready(function () {
bootstrapStylePandocTables();
});
</script>
<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
(function () {
var script = document.createElement("script");
script.type = "text/javascript";
script.src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
document.getElementsByTagName("head")[0].appendChild(script);
})();
</script>
</body>
</html>