写在前面
在中文社区,对于最优设计这块的讨论几近于无。写这篇博客也算是抛砖引玉。
论文:《Robust Experiment Design via Stochastic Approximation 》(1984)
作者:LUC PRONZATO AND ERIC WALTER
介绍了一种统计实验设计方案,利用Fisher矩阵的行列式的期望评判参数估计精度。然而,Fisher矩阵的计算需要知道未知参数的值,这实际上是不能获得的,因此需要一个标称值(nominal value of the parameter.)
ED-Optimal design
D-optimal design
在介绍ED-optimal之前,需要对D-optimal进行介绍。
最大似然估计写为:
θml=argmaxθ[p(y∣θ)]\theta_{ml}=argmax_\theta[p(y|\theta)]θml=argmaxθ[p(y∣θ)]
最大似然估计是渐进正态的
N(θ,M−1(θ))\mathcal{N}(\theta, M^{-1}(\theta))N(θ,M−1(θ))
其中,M(θ)M(\theta)M(θ)是Fisher信息矩阵,定义为
M(θ)=Ey∣θ∂lnp(y∣θ)∂θ∂lnp(y∣θ)∂θTdyM(\theta)=E_{y|\theta}{\frac{\partial lnp(y|\theta)}{\partial \theta}}{\frac{\partial lnp(y|\theta)}{\partial \theta^T}}dyM(θ)=Ey∣θ∂θ∂lnp(y∣θ)∂θT∂lnp(y∣θ)dy
也等于M(θ)=Ey∣θ[∂2lnp(y∣θ)∂θ∂θT]M(\theta)=E_{y|\theta}[\frac{\partial ^2 lnp(y|\theta)}{\partial \theta \partial \theta^T}]M(θ)=Ey∣θ[∂θ∂θT∂2lnp(y∣θ)]Fisher矩阵针对yyy求取了条件分布的期望,因此实际上是关于θ\thetaθ的函数。
在多元线性回归(高斯噪声)的情形下,Fisher信息矩阵的表达式是
M(θ)=XTΣ−1XM(\theta)=X^T\Sigma^{-1}XM(θ)=XTΣ−1X其中,X为数据矩阵(设计矩阵),或者说导数矩阵(对参数向量求导)。
可以发现,这个表达式与θ\thetaθ没有关系了。这个表达式的推导,我会另外写文章详细探讨。
然而,在非线性的情况下,上面的表达式仍然与θ\thetaθ有关系。其原因在于X是导数矩阵,线性情形下,任何点的导数都一致。然而,对于非线性问题,我们实际上是把模型先线性化,得到这个导数矩阵。而这个导数矩阵是在哪个点求得的呢?
实际上,应该从参数的真实值处展开,然而真值是不知道的。D-optimal的做法是,寻找一个标称值θ0\theta_0θ0,在该值处线性化正向模型,再利用线性回归的Fisher矩阵表达式。这个θ0\theta_0θ0则只能依靠先验知识或专家决策。
D-optimal就表达为:
e=argmaxe(det(M(θ0)))e=argmax_e(det(M(\theta_0)))e=argmaxe(det(M(θ0)))
eee表示实验配置。说白了就是找到一组设计变量,使得估计的Fisher矩阵的行列式最大。
ED-optimal
ED-optimal所作出的改进是考虑参数θ\thetaθ的不确定性,实际上也就是先验分布。我们在θ\thetaθ所有可能的值上求取期望,从而获得相比单取一个θ0\theta_0θ0更加鲁棒的结果。ED-optimal的表达式是:
e∗=argmaxe∫θdet(M(θ))p(θ)dθe^*=argmax_e\int_\theta det(M(\theta))p(\theta)d\thetae∗=argmaxe∫θdet(M(θ))p(θ)dθ其中,p(θ)p(\theta)p(θ)我个人理解是先验分布。
另一方面,从决策理论的角度看,这实际上是某种贝叶斯风险。风险函数是负的Fisher信息的行列式。
值得注意的是,上式要区别于下式:
e∗=argmaxedet(M(∫θθp(θ)dθ))e^*=argmax_e\det(M(\int_\theta\theta p(\theta)d\theta))e∗=argmaxedet(M(∫θθp(θ)dθ))即先求期望再计算Fisher信息。
一个例子
令y=e−θx(x≥0)+ϵy=e^{-\theta x} (x \ge0) +\epsilony=e−θx(x≥0)+ϵ
噪声服从ϵ∼N(0,σ2)\epsilon \sim N(0, \sigma^{2})ϵ∼N(0,σ2)
先求出Fisher信息的表达式:
∂2y∂θ2=x2σ2ye−θx−2x2σ2e−2θxM(θ,x)=−Ey∣θ[∂2y∂θ2]=−x2σ2e−θxEy∣θ[y]+2x2σ2e−2θx=−x2σ2e−2θx+2x2σ2e−2θx=x2σ2e−2θx
\begin{aligned}
\frac{\partial ^2y}{\partial \theta ^2}&=\frac{x^2}{\sigma^2}ye^{-\theta x}-2\frac{x^2}{\sigma^2}e^{-2 \theta x}\\
M(\theta,x)&=-E_{y|\theta}[\frac{\partial ^2y}{\partial \theta ^2}]\\
&=-\frac{x^2}{\sigma^2}e^{-\theta x}E_{y|\theta}[y]+2\frac{x^2}{\sigma^2}e^{-2 \theta x}\\
&=-\frac{x^2}{\sigma^2}e^{-2\theta x}+2\frac{x^2}{\sigma^2}e^{-2 \theta x}\\
&=\frac{x^2}{\sigma^2}e^{-2\theta x}
\end{aligned}
∂θ2∂2yM(θ,x)=σ2x2ye−θx−2σ2x2e−2θx=−Ey∣θ[∂θ2∂2y]=−σ2x2e−θxEy∣θ[y]+2σ2x2e−2θx=−σ2x2e−2θx+2σ2x2e−2θx=σ2x2e−2θx注意到和参数θ\thetaθ和设计变量xxx都有关系。
ED-optimal的表达式就是
∫θx2σ2e−2θxp(θ)dθ\int_\theta\frac{x^2}{\sigma^2}e^{-2\theta x}p(\theta)d\theta∫θσ2x2e−2θxp(θ)dθ设置θ\thetaθ为[a,b][a,b][a,b]上的均匀分布,则上面的积分为:
−x2σ2(b−a)(e−2xb−e−2xa)-\frac{x}{2\sigma^2(b-a)}(e^{-2xb}-e^{-2xa})−2σ2(b−a)x(e−2xb−e−2xa)
最优化的过程省略。