From d68894d9edb8de8d3cf9d9669bbb32f841f84505 Mon Sep 17 00:00:00 2001 From: ammeir2 Date: Wed, 17 Dec 2014 12:19:24 +0200 Subject: [PATCH 1/4] Replaced Line Search with Newton Raphson The previous version of the ML algorithm was implemented with a line search algorithm which works well for estimating correlations but may become unstable for general Normal mean estimation (especially of the mean to be estimated is outside of the line search range). The current implementation uses a Newton Raphson approach which works well for all sigma/cutoff/x values. --- selectiveCI/R/pointEstimates.R | 38 ++++++++++++++++++++++++++++++---- 1 file changed, 34 insertions(+), 4 deletions(-) diff --git a/selectiveCI/R/pointEstimates.R b/selectiveCI/R/pointEstimates.R index 34d7ccc..540fc97 100644 --- a/selectiveCI/R/pointEstimates.R +++ b/selectiveCI/R/pointEstimates.R @@ -1,12 +1,42 @@ -ConditionalPointML <- function(x,sigsq,cutoff) { +ConditionalPointML <- function(x,cutoff,sig,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(abs(x) Date: Wed, 17 Dec 2014 12:25:42 +0200 Subject: [PATCH 2/4] x can now be a vector made a small change in order to allow x to be a vector. In which case the mean will be computed and the standard deviation will be divided by sqrt(n). --- selectiveCI/R/pointEstimates.R | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/selectiveCI/R/pointEstimates.R b/selectiveCI/R/pointEstimates.R index 540fc97..942212d 100644 --- a/selectiveCI/R/pointEstimates.R +++ b/selectiveCI/R/pointEstimates.R @@ -15,6 +15,10 @@ ConditionalPointML <- function(x,cutoff,sig,init.mu=NA,max.iter=1000,tol=10^-5) 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 <- sig/sqrt(length(x)) + if(abs(x) Date: Wed, 17 Dec 2014 12:27:41 +0200 Subject: [PATCH 3/4] replaced sig with sigsq replaced the standard deviation with the variance in the function input. --- selectiveCI/R/pointEstimates.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/selectiveCI/R/pointEstimates.R b/selectiveCI/R/pointEstimates.R index 942212d..23c91f5 100644 --- a/selectiveCI/R/pointEstimates.R +++ b/selectiveCI/R/pointEstimates.R @@ -1,4 +1,4 @@ -ConditionalPointML <- function(x,cutoff,sig,init.mu=NA,max.iter=1000,tol=10^-5) { +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) @@ -17,7 +17,7 @@ ConditionalPointML <- function(x,cutoff,sig,init.mu=NA,max.iter=1000,tol=10^-5) ## If x is a vector then replace x with the mean and the variance with the variance of the mean x <- mean(x) - sig <- sig/sqrt(length(x)) + sig <- sqrt(sigsq/length(x)) if(abs(x) Date: Wed, 17 Dec 2014 12:39:17 +0200 Subject: [PATCH 4/4] Parameters and reference Added the extra parameters of the function and added the reference to the arxiv paper. --- selectiveCI/man/ConditionalPointML.Rd | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/selectiveCI/man/ConditionalPointML.Rd b/selectiveCI/man/ConditionalPointML.Rd index f08c083..b25ffaa 100644 --- a/selectiveCI/man/ConditionalPointML.Rd +++ b/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