Unexpectedly, the theoretically best reject-domain of T-test is bounded.

五月 11th, 2008 - 2 Responses

 f_{t}\left(x;\mu,df\right)\equiv C\left(df\right)\left(1+\frac{\left(x-\mu\right)^{2}}{df}\right)^{-\frac{df+1}{2}}

 \begin{eqnarray*}\lambda\left(x;\mu_{0},\mu_{1},df\right) & \equiv & \frac{f_{t}\left(x;\mu_{1},df\right)}{f_{t}\left(x;\mu_{0},df\right)}\\ & = & \left(\frac{v+\left(x-\mu_{1}\right)^{2}}{v+\left(x-\mu_{0}\right)^{2}}\right)^{-\frac{df+1}{2}}\\ & {\longrightarrow\atop x\rightarrow\infty} & 1\end{eqnarray*}

For NHST  H_{0}:T\sim t_{df} vs  H_{1}:T-1\sim t_{df}, theoretically,  p\left(t\right)=\int_{\left\{ x:\lambda\left(x\right)\ge\lambda\left(t\right)\right\} }f_{t}\left(x,\mu_{0},df\right)dx is s.t.  \lim_{t\rightarrow\infty}p\left(t\right)=\frac{1}{2} , rather than zero. Nevertheless, pratically a large t, rejecting both  H_{0} and  H_{1}, should not be counted as any evidence to retain or reject  H_{0}.

To verify the shape of  \lambda\left(x\right)

x<-.01*(-1000:1000);
plot(x,y=dt(x-1,df=5,ncp=0)/dt(x,df=5,ncp=0));
###Comparison the noncetrality case--
## plot(x,y=dt(x,df=5,ncp=1)/dt(x,df=5,ncp=0));

Confidence Domain and Not-reject Domain

五月 10th, 2008 - No Responses

Either Confidence Interval (CI) or Null Hypothesis Significance Test (NHST) has the same business, to advise whether some sample  X\equiv\left(X_{1},X_{2},\dots,X_{n}\right) is or is not disliked by some hypothesized parameter  \vartheta.

NHST.com manages a database. For each Miss  \vartheta, NHST spies out all she dislikes. Mr X logs in NHST.com and inputs a girl name and his credit card number, to bet his luck whispering– Does she dislikes me?

CI.com manages a database too. For each Mr X, CI only needs his credit card with his name X on it, then serves him a full list of available girl names.

NHST.com has been historically monopolizing the market. Nevertheless, somebody prefer the service of CI.com and find that the two may share database in most cases.

Not-reject Domain of  \vartheta is defined as  A\left(\vartheta\right)=\left\{ x:\vartheta\; doesn't\;dislike\;x\right\} .

Confidence Domain of x is defined as  S\left(x\right)\equiv\left\{ \vartheta:\vartheta\; doesn't\;dislike\;x\right\} .

 \begin{eqnarray*}\theta &amp; \in &amp; S\left(X\right)\\ &amp; \Updownarrow\\\theta\;doesn't\; &amp; dislikes &amp; X\\ &amp; \Updownarrow\\X &amp; \in &amp;A\left(\theta\right)\end{eqnarray*}

So,  Pr_{\vartheta}\left(\vartheta\in S\left(X\right)\right)\ge1-\alpha,\forall\vartheta\Longleftrightarrow Pr_{\vartheta}\left(X\notin A\left(\vartheta\right)\right)\le\alpha,\forall\vartheta

Automatize LISREL jobs

四月 28th, 2008 - One Response

LISREL routine can run in DOS or in command line mode of windows (windows-key + R -> CMD) . The command line is just like –

D:\My Documents>“C:\Program Files\lisrel87\lisrel87.exe” “C:\Program Files\lisrel87\LS8EX\EX61.LS8″ D:\myOutput.out

1. You only need edit and input the bold part.
2. Quotation marks are used wherever the paths or filenames include blanks.
3. The 2nd argument is the output file. You can still specify more output options in your .ls8 file.
4. A .bat file can automatize batches of such lisrel jobs.

Developing normal pdf from symmetry & independence

四月 26th, 2008 - No Responses

When I was in the 3rd grade of my middle school, I enjoyed my town bookstore as a standing library. There a series of six math-story books by Zhang Yuan-Nan impressed me a lot. I cited a case from one in my PPT when I taught the normal distribution — the normal pdf can be derived from simple symmetry & independence conditions.

Today I can even google out an illegal pdf of its new edition to verify the case (2005, pp. 89). Actually I have bought the new edition series (now 3 books) and lent them to students. Those conditions are as instinctive as–

1. For white noise errors on 2-D, the independence means pdf at  (x,y) is the product of 1-D pdf, that is,  \phi\left(x\right)\phi\left(y\right) .

2. The symmetry means pdf at  (x,y) is just a function of  x^{2}+y^{2}, nothing to do with direction. That is,  \phi\left(x\right)\phi\left(y\right)=f\left(x^{2}+y^{2}\right).

So,   f\left(x^{2}\right)f\left(y^{2}\right)=f\left(x^{2}+0\right)f\left(0+y^{2}\right)=\phi^{2}\left(0\right)f\left(x^{2}+y^{2}\right).

For middle school students, the book stated a gap here to arrive at the final result   f\left(x^{2}\right)=ke^{bx^{2}}, which is    \phi\left(x\right)=\frac{1}{\phi\left(0\right)}f\left(x^{2}+0\right)=\frac{k}{\phi\left(0\right)}e^{bx^{2}}.

I think non-math graduate students with interests can close the gap by themselves with following small hints.

Denote  \alpha=x^{2},\beta=y^{2}. We have  \log f\left(\alpha\right)+\log f\left(\beta\right)=\log\phi^{2}\left(0\right)+\log f\left(\alpha+\beta\right), or  \left[\log f\left(\alpha\right)-\log\phi^{2}\left(0\right)\right]+\left[\log f\left(\beta\right)-\log\phi^{2}\left(0\right)\right]=\left[\log f\left(\alpha+\beta\right)-\log\phi^{2}\left(0\right)\right]. Denote  g\left(\alpha\right)=\log f\left(\alpha\right)-\log\phi^{2}\left(0\right). That is,   g\left(\alpha\right)+g\left(\beta\right)=g\left(\alpha+\beta\right).

Now to prove  g\left(\frac{m}{n}\right)=\frac{m}{n}g\left(1\right),\forall m,n\in\mathbb{N}. With continuousness, it gets  g\left(\alpha\right)=\alpha g\left(1\right),\forall\alpha\in\mathbb{R}.

愉快地发现SciTE和LyX在WinXP下都支持中文

二月 21st, 2008 - No Responses

愉快地发现我在WinXP上最常用的编辑工具SciTE(当前版本1.75) 和 LyX(当前版本1.5.3) 都支持unicode(也就是说,支持中文)。之前不了解,只因为缺省设置不支持中文。需要手工操作修改设置。

SciTE的设置是Options->Open Global Options File,编辑SciTEGlobal.properties,找到如下段落

# Unicode
#code.page=65001
code.page=0
#character.set=204
# Required for Unicode to work on GTK+:
#LC_CTYPE=en_US.UTF-8
#output.code.page=65001

修改为

# Unicode
code.page=65001
#code.page=0
character.set=204
# Required for Unicode to work on GTK+:
LC_CTYPE=en_US.UTF-8
output.code.page=65001

保存。然后关闭再打开SciTE,就会发现不再出现中文被切一半的现象。如果编辑的文档格式不是utf-8而是ucs-2 ,还可以在File->Encoding 里临时选。

此外,我还推荐把SciTEGlobal.properties文件中的line.margin.visible=1 和 wrap=1 两处的注释#号去掉,效果是缺省显示行号,并使超长的行折行显示。SciTE的优点太多了–开源免费;轻巧,启动快;支持Ctrl+鼠标中轮滚动无级缩放;支持Ctrl+回车 出现单词补齐选项;支持Alt键方形选段;…

LyX(版本>=1.5.1)在winXP已经可以在.lyx文件正文和公式框中录入中文。麻烦的是输出中文的pdf。要实现这点,首先需要在LyX菜单Document->Settings->Language->选Chinese(simplified),并勾选 Use language’s default encoding设为English + GBK encoding组合(发现Language选项只影响F7拼写检查,而Encoding只影响pdf/dvi生成;还可以点Save as Document Defaults设为今后的缺省)。我试了若干方案,最后才试成功用南开网站上的MiKTeX+中文插件(Patches4miktex.exe),实现了输出正文的中文为pdf,美中不足是公式框中还不能实现中文公式输出pdf,甚至在Document中的中文设置选定后,公式框中的中文都不能在LyX环境中正确显示,期待后续的LyX版本处理好这个问题。

新版本LyX已经引入了文档版本控制,相当于word中的revision功能,有待深入试用。目前LyX仍不支持Ctrl+鼠标中轮滚动无级缩放,如果公式显得太小,需要在菜单设置中修改显示缩放比例:Tools->Preferences->Look and feel->Screen fonts->Zoom %。这可能是比较容易在后续版本中实现的功能。

相关网址:SciTE主页 http://www.scintilla.org/SciTE.html LyX主页 http://lyx.org/ 南开MiKTeX中文插件 http://miktex.math.nankai.edu.cn/

My first wp plugin work: LaTeX_Math_cgi 1.0

二月 8th, 2008 - No Responses


Install: Download LaTeX_Math_cgi.zip and unzip it to your wordpress’ /wp-content/plugins/ directory.

Activate it in Plugins menu. View its options in Options > LaTeX tab –

It is actually used as a mimeTeX plugin rather than just a  L^{A}T_{E}X plugin. My contribution is technically trivial. But, you must need it if you change from a light theme to a cool black one while find the default mimeTeX images are black in front and transparent in background. This plugin provides mimeTeX’s \reverse and \opaque options to tune your math forms to your wp theme without editing them one by one.

The 2nd function is the option for your cgi URL. The default is http://www.forkosh.dreamhost.com/mimetex.cgi You can DIY one following Forkosh’s instructions and then set it like /cgi-bin/mimetex.cgi It can help your visitors not to cost Forkosh’s unselfish bandwidth.

Understanding QQ plots

一月 30th, 2008 - No Responses

## Try distributions like rchisq, rt, runif, rf to view its heavy, or light, left, or right tail.

n <- 60;
ry <- rnorm(n);
qqnorm(ry);
max(ry)
min(ry)
##view and guess what are x(s) and y(s)
I <- rep(1,n);
qr <- ((ry%*%t(I) > I %*% t(ry))+.5*(ry %*% t(I) == I%*%t(ry)))%*%I *(1/n);
##qr are the sample quantiles
points(qr,ry,col="blue");
##to view the fact, try the following
points(qr,qr,col="green");
rx <- qnorm(qr);
points(rx,ry,col="red",pch="O");
##Red O(s) circle black o(s) exactly.

03DEC2007 R-workshop sponsored by dept of psy, ZSU(=SYSU, Guang-Zhou)

十二月 3rd, 2007 - No Responses

Here is the updated PPT for the talk in the afternoon–which includes the zipped example codes and set-up steps for the workshop in the evening within the 3rd page. The listed anonymous on-line test (result statistics) on p-value interpretation was cited indirectly From Gigerenzer, Krauss, & Vitouch (2004).

There is an advert on http://www.psy.sysu.edu.cn/detail_news.asp?id=258 and a formal CV of the speaker is available on http://lixiaoxu.googlepageS.com

Classic Neyman-Pearson approach demo

十一月 30th, 2007 - 3 Responses

It notes here that N-P approach does not utilize the information in the accurate p value. Actually, at the time N-P approach was firstly devised, the accurate p value was not available usually. Now almost all statistic softwares provide accurate p values and the N-P approach becomes obsolete. Wilkinson & APA TFSI (1999) recommended to report the accurate p value rather than just significance/insignificance, unless p is smaller than any meaningful precision.

p_prior<-0.22;
##try p_prior<-0.82 to demo significant repetition once more
p_alpha<-0.05;
p_power<-0.8;
##-----------------
op <- par(mfrow=c(1,2));
h_1_s<-p_power*p_prior/(p_alpha*(1-p_prior)+p_power*p_prior);
h_0_i<-(1-p_alpha)*(1-p_prior)/((1-p_alpha)*(1-p_prior)+(1-p_power)*p_prior);
##
D<-matrix(c(1-p_alpha,p_alpha,1-p_power,p_power),nrow=2);
rownames(D)<-c("Insig","Sig");
colnames(D)<-c(paste(1-p_prior,"H_0"),paste(p_prior,"H_1"));
barplot(D
,width=c(1-p_prior,p_prior)
,legend = rownames(D)
,space=0
,col=c(gray(0.6),gray(0.9)),asp=1
);
title(main=paste("Sig: H_1 from",p_prior,"to"
,round(h_1_s,digits=2)
,"\n","Insig: H_0 from",1-p_prior,"to"
,round(h_0_i,digits=2)
)
,sub=paste(round(h_1_s,digits=2),"=(",p_power,"*",p_prior,")/("
,p_alpha,"*",1-p_prior,"+",p_power,"*",p_prior,")"
,"\n",round(h_0_i,digits=2)
,"=(",1-p_alpha,"*",1-p_prior,")/("
,1-p_alpha,"*",1-p_prior,"+",1-p_power,"*",p_prior,")"
)
);
##
v_power<-c(0.1,(5:9)/10);
v_prior<-(1:1000)/1000;
n<-length(v_power);
v_c<-sample(colors())[1:n];
plot(v_prior,v_prior,type="l",col="black"
,xlab="a_prior",ylab="post_hoc",asp=1
,xlim=c(0,1),ylim=c(0,1),main="Power=.1,.5,.6,.7,.8, and .9");
lines(v_prior,1-v_prior);
for (i in 1:n){
p_power<-v_power[i];
v_1_s<-p_power*v_prior/(p_alpha*(1-v_prior)+p_power*v_prior);
v_0_i<-(1-p_alpha)*(1-v_prior)/((1-p_alpha)*(1-v_prior)+(1-p_power)*v_prior);
lines(v_prior,v_1_s,col=v_c[i])
lines(v_prior,v_0_i,col=v_c[i])
}
points(p_prior,h_1_s,col="red");
points(p_prior,p_prior,col="red");
points(p_prior,h_0_i,col="blue");
points(p_prior,1-p_prior,col="blue");
par(op);

Wilkinson, L. & APA TFSI (1999). Statistical methods in psychology journals: Guidelines and explanations. American Psychologist, 54, 594-604.

Different corr(s) of different IV scopes with same regression coef

十一月 22nd, 2007 - 2 Responses

 Y=\alpha+\beta X+\varepsilon,\;\varepsilon\sim N\left(0,\sigma^{2}\right)

With  \alpha,\beta, \sigma known in the linear relationship, can the correlation in the scatter plot of Y against X be estimated from the linear formula?

You may recall in Hierarchical Linear Model class, the scopes of the W dramatically impact the regression coefficients of F~W in the following R demo (hlm.jpg). While this time the regression coefficient has been fixed to a known  \beta. So the scopes of X would never impact the regression coefficient. However, it proved that the correlation r could range from zero to unit (or -1) according to the variance of X in the final close form  r=\frac{\beta\mbox{Var}\left(X\right)}{\mbox{Std}\left(Y\right)\mbox{Std}\left(X\right)}=\beta\frac{\mbox{Std}\left(X\right)}{\sqrt{\beta^{2}\mbox{Var}\left(X\right)+\sigma^{2}}}.

Let me quote as the final words from Cohen (1994; p.1001; Where the role of IV is replaced by that of DV within typical contexts like ANOVA) –

… standardized effect size measures, such as d and f, developed in power analysis (Cohen, 1988) are, like correlations, also dependent on population variability of the dependent variable and are properly used only when that fact is kept in mind.

Cohen, J. (1994). The earth is round (p<.05). American Psychologist, 49, 997-1003.

Compare to the following case: different corr(s) of different IV scopes with hierarchical regression coefficients –

### -----
m<-200; ## NO of Inc.
MeanW<-runif(m,max=10000,min=2000);
MeanF<-(MeanW-mean(MeanW))/sd(MeanW) *sqrt(0.8)+sqrt(0.2)*rnorm(m);
Beta<-(MeanW-mean(MeanW))/sd(MeanW)*sqrt(0.8)+sqrt(0.2)*rnorm(m);
Beta<- (-Beta)*2/(max(Beta)-min(Beta)); ## Beta ranges from 1 to 0
J<-c(1:m);
##
## to view Inc. level data
## plot(data.frame(MeanW,MeanF,Beta));
##
n<-6; ## Each Inc. includes n samples.
s<-c(rep(1,n));
MeanW<-c(t(MeanW%*%t(s)));
MeanF<-c(t(MeanF%*%t(s)));
Beta<-c(t(Beta%*%t(s)));
J<-c(t(J%*%t(s)));
###
CW<-rnorm(m*n,sd=100);
CF<-CW*Beta/(sd(MeanW)*sd(MeanF))+rnorm(m*n,sd=0.02);
W<-MeanW+CW;
F<-MeanF+CF;
###
df<-data.frame(W,F,CW,CF,MeanW,MeanF,J);
##
## to valify it in SAS Proc Mixed or SPSS Mixed Model
## write.csv(df,choose.files());
##
plot(W,F,col=sample(colors())[J]); ##
##
library(nlme);
df.int <- groupedData( F ~ CW | J, data=df) ;
## HLM
## F_{i,j}=\beta_{0,j}+\beta_{1,j}*CW + r_{i,j}
## \beta_{0,j}=\gamma_{0,0} + \gamma_{0,1}*MeanW + u_{0,j}
## \beta_{1,j}=\gamma_{1,0} + \gamma_{1,1}*MeanW + u_{1,j}
## Mixed Model
## F_{i,j}=\gamma_{0,0} + \gamma_{0,1}*MeanW + u_{0,j} +(\gamma_{1,0} + \gamma_{1,1}*MeanW + u_{1,j} ) *CW + r_{i,j}
##
lm_1 <- lme(F ~ MeanW+CW+MeanW:CW, data=df.int, random = ~1 + CW | J) ;
summary(lm_1);
intervals(lm_1);
anova(lm_1);
## plot(lm_1); ## view the residual plot
###################