drewlee
比如有一个样本,里面的Xi~N(u,sigma^2)
我们估计u的时候是说P(u1<u<u2)=1-alpha,那么是不是可以通过这个等式得出u是服从某个分布的随机变量的随机变量呢?以Bayesion的观点来看,这个分布应正好是F(u|x)对吗?
通过推导,u的分布是以{xi}为参数的t分布的变换,其中的变换包括平移和伸缩。
其实这里难以解释的是一旦有人通过选取u和sigma以后,模拟得到{xi},然后让我估计u,我说:“u现在是服从t分布的变换后得到的某个分布,我能告诉你它在任何一个区间内的概率值”,那人便说:“难道X={xi},就证明u服从的是你说的那个分布?证明给我看!”,我便哑巴了,因为此时我无法通过以X={xi}作为参数,模拟出u的分布来。所以我无法说服那个人。但是我也无法说服自己不相信自己的观点,因为当我知道某一个未知的量在任何区间的概率时,我就完全有理由说它服从某一个分布。
也许会有争议。Fisher应该不支持我,但是u是不是服从某个分布这个问题是无法证实和证伪的,因为P(X={xi}|X={xi})=0.
小木匠
如果你早出生几百年 贝叶斯学派就会被你取代
yihui
关于区间估计的问题,至少国内的数理统计教科书都把你的这个问题给埋没掉了。其实是一个很关键的问题,需要解释清楚。
P(u1<u<u2)=1-alpha究竟是什么意思?我记得历史上还有两位统计学家产生过争论。不幸我只记得其中一种观点了,就是说u不是随机变量,而是固定的真值。之所以有上面的式子,原因不在于u是随机变量,真正是随机变量的是u1和u2,关于这一点,我想只要明白统计量的来历,就会很容易理解。
为什么会造成“u是随机变量”的印象,原因仅仅在于我们的书写习惯而已。如果写成P(u1<u, u2>u)=1-alpha,我想效果会好一些。
leehman
频率派的观点是指 P(u1(X)<u<u2(X))=1-alpha,既随即区间(u1(X) , u2(X) )包含参数u的概率为1-alpha。
drewlee
[quote]
引用第1楼小木匠于2006-10-16 09:41发表的“”:
如果你早出生几百年 贝叶斯学派就会被你取代[/quote]
谢谢你的支持。
其实Bayes的方法我也很赞同的。这里我自己都觉得奇怪的是怎么能够不用假定先验分布,也就是关于参数u的边缘分布,而直接得出了参数u的条件分布。而且用程序模拟了以后也是对的。
我的程序如下:
<br />
#指定样本的size<br />
N1<-20<br />
#指定方差<br />
var<-1<br />
#存放模拟数<br />
data<-c()<br />
#存放IC包括了u的事件数<br />
k<-0<br />
#指定模拟的次数<br />
N<-1000<br />
for(i in 1:N)<br />
{<br />
#指定均值<br />
u<-runif(1,55,100)<br />
data<-rnorm(N1,u,var)<br />
<br />
#如果IC包括了u,则 k=k+1<br />
if(mean(data)-stdev(data)*qt(1.95/2,N1-1)/sqrt(N1)<=u&mean(data)+stdev(data)*qt(1.95/2,N1-1)/sqrt(N1)>=u)<br />
{k<-k+1}<br />
<br />
}<br />
k<br />
drewlee
[quote]引用第3楼leehman于2006-10-16 14:47发表的“”:
频率派的观点是指 P(u1(X)<u<u2(X))=1-alpha,既随即区间(u1(X) , u2(X) )包含参数u的概率为1-alpha。[/quote]
我认为这就足以得出u的条件分布形式了吧?
-(这一点我现在不这么明朗了)
yihui
(u1(X) , u2(X) )包含参数u的概率为1-alpha,说的是u1 u2在变化,而u是固定值。
drewlee
[quote]引用第2楼谢益辉于2006-10-16 13:36发表的“”:
关于区间估计的问题,至少国内的数理统计教科书都把你的这个问题给埋没掉了。其实是一个很关键的问题,需要解释清楚。
P(u1<u<u2)=1-alpha究竟是什么意思?我记得历史上还有两位统计学家产生过争论。不幸我只记得其中一种观点了,就是说u不是随机变量,而是固定的真值。之所以有上面的式子,原因不在于u是随机变量,真正是随机变量的是u1和u2,关于这一点,我想只要明白统计量的来历,就会很容易理解。
为什么会造成“u是随机变量”的印象,原因仅仅在于我们的书写习惯而已。如果写成P(u1<u, u2>u)=1-alpha,我想效果会好一些。[/quote]
我想听听你的意见。我打算把问题说的更清楚,矛盾更激化一点。(说的不对的地方请指正)
疑点一:我说:"某一个变量(先不说随机变量了)在某一个区间的概率是p(x)",那么如果证明我的这句话是对的还是错的呢?首先我们要确定的是这个变量指的是哪个变量。在这里,特殊的,我们指定这个变量是上帝为了模拟产生size=n的数据,通过电脑设定的那个u。
说到这里再插一点闲话,因为Frequentist和Bayesian跑过来了,Frequentist说:“u哪里是随机的,只有x才能是随机的,只有上帝都不能预料的东西才能是随机。”;Bayesian反驳:“哎,你错了,上帝那么懒惰,他哪里会亲自去设定u啊,为了迷惑我们不至于让所有现象都看上去服从一个规律,他的所有的u也都是模拟出来的随机数,今天我们看到的x是由无数个u里面随便抽出来一个而已。上帝自己也记不住他今天抽出来的是哪个u了。”
书归正传,确定对象了以后,我们开始证明了,我们想到原来的老办法,也就是让记录这个变量的每一个值,然后得到经验分布。可惜我们唯一的救命稻草救不了我们了。因为u不能被记录啊(这里应该是测度论里面的不可测了吧?没学过测度,以后我一定要学)。所以那句话的提出似乎是没有意义的,因为那跟我说上帝一小时前去了一趟伊甸园的概率是1/2是一回事。所以这句话应该没有什么争论才对。即使是Bayes,他在这里提条件概率分布也是没有意义的,u根本就不可观测。
疑点二:事情还有。Frequentist和Bayesian分别得出一套方法估计区间。我们再假想假如他们谁的区间没有包括住u,则谁便会挨一鞭子。习惯性的,他俩都取alpha=0.05,然后我们来料想一下,1000次上帝的玩弄之后,谁挨的鞭子比较多。我的程序显示,假如上帝选择u时用的是随机数,而Bayesian又恰好猜中了上帝的把戏,那么Frequentist和Bayesian挨的鞭子应该一样多,但是假如Bayesian没有正确的猜中u的分布,则很可能Bayesian会比较倒霉的多挨点鞭子。难道Frequentist的方法更好,他更恰如其分的描述了这个只有上帝才知道的量?
drewlee
[quote]引用第6楼谢益辉于2006-10-16 16:32发表的“”:
(u1(X) , u2(X) )包含参数u的概率为1-alpha,说的是u1 u2在变化,而u是固定值。[/quote]
假如把u1(X) , u2(X)看成是随机变量也行,它俩表示估计者的选定的估计函数,但是一旦估计者选定了估计函数u1,和u2的形式以后,又把X给观察到了,即X=x,那么u1(x),u2(x),u不都是常数了吗?P(u1(x)<u<u2(x))的意义就成了判断三个常数之间关系的概率了,就像求P(1<2<3)或者P(3<2<1)一样啊。要么是1要么是0。
yihui
有意思,这个问题有点哲学味道了。我们可以在前人的讨论上继续推进了:)不过今晚我没空,明晚再来说吧
drewlee
是有点哲学了:)
就是那种证实和证伪都不行的命题。
明天听你的高见!
micro@
the underlying question is "how the probability is defined". frequentists and bayesians are using different probability definitions.
cran
很好的问题,在mixed model里面,variance component就是我们要研究的东西。由于sigma^2是非负数,所以他们的MLE分布不是像一般的MLE那样的。
比如我们研究一个正态分布的数据,N(mu,sigma^2)的时候,我们要得到mu 和 sigma 的MLE,然后再研究VAR(mu^hat) 和VAR(sigma^hat)
drewlee
[quote]引用第12楼cran于2006-10-17 13:02发表的“”:
很好的问题,在mixed model里面,variance component就是我们要研究的东西。由于sigma^2是非负数,所以他们的MLE分布不是像一般的MLE那样的。
比如我们研究一个正态分布的数据,N(mu,sigma^2)的时候,我们要得到mu 和 sigma 的MLE,然后再研究VAR(mu^hat) 和VAR(sigma^hat)[/quote]
我只是本科,目前还没接触到mixed model. 不过听你的介绍挺有意思的。
不知道你们研究的问题和model selection 有关吗?我对这个比较感兴趣。里面也有大量的computing. 因为最近在探索cross-validation的新方法。
cran
[quote]引用第13楼drewlee于2006-10-17 17:16发表的“”:
我只是本科,目前还没接触到mixed model. 不过听你的介绍挺有意思的。
.......[/quote]
跟model selection没有关。关于variance component的解很复杂,他的likelihood function也很复杂,需要用积分做出来
比如
f(x;beta,sigma^2) = Int(f(x;beta,u)*g(u;sigma^2))
其中u~Normal(0,sigma^2)
rtist
[quote]引用第14楼cran于2006-10-18 06:46发表的“”:
跟model selection没有关。关于variance component的解很复杂,他的likelihood function也很复杂,需要用积分做出来
.......[/quote]
mixed model的选择更难一些,几乎都需要某种形式的近似。
cran
[quote]引用第15楼rtist于2006-10-19 21:19发表的“”:
mixed model的选择更难一些,几乎都需要某种形式的近似。[/quote]
Approximation是最大的问题,frequentist 方面的我已经接近放弃了,现在用MCMC
hangover
如果你想做模型比较的话,可以看一下Chib's Marginal Likelihood, 从Bayesian的角度去看。也可以用DIC,harmonic mean... cran对这些应该挺熟的。
如果做chib' marginal likelihood需要一个particle filter+ multivariate kernel estimation就完全解决了。做kernel estimation 我建议用matlab KDE的包,就我所知这个现在是最快的,免费的。particle filter就很简单了,随便找篇论文看看。 说到最后,难还是难在模型评估, 太费时间,不搞也罢。 我现在做的是Dynamic Trivariate Mixture model, 因为模型里面的log(time varying variance)是AR(1),所以用的Metropolis-Hastings跑起来很吃力, 收敛奇慢。最后觉得用laplace+kalman filter效果挺好的, 误差在0.01左右,然后结合简单的MCMC算法去做参数。即使这样也跑了将近一个小时,C++ code, 8000数据,3-dimensional obs. non-linear, non normal + 2 lantent processes AR(1) , 7000 MCMC iterations.
根导师研究了一下,他让我跑500-1000 套simulation数据看看算法怎么样,算了一下,500数据也要半个月才能跑完。时间到不是问题,我想的是,程序还没跑到一天random generator 也循环n个来回了,这样做跑出来的数据会不会有问题,问题多大? 有过这方面经验的兄弟们给点建议吧。
hangover
顺便加一句,DIC对mixtur model不准。结果都怪死了,呵呵,现在还没有人解释出来为什么,只是提出这个现象。
cran
[quote]引用第18楼hangover于2006-11-02 02:05发表的“”:
顺便加一句,DIC对mixtur model不准。结果都怪死了,呵呵,现在还没有人解释出来为什么,只是提出这个现象。[/quote]
在做Hierirchical model的时候,有负数的DIC和pD,不知道你有没有发现这个问题。
太好了,遇见同行了。