Every time I think I know what's going on, suddenly there's another layer of complications.
2015年4月29日星期三
2015年4月13日星期一
2015年4月11日星期六
2015年4月10日星期五
Use Newton-Raphson Method to calculate MLE of a Cauchy distribtuion
mlecauchy=function(x,toler=.001){ #x is a vector here
startvalue=median(x)
n=length(x);
thetahatcurr=startvalue;
# Compute first deriviative of log likelihood
firstderivll=2*sum((x-thetahatcurr)/(1+(x-thetahatcurr)^2))
# Continue Newton’s method until the first derivative
# of the likelihood is within toler of 0.001
while(abs(firstderivll)>toler){
# Compute second derivative of log likelihood
secondderivll=2*sum(((x-thetahatcurr)^2-1)/(1+(x-thetahatcurr)^2)^2);
# Newton’s method update of estimate of theta
thetahatnew=thetahatcurr-firstderivll/secondderivll;
thetahatcurr=thetahatnew;
# Compute first derivative of log likelihood
firstderivll=2*sum((x-thetahatcurr)/(1+(x-thetahatcurr)^2))
}
list(thetahat=thetahatcurr);
}
x<-c(-1.94,0.59,-5.98,-0.08,-0.77)
mlecauchy(x,0.0001)
###########result########'
#### $thetahat
#### [1] -0.5343968
#####use the R biuld in funcion ######
optimize(function(theta) -sum(dcauchy(x, location=theta, log=TRUE)), c(-100,100)) #we use negative sign here,
#by default it return a minimum
#optimize(function(theta) sum(dcauchy(x, location=theta, log=TRUE)),maximum=T, c(-100,100))
#########Result###############
### $minimum
### [1] -0.5343902
### $objective
### [1] 11.2959
Almost the same for the two methods.
2015年3月31日星期二
Two sided p value for one sample mean difference-bootstrap mehtod
boottestonemean<-function(x,theta0,b){
#
# x = sample
# theta0 is the null value of the mean
# b is the number of bootstrap resamples
#
# origtest contains the value of the test statistics
# for the original sample
# pvalue is the bootstrap p-value
# teststatall contains the b bootstrap tests
#
n<-length(x)
v<-mean(x)
z<-x-v+theta0
d<-mean(x)-theta0
counter<-0
teststatall<-rep(0,b)
for(i in 1:b){xstar<-sample(z,n,replace=T)
vstar<-mean(xstar)
#if(vstar-theta0 >= d){counter<-counter+1}# one sided test
if(abs(vstar-theta0) >= d){counter<-counter+1} #two sided
teststatall[i]<-vstar}
pvalue<-counter/b #if we do the count *2 it has the same effect as above
list(origtest=v,pvalue=pvalue,teststatall=teststatall)
#list(origtest=v,pvalue=pvalue)
}
x<-c(119.7,104.1,92.8,85.4,108.6,93.4,67.1,88.4,101.0,97.2,95.4,77.2,100.0,114.2,150.3,102.3,105.8,107.5,0.9,94.1)
R<-boottestonemean(x,90,3000)
R$pvalue
hist(R$teststatall)
T<-(mean(x)-90)/(sd(x)/20^0.5) #one sided t test
T
1-pt(T,19)
#
# x = sample
# theta0 is the null value of the mean
# b is the number of bootstrap resamples
#
# origtest contains the value of the test statistics
# for the original sample
# pvalue is the bootstrap p-value
# teststatall contains the b bootstrap tests
#
n<-length(x)
v<-mean(x)
z<-x-v+theta0
d<-mean(x)-theta0
counter<-0
teststatall<-rep(0,b)
for(i in 1:b){xstar<-sample(z,n,replace=T)
vstar<-mean(xstar)
#if(vstar-theta0 >= d){counter<-counter+1}# one sided test
if(abs(vstar-theta0) >= d){counter<-counter+1} #two sided
teststatall[i]<-vstar}
pvalue<-counter/b #if we do the count *2 it has the same effect as above
list(origtest=v,pvalue=pvalue,teststatall=teststatall)
#list(origtest=v,pvalue=pvalue)
}
x<-c(119.7,104.1,92.8,85.4,108.6,93.4,67.1,88.4,101.0,97.2,95.4,77.2,100.0,114.2,150.3,102.3,105.8,107.5,0.9,94.1)
R<-boottestonemean(x,90,3000)
R$pvalue
hist(R$teststatall)
T<-(mean(x)-90)/(sd(x)/20^0.5) #one sided t test
T
1-pt(T,19)
2015年3月30日星期一
one sided or two sided p value for permutation and bootstrap test
The one-sided p-value of the test is calculated as the proportion of sampled
permutations where the difference in means was greater than or equal to T(obs).
The two-sided p-value of the test is calculated as the proportion of sampled
permutations where the absolute difference was greater than or
equal to ABS(T(obs)).
From wiki
From wiki
订阅:
博文 (Atom)