假设我有2个data.frame对象:
df1 <- data.frame(x = 1:100)
df1$y <- 20 + 0.3 * df1$x + rnorm(100)
df2 <- data.frame(x = 1:200000)
df2$y <- 20 + 0.3 * df2$x + rnorm(200000)
我想做MLE.随着df1一切正常:
LL1 <- function(a, b, mu, sigma) {
R = dnorm(df1$y - a- b * df1$x, mu, sigma)
-sum(log(R))
}
library(stats4)
mle1 <- mle(LL1, start = list(a = 20, b = 0.3, sigma=0.5),
fixed = list(mu = 0))
> mle1
Call:
mle(minuslogl = LL1, start = list(a …我有以下代码，希望估算自定义分发的参数。有关发行的更多详细信息。然后，使用估计的参数，我想查看估计的PDF是否类似于给定数据的分布（应该与给定数据的分布匹配）。
[编辑]：“ x”现在保存数据样本，而不是PDF
主要代码是：
x = [0.0320000000000000 0.0280000000000000 0.0280000000000000 0.0270000000000000 0.0320000000000000 0.0320000000000000 0.0480000000000000 0.0890000000000000 0.0500000000000000 0.0620000000000000 0.0480000000000000 0.0300000000000000 0.0520000000000000 0.0460000000000000 0.0540000000000000 0.0520000000000000 0.0510000000000000 0.0310000000000000 0.0330000000000000 0.0330000000000000 0.0380000000000000 0.0850000000000000 0.102000000000000 0.0290000000000000 0.0530000000000000 0.0590000000000000 0.0320000000000000 0.0800000000000000 0.0410000000000000 0.0280000000000000 0.0670000000000000 0.0350000000000000 0.0420000000000000 0.0280000000000000 0.0370000000000000 0.0480000000000000 0.0330000000000000 0.101000000000000 0.0420000000000000 0.0840000000000000 0.0340000000000000 0.0900000000000000 0.0900000000000000 0.0460000000000000 0.0290000000000000 0.0330000000000000 0.0350000000000000 0.0330000000000000 0.0320000000000000 0.0420000000000000 0.0600000000000000 0.0500000000000000 0.0390000000000000 0.0480000000000000 0.0680000000000000 0.0330000000000000 0.0510000000000000 0.0430000000000000 0.0270000000000000 0.0330000000000000 0.0590000000000000 0.0380000000000000 0.0270000000000000 0.0600000000000000 0.0310000000000000 0.0520000000000000 0.0350000000000000 0.0640000000000000 0.0570000000000000 0.0520000000000000 0.0330000000000000 0.0480000000000000 …我有一个简单的x, ydata.frame.
mydata <- data.frame(days = 1:96, risk = c(5e-09, 5e-09, 5e-09, 1e-08, 4e-08, 6e-08, 9e-08, 1.5e-07, 4.2e-07,
7.2e-07, 1.02e-06, 1.32e-06, 1.66e-06, 2.19e-06, 2.76e-06, 3.32e-06,
3.89e-06, 4.55e-06, 5.8e-06, 7.16e-06, 8.51e-06, 9.85e-06, 1.138e-05,
1.396e-05, 1.672e-05, 1.947e-05, 2.222e-05, 2.521e-05, 2.968e-05,
3.439e-05, 3.909e-05, 4.378e-05, 4.894e-05, 5.697e-05, 6.546e-05,
7.392e-05, 8.236e-05, 9.16e-05, 0.00010573, 0.00012063, 0.00013547,
0.00015025, 0.00016642, 0.00019127, 0.00021743, 0.00024343, 0.00026924,
0.00029818, 0.00034681, 0.00039832, 0.00044932, 0.00049976, 0.0005451,
0.00056293, 0.00057586, 0.00058838, 0.0006005, 0.00061562, 0.00065079,
0.00068845, 0.00072508, 0.00076062, 0.00079763, 0.00084886, 0.00090081,
0.0009507, 0.00099844, 0.00104427, 0.00108948, …我正在对气象数据进行极值分析，以准确计算以 mm/d 为单位的可用降水数据。我正在使用阈值过量方法来估计具有最大似然法的广义帕累托分布的参数。
目的是计算日降水量的几个回归水平（即 2、5、10、20、50、100 年事件）。
虽然 R 代码工作正常，但我想知道为什么在根据具有不同包的拟合 GPD 的分位数计算回报水平时，我得到明显不同的结果。尽管每个包中 GPD 的估计参数几乎相同，但分位数差异很大。
我使用的软件包是：ismev、extRemes、evir 和 POT。
我猜对 GPD 参数的不同估计是由于不同的计算程序，但我不明白为什么分位数的计算会因不同的包而有很大差异。
虽然 lmom、evir 和 POT 返回相同的分位数值，但从 extRemes 包导出的返回期与其他结果不同。
# packages
library(ismev)
library(extRemes)
library(evir)
library(POT)
library(lmom)
th <- 50
# sample data:
potvalues <- c(
58.5,44.2,49.6,59.3,48.3,60.9,94.5,47.1,45.3,57.6,48.2,46.2,44.2,50.6,42.1,52.7,80.9,
58.5,51.3,48.4,51.7,71.9,60.1,64.4,43.5,55.5,49.3,58.2,47.5,43.7,45.2,52.8,42.2,46.4,
96.1,47.5,50.1,42.4,60.9,72.6,51.6,59.4,80.5,63.7,59.9,45.0,66.7,47.6,53.3,43.1,51.0,
46.2,53.6,59.8,51.7,46.7,42.6,44.5,45.0,50.0,44.0,89.9,44.2,47.8,53.3,43.0,55.7,44.6,
44.6,54.9,45.1,43.9,78.7,45.5,64.0,42.7,47.4,57.0,105.4,64.3,43.2,50.4,80.2,49.9,71.6,
47.4,44.1,47.6,55.2,44.4,78.6,50.8,42.4,47.1,43.5,51.4)
#------------------------------------------------------------------------------------------#
# MLE Fitting of GPD - package extRemes
# fit gpd
pot.ext <- fevd(potvalues, method = "MLE", type="GP", threshold=th)
# return levels:
rl.extremes <- return.level(pot.ext, conf = 0.05,
return.period= c(2,5,10,20,50,100))
rl.extremes …我的朋友目前正致力于他的关于时间序列模型参数估计的任务,SARIMAX(季节性ARIMA外生),具有最大似然估计(MLE)方法.他使用的数据是关于2000 - 2012年的月降雨量,印度洋偶极子(IOD)指数作为外生变量.这是数据:
MONTH YEAR RAINFALL IOD
1 1 2000 15.3720526 0.0624
2 2 2000 10.3440804 0.1784
3 3 2000 14.6116392 0.3135
4 4 2000 18.6842179 0.3495
5 5 2000 15.2937896 0.3374
6 6 2000 15.0233152 0.1946
7 7 2000 11.1803399 0.3948
8 8 2000 11.0589330 0.4391
9 9 2000 10.1488916 0.3020
10 10 2000 21.1187121 0.2373
11 11 2000 15.3980518 -0.0324
12 12 2000 18.9393770 -0.0148
13 1 2001 19.1075901 -0.2448
14 2 2001 14.9097284 0.1673
15 3 2001 …我需要使用MLE从样本数据中估计Weibull分布的形状和比例参数.我检查了Apache commons Math的WeibullDistribution类,但它没有这样的功能.有什么建议吗？我需要在我的java应用程序中使用该类.
我想知道是否有人知道一个 python 包，它实现 MLE 来估计狄利克雷分布的参数。
我试图使用MLE拟合多元正态分布的参数.
import numpy as np
from scipy.stats import norm
from scipy.optimize import minimize
from scipy.stats import multivariate_normal as mnorm
def estimation(obs,fun,init,method='Nelder-Mead'):
mle = lambda param: -np.sum(fun(*[obs,param])) ## negate since we will minimize
result = minimize(mle,init,method=method)
return result.x
拟合单变量正态分布是好的:
obs = np.random.normal(1,4,50000)
ini = [0,1]
print(estimation(obs,lambda ob,p:norm.logpdf(ob,p[0],p[1]),ini))
但是遇到了多变量的一些问题(错误将数组赋值给变量):
obs_m = np.random.multivariate_normal([0,0],[[1,0],[0,100]],50000)
ini_m = [[0,0],[[1,0],[0,100]]]
print(estimation(obs_m,lambda ob,p:mnorm.logpdf(ob,p[0],p[1],ini_m))
似乎优化算法不适用于任意数组/矩阵.我必须将平均数组和协方差矩阵打包成一个平面阵列,以便最小化.
ini_m = [0,0,1,0,0,100]
print(estimation(obs_m,lambda ob,p:mnorm.logpdf(ob,[p[0],p[1]],[[p[2],p[3]],[p[4],p[5]]]),ini_m))
显然,当尺寸增加时,这将很快失控,或者在没有封闭形式解决方案的情况下,这将更加复杂.这里最好做什么？谢谢.
我正在努力解决以下问题。简而言之：两个不同的软件包（Aptech 的 Gauss 和 R）在最大似然法中产生完全不同的 Hessian 矩阵。我使用相同的程序（BFGS），完全相同的数据，相同的最大似然公式（这是一个非常简单的 logit 模型），具有完全相同的起始值，令人困惑的是，我得到了相同的参数和 log-结果可能性。两个程序中只有 Hessian 矩阵不同，因此标准误差的估计和统计推断不同。
在这个特定的例子中，它并没有出现太大的偏差，但是模型的每一个增加的复杂性都会增加差异，所以如果我尝试估计我的最终模型，两个程序都会产生完全错误的结果。
有谁知道，这两个程序在计算 Hessian 矩阵的方式上有何不同，以及获得相同结果的正确方法可能是什么？
编辑：在 R （高斯）代码中，向量X ( alt ) 是自变量，由两列向量组成，第一列完全是 1，第二列是受试者的响应。向量y ( itn ) 是因变量，由一列和受试者的响应组成。该示例（R 代码和数据集）取自http://www.polsci.ucsb.edu/faculty/glasgow/ps206/ps206.html，只是作为重现和隔离问题的示例。
我已附上代码（高斯和 R 语法）和输出。
任何帮助将不胜感激。谢谢 ：）
高斯：
start={ 0.95568840 , -0.20459156 };
library maxlik,pgraph;
maxset;
_max_Algorithm = 2;
_max_Diagnostic = 1;
{betaa,f,g,cov,ret} = maxlik(XMAT,0,&ll,start);
call maxprt(betaa,f,g,cov,ret);
print _max_FinalHess;
proc ll(b,XMAT);
local exb, probo, logexb, yn, logexbn, yt, ynt, logl;
exb = EXP(alt*b);
//print exb;
probo = exb./(1+exb); …我想将权重合并到做权重的可能性中svyglm。
据杰里米万里和其他地方，svyglm功能使用权，以“体重每起案件的重要性，使他们代表（给对方，...）”。
这是我的数据：
(dat <- data.frame(
A = c(1, 1, 0, 0), B = c(1, 0, 1, 0),
Pass = c(278, 100, 153, 79), Fail = c(743, 581, 1232, 1731), Weights= c(3, 1, 12, 3)
))
这是我的likelihood功能：
ll <- function (b0, b1, b2, b3) {
odds <- exp(b0) * (1 + b1 * dat$A + b2 * dat$B + b3 * dat$A * dat$B)
-sum(dbinom(
x = dat$Pass, size = …mle ×10
r ×6
estimation ×2
python ×2
statistics ×2
dirichlet ×1
distribution ×1
fminsearch ×1
gauss ×1
java ×1
matlab ×1
quantile ×1
regression ×1
scipy ×1
time-series ×1
weibull ×1