Numerically Solve for Parameters Characterizing Lognormal RV that's Truncated from Above

Question

I am trying to numerically solve for parameters characterizing a lognormal distribution truncated from above with first moment = mean, second moment = moment_2, and upper = 99th percentile of the untruncated distibution. This involves solving a constrained non-linear system of 2 equations and 2 unknowns. Currently I am having issues solving these equations numerically, and would appreciate any suggestions on the best solving technique to use.

I am trying to do this in R, and hence have expressed the equations as R functions. The desired parameters logmean and logstd can be recovered by solving the following system of 2 equations and 2 unknowns, which relate the first and second moments of the truncated distribution to the parameters logmean and logstd. Note that these equations have the PDF and CDF for a standard normal random variable in them written as dnorm and pnorm respectively.

mean_fun <-function(mean, logmean, logstd, upper = .99){

    upper_lim <- qlnorm(upper, logmean, logstd)
    b_0 <- (log(upper_lim) - logmean)/logstd
    return(mean - exp(logmean+logstd^2/2)*pnorm(-logstd+b_0)/pnorm(b_0))
}

second_moment_fun <-function(moment_2, logmean, logstd, upper = .99){
     upper_lim <- qlnorm(upper, logmean, logstd)
     b_0 <- (log(upper_lim) - logmean)/logstd
     return(moment_2 -exp(2*logmean+2*logstd^2)*pnorm(-2*logstd+b_0)/pnorm(b_0))
}

UPDATE: This can be consistently solved using R's nleqslv and by providing the jacobian. These equations and their derivation can also be found at this paper by searching for Example 21.75.

Answer 1

Here is code that uses nleqslv to solve this problem. It requires defining the jacobian, which is why the code is somewhat lengthy. Here is the code in full, which includes a test example at the bottom.

library("nleqslv")



mean_fun <-function(mean, logmean, logstd, upper = .99){

    upper_lim <- qlnorm(upper, logmean, logstd)
    b_0 <- (log(upper_lim) - logmean)/logstd
    return(mean - exp(logmean+logstd^2/2)*pnorm(-logstd+b_0)/pnorm(b_0))
}

moment2_fun <-function(moment_2, logmean, logstd, upper = .99){
    upper_lim <- qlnorm(upper, logmean, logstd)
    b_0 <- (log(upper_lim) - logmean)/logstd
    return(moment_2 -exp(2*logmean+2*logstd^2)*pnorm(-2*logstd+b_0)/pnorm(b_0))
}

d_mean_du <- function(mean, logmean, logstd, upper = .99){
    upper_lim <- qlnorm(upper, logmean, logstd)
    b_0 <- (log(upper_lim) - logmean)/logstd
    -exp(logmean+logstd^2/2)*(pnorm(-logstd+b_0)/pnorm(b_0)+
                           -dnorm(-logstd+b_0)/(logstd*pnorm(b_0))+
                           pnorm(-logstd+b_0)*dnorm(b_0)/(logstd*pnorm(b_0)^2)
    )
}

d_mean_dstd <- function(mean, logmean, logstd, upper = .99){
    upper_lim <- qlnorm(upper, logmean, logstd)
    b_0 <- (log(upper_lim) - logmean)/logstd
    -exp(logmean+logstd^2/2)*(pnorm(-logstd+b_0)*logstd/pnorm(b_0)+
                           dnorm(-logstd+b_0)*(-1-b_0/logstd)/pnorm(b_0)+
                           pnorm(-logstd+b_0)*dnorm(b_0)*b_0/(logstd*pnorm(b_0)^2))
}


d_moment2_du <- function(moment2, logmean, logstd, upper = .99){
    upper_lim <- qlnorm(upper, logmean, logstd)
    b_0 <- (log(upper_lim) - logmean)/logstd
    -exp(2*logmean+2*logstd^2)*(pnorm(-2*logstd+b_0)*2/pnorm(b_0)+
                             -dnorm(-2*logstd+b_0)/(logstd*pnorm(b_0))+
                             pnorm(-2*logstd+b_0)*dnorm(b_0)/(logstd*pnorm(b_0)^2)
    )
}


d_moment2_dstd <- function(moment2, logmean, logstd, upper = .99){
    upper_lim <- qlnorm(upper, logmean, logstd)
    b_0 <- (log(upper_lim) - logmean)/logstd
    -exp(2*logmean+2*logstd^2)*(pnorm(-2*logstd+b_0)*4*logstd/pnorm(b_0)+
                             -dnorm(-2*logstd+b_0)*(-2-b_0/logstd)/pnorm(b_0)+
                             pnorm(-2*logstd+b_0)*dnorm(b_0)*b_0/(logstd*pnorm(b_0)^2)
    )
}


trunc_log_norm_mean_moment2 <-function(x, mean, moment2, upper){
    c(mean_fun(mean, x[1], x[2], upper), moment2_fun(moment2, x[1], x[2], upper))
}


trunc_log_norm_jacob <- function(x, mean, moment2, upper){
    matrix(c(d_mean_du(mean, x[1], x[2], upper),
             d_moment2_du(moment2, x[1], x[2], upper),
             d_mean_dstd(mean, x[1], x[2], upper),
             d_moment2_dstd(moment2, x[1], x[2], upper)),2)

}




mean1 <- 10
moment21 <- 20+10^2
log_normal_m_upper_param_gen <-function(mean, std){
    # takes vector of means and stds and returns inputs needed
    # to enter into R's plnorm and dlnorm functions so that
    # these have the desired means. Returns data.frame with named columns
    # as output
    logstd <- sqrt(log((std/mean)^2+1))
    logmean <- log(mean)-.5*logstd^2
    c(logmean, logstd)
}



initial_guess <- log_normal_m_upper_param_gen(mean1, moment21-mean1^2)
upper <- .99

w <- nleqslv(x=initial_guess, fn=trunc_log_norm_mean_moment2, jac =trunc_log_norm_jacob,
        mean = mean1, moment2 = moment21, upper = upper)

p <- c(qlnorm(.99,w$x[1], w$x[2] ), w$x[1], w$x[2])
mean <- integrate(function(m) m*dlnorm(m, p[2], p[3]/plnorm(p[1], p[2], p[3])), 0, p[1])
moment2 <- integrate(function(m) m^2*dlnorm(m, p[2], p[3]/plnorm(p[1], p[2], p[3])), 0, p[1])