Commit 2009be4e authored by Felix Fauer's avatar Felix Fauer
Browse files

first try, handle eta2~~0

parent bcc3dd8e
...@@ -39,8 +39,8 @@ ...@@ -39,8 +39,8 @@
#' lines(x,dens[[i]],lty=i) #' lines(x,dens[[i]],lty=i)
#' } #' }
#' legend('topright',title = 'Duration',legend = 1:4,lty=1:4) #' legend('topright',title = 'Duration',legend = 1:4,lty=1:4)
dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,eta2=NULL,tau=0,...) { dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,eta2=0,tau=0,...) {
if(is.null(eta2)){eta2=eta} #if(is.null(eta2)){eta2=eta}
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta),length(eta2),length(tau))>1)){ if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta),length(eta2),length(tau))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta, tau is a vector. ', message('One of the parameters mut, sigma0, xi, theta, eta, tau is a vector. ',
'This is not intended and might cause an error.')} 'This is not intended and might cause an error.')}
...@@ -52,7 +52,7 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,eta2=NULL,tau=0,...) { ...@@ -52,7 +52,7 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,eta2=NULL,tau=0,...) {
if(d+theta<=0){return(rep(NA,length(q)))}else{ if(d+theta<=0){return(rep(NA,length(q)))}else{
#sigma.d <-sigma0/(d+theta)^eta+ tau # old #sigma.d <-sigma0/(d+theta)^eta+ tau # old
sigma.d <- sigma0/(d+theta)^eta2 +tau sigma.d <- sigma0/(d+theta)^(eta+eta2) +tau
mu.d <- mut*(sigma0/(d+theta)^eta +tau) mu.d <- mut*(sigma0/(d+theta)^eta +tau)
return(dgev(q,loc=mu.d,scale=sigma.d,shape=xi,...))} return(dgev(q,loc=mu.d,scale=sigma.d,shape=xi,...))}
...@@ -92,8 +92,8 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,eta2=NULL,tau=0,...) { ...@@ -92,8 +92,8 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,eta2=NULL,tau=0,...) {
#' @examples #' @examples
#' x <- seq(4,20,0.1) #' x <- seq(4,20,0.1)
#' prob <- pgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1) #' prob <- pgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1)
pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL, ...) { pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,eta2=0, ...) {
if(is.null(eta2)){eta2=eta} #if(is.null(eta2)){eta2=eta}
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta),length(eta2),length(tau))>1)){ if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta),length(eta2),length(tau))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta, tau is a vector. ', message('One of the parameters mut, sigma0, xi, theta, eta, tau is a vector. ',
'This is not intended and might cause an error.')} 'This is not intended and might cause an error.')}
...@@ -105,7 +105,7 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL, ...) { ...@@ -105,7 +105,7 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL, ...) {
if(d+theta<=0){return(rep(NA,length(q)))}else{ if(d+theta<=0){return(rep(NA,length(q)))}else{
#sigma.d <-sigma0/(d+theta)^eta+tau # old #sigma.d <-sigma0/(d+theta)^eta+tau # old
sigma.d <- sigma0/(d+theta)^eta2 +tau sigma.d <- sigma0/(d+theta)^(eta+eta2) +tau
mu.d <- mut*(sigma0/(d+theta)^eta +tau) mu.d <- mut*(sigma0/(d+theta)^eta +tau)
return(pgev(q,loc=mu.d,scale=sigma.d,shape=xi,...))} return(pgev(q,loc=mu.d,scale=sigma.d,shape=xi,...))}
...@@ -158,8 +158,8 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL, ...) { ...@@ -158,8 +158,8 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL, ...) {
#' } #' }
#' legend('topright',title = 'p-quantile', #' legend('topright',title = 'p-quantile',
#' legend = p,lty=1:3,bty = 'n') #' legend = p,lty=1:3,bty = 'n')
qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL, ...) { qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0,eta2=0, ...) {
if (is.null(eta2)){eta2=eta} #if (is.null(eta2)){eta2=eta}
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta), length(eta2), length(tau))>1)){ if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta), length(eta2), length(tau))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta, eta2, tau is a vector. ', message('One of the parameters mut, sigma0, xi, theta, eta, eta2, tau is a vector. ',
'This is not intended and might cause an error.')} 'This is not intended and might cause an error.')}
...@@ -172,7 +172,7 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL, ...) { ...@@ -172,7 +172,7 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL, ...) {
#sigma.d <-sigma0/(d+theta)^eta #sigma.d <-sigma0/(d+theta)^eta
#sigma.d <-sigma0/(d+theta)^eta+tau #sigma.d <-sigma0/(d+theta)^eta+tau
sigma.d <- sigma0/(d+theta)^eta2 +tau sigma.d <- sigma0/(d+theta)^(eta+eta2) +tau
mu.d <- mut*(sigma0/(d+theta)^eta +tau) mu.d <- mut*(sigma0/(d+theta)^eta +tau)
return(qgev(p,loc=as.numeric(mu.d) return(qgev(p,loc=as.numeric(mu.d)
...@@ -220,8 +220,8 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL, ...) { ...@@ -220,8 +220,8 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL, ...) {
#' hist(samp[[2]],breaks = 10,add=TRUE,col=rgb(0,0,1,0.5),freq = FALSE) #' hist(samp[[2]],breaks = 10,add=TRUE,col=rgb(0,0,1,0.5),freq = FALSE)
#' legend('topright',fill = c(rgb(1,0,0,0.5),rgb(0,0,1,0.5)), #' legend('topright',fill = c(rgb(1,0,0,0.5),rgb(0,0,1,0.5)),
#' legend = paste('d=',1:2,'h'),title = 'Duration') #' legend = paste('d=',1:2,'h'),title = 'Duration')
rgev.d <- function(n,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL) { rgev.d <- function(n,mut,sigma0,xi,theta,eta,d,tau=0,eta2=0) {
if (is.null(eta2)){eta2=eta} #if (is.null(eta2)){eta2=eta}
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta),length(eta2),length(tau))>1)){ if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta),length(eta2),length(tau))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta, tau is a vector. ', message('One of the parameters mut, sigma0, xi, theta, eta, tau is a vector. ',
'This is not intended and might cause an error.')} 'This is not intended and might cause an error.')}
...@@ -233,7 +233,7 @@ rgev.d <- function(n,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL) { ...@@ -233,7 +233,7 @@ rgev.d <- function(n,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL) {
if(d+theta<=0){return(rep(NA,n))}else{ if(d+theta<=0){return(rep(NA,n))}else{
#sigma.d <-sigma0/(d+theta)^eta+tau # old #sigma.d <-sigma0/(d+theta)^eta+tau # old
sigma.d <- sigma0/(d+theta)^eta2 +tau sigma.d <- sigma0/(d+theta)^(eta+eta2) +tau
mu.d <- mut*(sigma0/(d+theta)^eta +tau) mu.d <- mut*(sigma0/(d+theta)^eta +tau)
return(rgev(n,loc=mu.d,scale=sigma.d,shape=xi))} return(rgev(n,loc=mu.d,scale=sigma.d,shape=xi))}
......
...@@ -672,8 +672,10 @@ gev.d.params <- function(fit,ydat=NULL){ ...@@ -672,8 +672,10 @@ gev.d.params <- function(fit,ydat=NULL){
eta <- etalink(etmat %*% (fit$mle[seq(npmu + npsc + npsh + npth + 1, length = npet)])) eta <- etalink(etmat %*% (fit$mle[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
if(!fit$eta2_zero){ if(!fit$eta2_zero){
eta2 <- eta2link(e2mat %*% (fit$mle[seq(npmu + npsc + npsh + npth + npet + 1, length = npe2)])) eta2 <- eta2link(e2mat %*% (fit$mle[seq(npmu + npsc + npsh + npth + npet + 1, length = npe2)]))
# transform eta2 from eta2~~eta to eta2~~0
eta2 <- eta2-eta
}else{ }else{
eta2 <- eta #rep(0,dim(ydat)[1]) eta2 <- rep(0,dim(ydat)[1]) #eta
} }
if(!fit$tau_zero){ if(!fit$tau_zero){
tau <- taulink(tamat %*% (fit$mle[seq(npmu + npsc + npsh + npth + npet + npe2 + 1, length = npta)])) tau <- taulink(tamat %*% (fit$mle[seq(npmu + npsc + npsh + npth + npet + npe2 + 1, length = npta)]))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment