【发布时间】：2022-01-22 10:46:17
【问题描述】：
我正在尝试将用于从 stats 包中进行非线性最小化的 nlm 函数从 R 转换为 Python。 在 R 的帮助菜单中，它说“函数使用牛顿型算法执行函数 f 的最小化”
R中原正确的代码如下：
m=3
la1 = 1
la2 = 3
la3 = 7
del = matrix(0,1,m-1);del
del[1] = 0.5
del[2] = 0.2
del[3] = 1 - (del[1]+del[2])
del
cd = cumsum(del);cd
# Generate random data for a m mixture model
N = 100 # number of data points
unifvec = runif(N)
d1 = rpois(sum(unifvec < cd[1]),la1);d1
d2 = rpois(sum(unifvec > cd[1] & unifvec < cd[2]),la2)
d3 = rpois(sum(unifvec > cd[2]),la3);d3
data = c(d1,d2,d3);data # Data vector
# Functions for parameter transformation
logit <- function(vec) log(vec/(1-sum(vec)))
invlogit <- function(vec) exp(vec)/(1+sum(exp(vec)))
# Make function for the negative log-likelihood
f <- function(PAR) {
M = length(PAR)
m = ceiling(M/2)
LA = exp(PAR[1:m]) # transform lambdas
DELs = invlogit(PAR[(m+1):M]) # transform deltas
DEL = c(DELs,1-sum(DELs))
# Equation (1.1) on p. 9
L = DEL[1]*dpois(data,LA[1])
for(i in 2:m){
L = L+DEL[i]*dpois(data,LA[i])
}
-sum(log(L))
}
# Define starting guess for optimization
par = c(2,4,7,0.5,0.2)
PAR = par
PAR[1:3] = log(par[1:3])
PAR[4:5] = logit(par[4:5])
PAR
# Optimize using nlm
res = nlm(f,PAR);res
结果：
$minimum
[1] 211.2481
$estimate
[1] -0.2827635 1.1449476 1.8346681 0.7380511 -0.3305408
$gradient
[1] -1.185185e-05 -1.372745e-05 -3.084352e-05 4.720846e-05
[5] 1.506351e-06
$code
[1] 1
$iterations
[1] 22
在 Python 中这样做我发现 scipy.optimize.minimize 我认为它做同样的事情。
m=3
la1 = 1
la2 = 3
la3 = 7
de = np.zeros(m);de
de[0] = 0.5
de[1] = 0.2
de[2] = 1 - (de[0]+de[1])
de
cd = np.cumsum(de);cd
N = 100 # number of data points
unifvec = np.random.uniform(0,1,N)
d1 = np.random.poisson(la1, sum(unifvec < cd[0], la1));d1
d2 = np.random.poisson(la2,sum( (unifvec > cd[0]) & (unifvec <cd[1]) ) );d2
d3 = np.random.poisson(la1, sum(unifvec < cd[1], la1));d3
data = np.concatenate((d1,d2,d3), axis=None);(data)
# Functions for parameter transformation
def logit(vec):
return(np.log(vec/(1-sum(vec))))
def invlogit(vec):
return(np.exp(vec)/(1+sum(np.exp(vec))))
# Make function for the negative log-likelihood
def f(PAR):
M = len(PAR)
m = int(np.ceil(M/2))
LA = np.exp(PAR[:m]) # transform lambdas
DELs = invlogit(PAR[m:M]) # transform deltas
DEL = np.concatenate((DELs, 1-sum(DELs)), axis=None)
# Equation (1.1) on p. 9
from scipy.stats import poisson
L = DEL[0]*poisson.pmf(data,LA[0])
for i in range(1,m):
L = L+DEL[i]*poisson.pmf(data,LA[i])
return(-sum(np.log(L)))
# Define starting guess for optimization
par = np.array([2,4,7,0.5,0.2])
PAR = par
PAR[0:3] = np.log(par[:3])
PAR[3:5] = logit(par[3:5])
from scipy.optimize import minimize
res = minimize(f, PAR, method='Nelder-Mead', tol=1e-6)
res
但结果与R中的不相似
final_simplex: (array([[ 0.29287179, 1.87973654, -175.80084603, 3.0604109 ,
-0.95479 ],
[ 0.29287179, 1.87973654, -175.80084639, 3.06041089,
-0.95479 ],
[ 0.29287179, 1.87973654, -175.80084608, 3.0604109 ,
-0.95479 ],
[ 0.29287179, 1.87973654, -175.80084587, 3.0604109 ,
-0.95479 ],
[ 0.29287179, 1.87973654, -175.80084553, 3.0604109 ,
-0.95478999],
[ 0.29287179, 1.87973654, -175.80084576, 3.0604109 ,
-0.95478999]]), array([211.89342437, 211.89342437, 211.89342437, 211.89342437,
211.89342437, 211.89342437]))
fun: 211.89342437364002
message: 'Optimization terminated successfully.'
nfev: 782
nit: 457
status: 0
success: True
x: array([ 0.29287179, 1.87973654, -175.80084603, 3.0604109 ,
-0.95479 ])
最小值相似，但估计值不同。特别是第三个估计值是完全错误的
【问题讨论】：
-
1.当我运行
res = nlm(f,PAR)时，它返回 "Error in dpois(data, LA[1]) : Non-numeric argument to math function" 。 2. 您还应该列出您在示例中使用的包。 -
@PeaceWang 我已经用我一开始没有插入的缺失函数更新了这个问题。你现在可以再次运行它，首先运行两个添加的函数
-
不幸的是，我在运行
dpois(data,LA[1])（数学函数的非数字参数）时仍然遇到问题。你的data是什么？ -
抱歉我会再次更新
-
很好。为了在
R和Python（使用np.random）中保持相同的data，您可以手动为data分配一个值，即data = ...（删除如何生成它的过程） ?
标签： python r minimization