A numerical approach must be taken to solve for the steepness parameter of a logistic decay when minmax normalization is being applied. The provided function takes the following inputs:
These inputs are then used to solve numerically for the steepness \(k\). We coerce point \((n_c)^{th}\) or \((n_c + 1)^{th}\) to weight \(w_b\) to ensure accuracy of the model at that specific point. The rest of the model then follows the existing logistic pattern subject to this specified constraint.
Consider the input parameters \(n_c \ \epsilon \ [1, 2, 3, ...]\) and \(w_b \ \epsilon \ [0,1] \ni w_b \neq 0.5\). To solve for k, we’ll need to consider three cases.
1. The Case: \(w_b < 0.5 \ , \ n_c = 1, 2, ...\)
Selecting \(w_b < 0.5\) indicates that we are specifying the boundary weight with regard to the \((n_c + 1)^{th}\) point of DUC \(c\) in a service area. Solving will be completed by coercing the \((n_c + 1)^{th}\) point to have weight \(w_b\). Because \(n_c \nless 1\), we can always achieve a reasonable solution for \(k\) at the \((n_c + 1)^{th}\) point despite the use minmax normalization, because the \((n_c + 1)^{th}\) rank will never be a boundary point (i.e., will never be forced to 0 or 1). We solve the following equation numerically to find \(k\).
\(w_b = \frac{\frac{1}{1+e^{k(0.5)}} - \frac{1}{1+e^{k(z - n_c - 0.5)}}}{\frac{1}{1+e^{k(0.5 - n_c)}} - \frac{1}{1+e^{k(z - n_c - 0.5)}}}\)
Where \(w_b\), \(k\), and \(n_c\) are defined as above, and \(z\) is the number of points on which the function is built (in order to guarantee that all points in a service area will have a non-zero weight; \(z = 100\) is the default value, but any large number will suffice and have minimal impact on the model)
Regardless of the value of \(n_c\), when a user selects \(w_b : w_b < 0.5\), output will indicate that \(k\) will be estimated by coercing the \((n_c + 1)^{th}\) weight to \(w_b\).
2. The Case: \(w_b > 0.5 \ , \ n_c = 1\)
Selecting \(w_b > 0.5\) indicates that we are specifying the boundary weight with regard to the \((n_c)^{th}\) point of DUC \(c\) in a service area. Now that we are dealing with the \((n_c)^{th}\) rank itself, note that we cannot always achieve a reasonable solution at rank \(n_c\) because the case of minmax normalization; specifically, we cannot achieve weight \(w_b\) when \(n_c = 1\), because the minmax normalization guarantees that the closest point will always have point closeness weight 1!
When \(w_b > 0.5\) and \(n_c = 1\) are input together, it will be considered an error – output will indicate that this sort of coercion is impossible due to the minmax normalization, and that solving will continue by coercing point \(n = 2\) to weight \(w_n = 1 - w_b\). So, the function automatically turns this case into Case 1, and the same equation as above will be solved.
\(w_b = \frac{\frac{1}{1+e^{k(0.5)}} - \frac{1}{1+e^{k(z - n_c - 0.5)}}}{\frac{1}{1+e^{k(0.5 - n_c)}} - \frac{1}{1+e^{k(z - n_c - 0.5)}}}\)
3. The Case: \(w_b > 0.5 \ , \ n_c > 2\)
When \(n_c > 2\), we are guaranteed that \(n_c\) will not be a boundary point, so we can achieve a reasonable solution under minmax normalization by coercing the point \(n_c\) to have weight \(w_b\). The equation to solve now changes slightly since we are not considering the \((n_c+1)\) case.
\(w_b = \frac{\frac{1}{1+e^{k(-0.5)}} - \frac{1}{1+e^{k(z - n_c - 0.5)}}}{\frac{1}{1+e^{k(0.5 - n_c)}} - \frac{1}{1+e^{k(z - n_c - 0.5)}}}\)
Output will indicate that \(k\) will be estimated by coercing the \((n_c)^{th}\) weight to \(w_b\).
In practice, the Steepness Parameter Function will be wrapped within the Point Closeness Weighting Function. In addition to the information already provided by the Point Closeness Weighting Function, the Steepness Parameter Function will print information on:
The function itself is written using the package rPython in R, which allows interaction between Python and R. This function is coded in Python language due to enhanced precision - the fsolve function in Python’s scipy.optimize package is able to achieve significantly higher precision than any equation solving package in R. Still, it is ultimately implemented in R to maintain the already-established structure of the Category Weighting Function.
If interested, the function code is provided below:
k_est = function(fl = 100, dn = 1, bw = 0.1, ...){
# fl is the same as z above. dn is the same as n_c above. bw is the same as w_b above.
# Assign the variables to a Python environment
python.assign("l", fl)
python.assign("n_c", dn)
python.assign("p", bw)
# Solve for k numerically in Python
python.exec('
import numpy
from scipy.optimize import fsolve
if p < 0.5:
print("Solve for k by coercing the weight for point with closeness rank " + str(n_c + 1) + " to " + str(p))
func = lambda k : (1/(1+numpy.exp(k*(n_c + 1 - n_c - 0.5))) - 1/(1+numpy.exp(k*(l - n_c - 0.5)))) / (p*(1/(1+numpy.exp(k*(1 - n_c - 0.5))) - 1/(1+numpy.exp(k*(l - n_c - 0.5))))) - 1
k_init = 2*numpy.log((1 - p)/p)
k_sol = fsolve(func, k_init)
k_sol = k_sol[0]
elif p > 0.5 and n_c != 1:
print("Solve for k by coercing the weight for point with closeness rank " + str(n_c) + " to " + str(p))
func = lambda k : (1/(1+numpy.exp(k*(n_c - n_c - 0.5))) - 1/(1+numpy.exp(k*(l - n_c - 0.5)))) / (p*(1/(1+numpy.exp(k*(1 - n_c - 0.5))) - 1/(1+numpy.exp(k*(l - n_c - 0.5))))) - 1
k_init = 2*numpy.log((1 - p)/p)
k_sol = fsolve(func, k_init)
k_sol = -k_sol[0]
elif p > 0.5 and n_c == 1:
print("Cannot coerce normalization for dn = 1 case to a weight != 1. By default, the 2nd point will be coerced to a weight = " + str(1-p) + ". Please re-check the input parameters if this is an issue")
p = 1-p
func = lambda k: (1 / (1 + numpy.exp(k * (n_c + 1 - n_c - 0.5))) - 1 / (1 + numpy.exp(k * (l - n_c - 0.5)))) / (p * (1 / (1 + numpy.exp(k * (1 - n_c - 0.5))) - 1 / (1 + numpy.exp(k * (l - n_c - 0.5))))) - 1
k_init = 2*numpy.log((1 - p)/p)
k_sol = fsolve(func, k_init)
k_sol = k_sol[0]
')
# Return the solved value of k to R for use in other functions
assign("k", python.get('k_sol'), envir = .GlobalEnv)
# Print the solved value of k
print(paste("Steepness parameter =", k))
}