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)




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

Hang-in and Have Smart Friends - The Road to HIV Resistance-by Professor Roger Detels