Replaced Line Search with Newton Raphson

selectiveCI/R/pointEstimates.R

@@ -1,12 +1,47 @@
-ConditionalPointML <- function(x,sigsq,cutoff) {
+ConditionalPointML <- function(x,cutoff,sigsq,init.mu=NA,max.iter=1000,tol=10^-5) {
+  comp.power <- function(mu,sig,cutoff) 1-pnorm(cutoff,mean=mu,sd=sig)+pnorm(-cutoff,mean=mu,sd=sig)
+  phi <- function(c) dnorm(c,mean=muhat,sd=sig)
+
+  comp.grad <- function(x,mu,sig,cutoff) {
+    Qc <- comp.power(mu,sig,cutoff)
+    (x-mu)/sig^2 - (phi(cutoff)-phi(-cutoff))/Qc
+  }
+  comp.hess <- function(x,mu,sig,cutoff) {
+    Qc <- comp.power(mu,sig,cutoff)
+    hess <- Qc*(phi(cutoff)*(cutoff-mu)/sig^2+phi(-cutoff)*(cutoff+mu)/sig^2) + 
+      (phi(cutoff)-phi(-cutoff))^2
+    hess <- -hess/Qc^2
+    hess <- hess-1/sig^2
+    return(hess)
+  }
+
+  ## If x is a vector then replace x with the mean and the variance with the variance of the mean
+  x <- mean(x)
+  sig <- sqrt(sigsq/length(x))
+
   if(abs(x)<cutoff) {
-    warning("abs(x) is less than the cutoff value!")
+    warning("abs(x) is less than the cutoff value")
+  }
+
+  if(is.na(init.mu)) init.mu <- x
+  muhat <- init.mu
+
+  for(i in 1:max.iter) {
+    old.muhat <- muhat
+
+    grad <- comp.grad(x,muhat,sig,cutoff)
+    hess <- comp.hess(x,muhat,sig,cutoff)
+
+    muhat <- muhat - grad/hess
+
+    if(abs(muhat-old.muhat)<tol) break
   }
-  f <- function(mu,x,cutoff) dnorm(x,mean=mu,sd=sqrt(sigsq))/(1-pnorm(cutoff,mean=mu,sd=sqrt(sigsq))+pnorm(-cutoff,mean=mu,sd=sqrt(sigsq)))
-  muhat <- suppressWarnings(optimize(f,interval=c(-6,6),x,cutoff,maximum=TRUE)[[1]])
+
   return(muhat)
 }
 
+
+
 ##Testing
 # ConditionalPointML(0,1,1)
 

selectiveCI/man/ConditionalPointML.Rd

@@ -8,13 +8,16 @@ Selective point estimation using maximum likelihood.
 Returns a selective maximum likelihood point estimate for the mean of a Gaussian population.
 }
 \usage{
-  ConditionalPointML(x, sigsq, cutoff)
+  ConditionalPointML(x, sigsq, cutoff,init.mu=NA,max.iter=1000,tol=10^-5)
 }
 
 \arguments{
   \item{x}{The observed value from the population.}
   \item{sigsq}{The (known) standard deviation of the population.}
   \item{cutoff}{The selection cutoff.}
+  \item{init.mu}{Starting point for the Newton Raphson algorithm. If NA then the algorithm initializes from the observed value}
+  \item{max.iter}{Maximal number of iteration for the Newton Raphson algorithm.}
+  \item{tol}{The desired accuracy}
 }
 
 %\details{
@@ -23,9 +26,9 @@ Returns a selective maximum likelihood point estimate for the mean of a Gaussian
 
 \value{Returns a numeric value of the estimated mean.}
 
-%\references{
-%% ~put references to the literature/web site here ~
-%}
+\references{
+Benjamini, Yoav, and Amit Meir. "Selective Correlations-the conditional estimators." arXiv preprint arXiv:1412.3242 (2014).‏
+}
 
 \author{
   Amit Meir