有winbugs的高手，帮忙看看这个文献，感激不尽

> ## Avoid printing to unwarranted accuracy
> od <- options(digits = 5)
> x <- 0:10
> y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
>
> ## Easy one-dimensional MLE:
> nLL <- function(lambda) -sum(stats::dpois(y, lambda, log=TRUE))
> fit0 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y))
错误: 没有"mle"这个函数
> # For 1D, this is preferable:
> fit1 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y),
+             method = "Brent", lower = 1, upper = 20)
错误: 没有"mle"这个函数
> stopifnot(nobs(fit0) == length(y))
错误于nobs(fit0) : 找不到对象'fit0'
>
> ## This needs a constrained parameter space: most methods will accept NA
> ll <- function(ymax = 15, xhalf = 6) {
+     if(ymax > 0 && xhalf > 0)
+       -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE))
+     else NA
+ }
> (fit <- mle(ll, nobs = length(y)))
错误: 没有"mle"这个函数
> mle(ll, fixed = list(xhalf = 6))
错误: 没有"mle"这个函数
> ## alternative using bounds on optimization
> ll2 <- function(ymax = 15, xhalf = 6)
+     -sum(stats::dpois(y, lambda=ymax/(1+x/xhalf), log = TRUE))
> mle(ll2, method = "L-BFGS-B", lower = rep(0, 2))
错误: 没有"mle"这个函数
>
> AIC(fit)
错误于AIC(fit) : 找不到对象'fit'
> BIC(fit)
错误于BIC(fit) : 找不到对象'fit'
>
> summary(fit)
错误于summary(fit) : 找不到对象'fit'
> logLik(fit)
错误于logLik(fit) : 找不到对象'fit'
> vcov(fit)
错误于vcov(fit) : 找不到对象'fit'
> plot(profile(fit), absVal = FALSE)
错误于profile(fit) : 找不到对象'fit'
> confint(fit)
错误于confint(fit) : 找不到对象'fit'
>
> ## Use bounded optimization
> ## The lower bounds are really > 0,
> ## but we use >=0 to stress-test profiling
> (fit2 <- mle(ll, method = "L-BFGS-B", lower = c(0, 0)))
错误: 没有"mle"这个函数
> plot(profile(fit2), absVal=FALSE)
错误于profile(fit2) : 找不到对象'fit2'
>
> ## a better parametrization:
> ll3 <- function(lymax = log(15), lxhalf = log(6))
+     -sum(stats::dpois(y, lambda=exp(lymax)/(1+x/exp(lxhalf)), log=TRUE))
> (fit3 <- mle(ll3))
错误: 没有"mle"这个函数
> plot(profile(fit3), absVal = FALSE)
错误于profile(fit3) : 找不到对象'fit3'
> exp(confint(fit3))
错误于confint(fit3) : 找不到对象'fit3'
>
> options(od)
>
epoh老师，您好·！
这个MLE函数在R的什么地方？如何安装？

epoh 发表于 2011-12-31 08:11
R package "POT"
提供了多种估计 GPD 参数的方法
method of moments moments,

好的，实在太感谢了，我先仔细看看，有问题恐怕还得麻烦你喽

zhangtao 发表于 2011-12-31 17:09
> ## Avoid printing to unwarranted accuracy
> od  x  y  
> ## Easy one-dimensional MLE:

library(stats4)
http://stat.ethz.ch/R-manual/R-devel/library/stats4/html/mle.html
请注意短信息

epoh 发表于 2011-12-31 08:11
R package "POT"
提供了多种估计 GPD 参数的方法
method of moments moments,

您好！首先很感谢您的指导。不过我还存在一些问题，因为我主要是想用文献中MCMC对GPD进行参数估计，但是R中的没有直接用此方法的，所以我看到用Winbugs,不知道openbugs是不是用的MCMC法，openbugs和winbugs一样的吗？实在不好意思，问这样弱智问题，因为我从来没接触过这样的软件，让您笑话啦

sky-xiaotian 发表于 2011-12-31 19:40
您好！首先很感谢您的指导。不过我还存在一些问题，因为我主要是想用文献中MCMC对GPD进行参数估计，但是R ...

Openbugs & winbugs
两者极为相同,
但openbugs含GPD
下载了你就知道
http://openbugs.info/w/Downloads

epoh 发表于 2011-12-31 20:29
Openbugs & winbugs
两者极为相同,
但openbugs含GPD

谢谢，新年快乐哈

sky-xiaotian 发表于 2011-12-31 20:45
谢谢，新年快乐哈

安装完毕
请看一下help\Distributions
Generalized Pareto
x ~ dgpar(mu,sigma,eta)
你只要设下sigma & eta prior
即可运行
也祝你 新年快乐!

补上以你的数据做出的结果

            mean    sd  MC_errorval2.5pcmedian  val97.5pc    start  sample

sigma   99.69   22.6    0.5598    60.89      97.43   149.9    10001   10000

xi       0.3189  0.2027  0.006802 6.816E-40.2985  0.7821  10001   10000

epoh 发表于 2011-12-31 20:55
安装完毕
请看一下help\Distributions
Generalized Pareto

您好！已经运行了一下，有结果出来，但是结果跟事实有点不符
        mean        sd        MC_error        val2.5pc        median        val97.5pc        starts       ample

eta        1.182E-10        3.01E-9        8.422E-11        1.0E-20        1.013E-20        1.086E-20        501        29500
sigma        0.0789        0.002152        3.72E-5        0.0763        0.07838        0.08449        10001        20000

        mean        sd        MC_error        val2.5pc        median        val97.5pc        starts       ample
sigma        1.948        15.65        1.19        0.07631        0.07841        0.08609        501        24500

每一次试验结果还都不一样，我不知道应该如何确定，还请老师指导下；
还有一个问题，sigma和eta的初始分布、初始值是随意设置的吗？不同的设置会不会影响结果呢？

附我给出的初始分布：

model
{
        for (i in 1:N)
                {
                        x ~ dgpar(mu,sigma,eta)
                }
        eta~ dgamma(0.01,0.01)
               sigma~ dpar(0.01,0.01)
}
data
list(N=53,mu=308,x=c(309,310,311,313,313,316,319,320,328,330,332,332,338,341,341,349,351,354,356,365,368,370,371,375,384,384,385,389,389,391,396,408,417,418,418,435,447,448,460,463,465,473,500,507,521,522,554,618,702,754,842,858,1152))
        initial values
        list(eta=1,sigma=1)

谢谢！又麻烦您了

sky-xiaotian 发表于 2012-1-1 09:41
您好！已经运行了一下，有结果出来，但是结果跟事实有点不符
        mean        sd        MC_error        val2.5pc        median        val97 ...

哈哈,答案昨天就给你了
请注意短信息

老师，您好，关于我那篇文献，第一步（Generic Leak Frequencies (Phase 1)）已经求出，那如何再求Hydrogen Leak Frequencies (Phase 2)呢，需要重新建方程吗，还是利用Equation 4.3，我看了几天，还是没有思路，麻烦老师再指点，谢谢！