mpsbetaexpg<-function(mydata, g, location=TRUE, method , sig.level){
min.mydata<-min(mydata)
sort.mydata<-sort(mydata)
n<-length(mydata)
y<-(sort.mydata[2:n]-min.mydata)
y<-(sort.mydata[2:n]-min.mydata)[is.finite(log(sort.mydata[2:n]-min.mydata))]
n<-length(mydata)
qp1<-y[floor(.25*n)]
qp3<-y[floor(.75*n)]
median.mydata<-y[floor(.5*n)]
mean.mydata<-mean(y)
inv.mydata<-1/y[is.finite(1/y)]
inv.mean<-mean(inv.mydata)
std.mydata<-sd(y)
if(g!="birnbaum-saunders" & g!="exp" & g!="rayleigh" & g!="weibull" & g!="gompertz" & g!="gamma"
& g!="log-normal" & g!="chisq" & g!="f" & g!="burrxii" & g!="frechet"
& g!="lomax" & g!="log-logistic" & g!="lfr" & g!="chen")
{ stop ("Baseline distribution not implemented or misspelled. Please check the manual for guidelines.") }
if(g=="log-normal"){
den=function(par,x){meanlog=par[4]; sdlog=par[5]; loc=par[6]; dlnorm(x-loc,meanlog,sdlog)}
cum=function(par,x){meanlog=par[4]; sdlog=par[5]; loc=par[6]; plnorm(x-loc,meanlog,sdlog)}
starts<-c(1,1,1,log(median.mydata),sqrt(2*abs(log(mean.mydata/median.mydata))),(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){meanlog=par[4]; sdlog=par[5]; dlnorm(x,meanlog,sdlog)}
cum=function(par,x){meanlog=par[4]; sdlog=par[5]; plnorm(x,meanlog,sdlog)}
starts<-c(1,1,1,log(median.mydata),sqrt(2*abs(log(mean.mydata/median.mydata))))
}}
if(g=="chisq"){
den=function(par,x){df=par[4]; loc=par[5]; dchisq(x-loc,df)}
cum=function(par,x){df=par[4]; loc=par[5]; pchisq(x-loc,df)}
starts<-c(1,1,1,mean.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){df=par[4]; dchisq(x,df)}
cum=function(par,x){df=par[4]; pchisq(x,df)}
starts<-c(1,1,1,mean.mydata)
}}
if(g=="f"){
den=function(par,x){df1=par[4]; df2=par[5]; loc=par[6]; df(x-loc,df1,df2)}
cum=function(par,x){df1=par[4]; df2=par[5]; loc=par[6]; pf(x-loc,df1,df2)}
df.1<-ifelse(inv.mean>1,2*inv.mean/(inv.mean-1),1)
df.2<-ifelse(mean.mydata>1,2*mean.mydata/(mean.mydata-1),1)
starts<-c(1,1,1,df.1,df.2,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){df1=par[4]; df2=par[5]; df(x,df1,df2)}
cum=function(par,x){df1=par[4]; df2=par[5]; pf(x,df1,df2)}
starts<-c(1,1,1,df.1,df.2)
}}
if(g=="chen"){
den=function(par,x){a=par[4]; b=par[5]; loc=par[6]; a*b*(x-loc)^(a-1)*exp((x-loc)^a)*exp(-b*(exp((x-loc)^a)-1))}
cum=function(par,x){a=par[4]; b=par[5]; loc=par[6]; 1-exp(-b*(exp((x-loc)^a)-1))}
standard.mydata<-(y-mean.mydata)/std.mydata
s.hat<-ifelse(max(standard.mydata)<=1,-log(1-n/(n+.4))/(exp(1)-1),-
log(1-length(standard.mydata[standard.mydata<=1])/length(standard.mydata))/(exp(1)-1))
r.hat<-abs(log(log(1-log(1-0.5)/s.hat))/log(median.mydata))
starts<-c(1,1,1,r.hat,s.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[4]; b=par[5]; a*b*(x)^(a-1)*exp((x)^a)*exp(-b*(exp((x)^a)-1))}
cum=function(par,x){a=par[4]; b=par[5]; 1-exp(-b*(exp((x)^a)-1))}
starts<-c(1,1,1,r.hat,s.hat)
}}
if(g=="lfr"){
den=function(par,x){a=par[4]; b=par[5]; loc=par[6]; (a+b*(x-loc))*exp(-a*(x-loc)-(b*(x-loc)^2)/2)}
cum=function(par,x){a=par[4]; b=par[5]; loc=par[6]; 1-exp(-a*(x-loc)-(b*(x-loc)^2)/2)}
ii<-seq(1,(n-1))
yy<--log(1-(ii+0.3)/(n-1+0.4))
res<-summary(lm(yy ~ -1+y+I(y^2)))$coefficient
coeff<-c(abs(res[1,1]),abs(res[2,1]))
z1=NULL;
lfr.log<-function(p) {
z1<--log(sum(p[1]+p[2]*y))+p[1]*sum(y)+p[2]/2*sum(y^2)
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(abs(coeff[1]),2*abs(coeff[2])),lfr.log)$par)
starts<-c(1,1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[4]; b=par[5]; (a+b*(x))*exp(-a*(x)-(b*(x)^2)/2)}
cum=function(par,x){a=par[4]; b=par[5]; 1-exp(-a*(x)-(b*(x)^2)/2)}
starts<-c(1,1,1,p.hat[1],p.hat[2])
}}
if(g=="log-logistic"){
den=function(par,x){a=par[4]; b=par[5]; loc=par[6]; (a*b^(-a)*(x-loc)^(a-1))/((((x-loc)/b)^a +1)^2)}
cum=function(par,x){a=par[4]; b=par[5]; loc=par[6]; 1/(((x-loc)/b)^(-a)+1)}
starts<-c(1,1,1,log(0.75/(1-0.75))/log(qp3/median.mydata),median.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[4]; b=par[5]; (a*b^(-a)*(x)^(a-1))/((((x)/b)^a +1)^2)}
cum=function(par,x){a=par[4]; b=par[5]; 1/(((x)/b)^(-a)+1)}
starts<-c(1,1,1,log(0.75/(1-0.75))/log(qp3/median.mydata),median.mydata)
}}
if(g=="lomax"){
den=function(par,x){a=par[4]; b=par[5]; loc=par[6]; (a*b)/((1+a*(x-loc))^(b+1))}
cum=function(par,x){a=par[4]; b=par[5]; loc=par[6]; 1-(1+a*(x-loc))^(-b)}
a.hat<-1
b.hat<-(.5^(-1/a.hat)-1)/median.mydata
z1=NULL;
lomax.log<-function(p) {
z1<--n*log(p[1])-n*log(p[2])+(p[2]+1)*sum(log(1+p[1]*y))
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(1,b.hat),lomax.log, method)$par)
starts<-c(1,1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[4]; b=par[5]; (a*b)/((1+a*(x))^(b+1))}
cum=function(par,x){a=par[4]; b=par[5]; 1-(1+a*(x))^(-b)}
starts<-c(1,1,1,p.hat[1],p.hat[2])
}}
if(g=="gompertz"){
den=function(par,x){a=par[4]; b=par[5]; loc=par[6]; b*exp(a*(x-loc))*exp(-(exp(a*(x-loc))-1)*b/a)}
cum=function(par,x){a=par[4]; b=par[5]; loc=par[6]; 1-exp(-(exp(a*(x-loc))-1)*b/a)}
a.hat<-log(log(1-.75)/log(1-.5))/(qp3-median.mydata)
b.hat<--a.hat*log(0.5)/(exp(a.hat*median.mydata)-1)
z1=NULL;
gompertz.log<-function(p){
z1<--n*log(p[2])-p[1]*sum(y)+p[2]/p[1]*sum(exp(p[1]*y)-1)
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(a.hat,b.hat),gompertz.log, method)$par)
starts<-c(1,1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[4]; b=par[5]; b*exp(a*(x))*exp(-(exp(a*(x))-1)*b/a)}
cum=function(par,x){a=par[4]; b=par[5]; 1-exp(-(exp(a*(x))-1)*b/a)}
starts<-c(1,1,1,p.hat[1],p.hat[2])
}}
if(g=="exp"){
den=function(par,x){rate=par[4]; loc=par[5]; dexp(x-loc,rate)}
cum=function(par,x){rate=par[4]; loc=par[5]; pexp(x-loc,rate)}
starts<-c(1,1,1,1/mean.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){rate=par[4]; dexp(x,rate)}
cum=function(par,x){rate=par[4]; pexp(x,rate)}
starts<-c(1,1,1,1/mean.mydata)
}}
if(g=="rayleigh"){
den=function(par,x){a=par[4]; loc=par[5]; 2*(x-loc)/a*exp(-((x-loc)/a)^2)}
cum=function(par,x){a=par[4]; loc=par[5]; 1-exp(-((x-loc)/a)^2)}
starts<-c(1,1,1,log(2)/(median.mydata)^2,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[4]; 2*x/a*exp(-(x/a)^2)}
cum=function(par,x){a=par[4]; 1-exp(-(x/a)^2)}
starts<-c(1,1,1,log(2)/(median.mydata)^2)
}}
if(g=="burrxii"){
den=function(par,x){a=par[4]; d=par[5]; loc=par[6]; d*a*(1+(x-loc)^d)^(-a-1)*((x-loc)^(d-1))}
cum=function(par,x){a=par[4]; d=par[5]; loc=par[6]; 1-(1+(x-loc)^d)^(-a)}
z1=NULL;
burrxii.log<-function(p){
z1<--n*log(p[1])-n*log(p[2])-(p[2]-1)*sum(log(y))+(p[1]+1)*sum(log(1+y^p[2]))
}
p.hat<-suppressWarnings(optim(c(1,1),burrxii.log, method="BFGS")$par)
starts<-c(1,1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[4]; d=par[5]; d*a*(1+(x)^d)^(-a-1)*((x)^(d-1))}
cum=function(par,x){a=par[4]; d=par[5]; 1-(1+(x)^d)^(-a)}
starts<-c(1,1,1,p.hat[1],p.hat[2])
}}
if(g=="frechet"){
den=function(par,x){a=par[4]; b=par[5]; loc=par[6]; (a*exp(-((x-loc)/b)^(-a))*((x-loc)/b)^(-a-1))/(b)}
cum=function(par,x){a=par[4]; b=par[5]; loc=par[6]; exp(-((x-loc)/b)^(-a))}
cons<-cor(inv.mydata,rank(inv.mydata))*(sd(inv.mydata)/mean(inv.mydata))*sqrt((n+1)/(n-1))/sqrt(3)
a.hat<-abs(ifelse(cons<1,-log(2)/log(1-cons),-log(2)/log(1-.98)))
b.hat<-median(inv.mydata)/(log(2))^(1/a.hat)
starts<-c(1,1,1,a.hat,1/b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[4]; b=par[5]; (a*exp(-((x)/b)^(-a))*((x)/b)^(-a-1))/(b)}
cum=function(par,x){a=par[4]; b=par[5]; exp(-((x)/b)^(-a))}
starts<-c(1,1,1,a.hat,1/b.hat)
}}
if(g=="birnbaum-saunders"){
den=function(par,x){a=par[4]; b=par[5]; loc=par[6]; (sqrt((x-loc)/b)+sqrt(b/(x-loc)))/(2*a*(x-loc))*dnorm((sqrt((x-loc)/b)-sqrt(b/(x-loc)))/a)}
cum=function(par,x){a=par[4]; b=par[5]; loc=par[6]; pnorm((sqrt((x-loc)/b)-sqrt(b/(x-loc)))/a)}
r<-(mean(inv.mydata))^(-1)
b.hat<-sqrt(mean.mydata*r)
a.hat<-sqrt(2*(sqrt(mean.mydata/r)))
starts<-c(1,1,1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[4]; b=par[5]; (sqrt((x)/b)+sqrt(b/(x)))/(2*a*(x))*dnorm((sqrt((x)/b)-sqrt(b/(x)))/a)}
cum=function(par,x){a=par[4]; b=par[5]; pnorm((sqrt((x)/b)-sqrt(b/(x)))/a)}
starts<-c(1,1,1,a.hat,b.hat)
}}
if(g=="weibull"){
den=function(par,x){shape=par[4]; scale=par[5]; loc=par[6]; dweibull(x-loc,shape,scale)}
cum=function(par,x){shape=par[4]; scale=par[5]; loc=par[6]; pweibull(x-loc,shape,scale)}
cons<-cor(mydata,rank(mydata))*(std.mydata/mean.mydata)*sqrt((n+1)/(n-1))/sqrt(3)
a.hat<-abs(ifelse(cons<1,-log(2)/log(1-cons),-log(2)/log(1-.98)))
b.hat<-median.mydata/(log(2))^(1/a.hat)
starts<-c(1,1,1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){shape=par[4]; scale=par[5]; dweibull(x,shape,scale)}
cum=function(par,x){shape=par[4]; scale=par[5]; pweibull(x,shape,scale)}
starts<-c(1,1,1,a.hat,b.hat)
}}
if(g=="gamma"){
den=function(par,x){shape=par[4]; scale=par[5]; loc=par[6]; dgamma(x-loc,shape,scale)}
cum=function(par,x){shape=par[4]; scale=par[5]; loc=par[6]; pgamma(x-loc,shape,scale)}
a.hat<-uniroot(function(ss) trigamma(ss)-var(log(y)[is.finite(log(y))]),c(0,1000000))$root
b.hat<-mean.mydata/a.hat
starts<-c(1,1,1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){shape=par[4]; scale=par[5]; dgamma(x,shape,scale)}
cum=function(par,x){shape=par[4]; scale=par[5]; pgamma(x,shape,scale)}
starts<-c(1,1,1,a.hat,b.hat)
}}
pdf0<-function(par,x){
f0=den(par,x)
c0=cum(par,x)
a=par[1]
b=par[2]
d=par[3]
d*f0*beta(a,b)^(-1)*(1-c0)^(d*a-1)*(1-(1-c0)^d)^(b-1)}
cdf0<-function(par,x){
f0=den(par,x)
c0=cum(par,x)
a=par[1]
b=par[2]
d=par[3]
1-pbeta((1-c0)**d,shape1=a,shape2=b)}
mpsw<-function(par,x){
z=NULL
z=diff(c(0,cdf0(par,x),1))
for (j in 2:n){if((x[j]-x[j-1])==0){z[j]=pdf0(par,x[j])
}}
z[z<1e-16]=1e-16
-sum(log(z))
}
out<-suppressWarnings(optim(starts,fn=mpsw,x=sort.mydata,method=method))
n.p<-length(out$par)
if (location==TRUE){
out$par[1:(n.p-1)]<-abs(out$par[1:(n.p-1)])}
else{
out$par[1:n.p]<-abs(out$par[1:n.p])
}
gammam<-(n+1)*(log(n+1)-digamma(1))-1/2-1/(12*(n+1))
sigma2m<-(n+1)*(pi^2/6-1)-1/2-1/(6*(n+1))
c1<-gammam-sqrt(n*sigma2m/2)
c2<-sqrt(sigma2m/(2*n))
stat.chisquare<-(mpsw(out$par,sort.mydata)+n.p/2-c1)/c2
pvalue.chisquare<-pchisq(stat.chisquare,df=n,lower.tail=FALSE)
Moran<-out$value
log.likelihood=sum(log(pdf0(out$par,sort.mydata)))
u=cdf0(out$par,sort.mydata)
von<-c()
anderson<-c()
for(i in 1:n){
u[i]<-ifelse(u[i]==1,0.999999999,u[i])
von[i]=(u[i]-(2*i-1)/(2*n))^2
anderson[i]=(2*i-1)*log(u[i])+(2*n+1-2*i)*log(1-u[i])
}
anderson.stat=-n-mean(anderson)
von.stat=sum(von)+1/(12*n)
CAIC=-2*log.likelihood + 2*n.p + 2*(n.p*(n.p+1))/(n-n.p-1)
AIC=-2*log.likelihood + 2*n.p
BIC=-2*log.likelihood + n.p*log(n)
HQIC=-2*log.likelihood + 2*log(log(n))*n.p
ks.stat=suppressWarnings(ks.test(mydata, "cdf0", par=out$par))
aux2=cbind(AIC, CAIC, BIC, HQIC, von.stat, anderson.stat,log.likelihood,Moran)
colnames(aux2)=c("AIC","CAIC","BIC","HQIC","CM","AD", "log",
"Moran")
rownames(aux2)=c("")
aux3=cbind(ks.stat$statistic,ks.stat$p.value)
colnames(aux3)=c("statistic","p-value")
rownames(aux3)=c("")
aux4=cbind(stat.chisquare,qchisq(sig.level,df=n,lower.tail=FALSE),1-pchisq(stat.chisquare,df=n))
colnames(aux4)=c("statistic","chi-value","p-value")
rownames(aux4)=c("")
aux5=cbind(if(out$convergence==0){"Algorithm Converged"} else {"Algorithm Not Converged"})
list("MPS"=out$par,"Measures"=aux2,"KS"=aux3,"chi-square"=aux4,"Convergence Status"=aux5)
}
mpsbetag<-function(mydata, g, location=TRUE, method, sig.level){
min.mydata<-min(mydata)
sort.mydata<-sort(mydata)
n<-length(mydata)
y<-(sort.mydata[2:n]-min.mydata)
y<-(sort.mydata[2:n]-min.mydata)[is.finite(log(sort.mydata[2:n]-min.mydata))]
n<-length(mydata)
qp1<-y[floor(.25*n)]
qp3<-y[floor(.75*n)]
median.mydata<-y[floor(.5*n)]
mean.mydata<-mean(y)
inv.mydata<-1/y[is.finite(1/y)]
inv.mean<-mean(inv.mydata)
std.mydata<-sd(y)
if(g!="birnbaum-saunders" & g!="exp" & g!="rayleigh" & g!="weibull" & g!="gompertz" & g!="gamma"
& g!="log-normal" & g!="chisq" & g!="f" & g!="burrxii" & g!="frechet"
& g!="lomax" & g!="log-logistic" & g!="lfr" & g!="chen")
{ stop ("Baseline distribution not implemented or misspelled. Please check the manual for guidelines.") }
if(g=="log-normal"){
den=function(par,x){meanlog=par[3]; sdlog=par[4]; loc=par[5]; dlnorm(x-loc,meanlog,sdlog)}
cum=function(par,x){meanlog=par[3]; sdlog=par[4]; loc=par[5]; plnorm(x-loc,meanlog,sdlog)}
starts<-c(1,1,log(median.mydata),sqrt(2*abs(log(mean.mydata/median.mydata))),(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){meanlog=par[3]; sdlog=par[4]; dlnorm(x,meanlog,sdlog)}
cum=function(par,x){meanlog=par[3]; sdlog=par[4]; plnorm(x,meanlog,sdlog)}
starts<-c(1,1,log(median.mydata),sqrt(2*abs(log(mean.mydata/median.mydata))))
}}
if(g=="chisq"){
den=function(par,x){df=par[3]; loc=par[4]; dchisq(x-loc,df)}
cum=function(par,x){df=par[3]; loc=par[4]; pchisq(x-loc,df)}
starts<-c(1,1,mean.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){df=par[3]; dchisq(x,df)}
cum=function(par,x){df=par[3]; pchisq(x,df)}
starts<-c(1,1,mean.mydata)
}}
if(g=="f"){
den=function(par,x){df1=par[3]; df2=par[4]; loc=par[5]; df(x-loc,df1,df2)}
cum=function(par,x){df1=par[3]; df2=par[4]; loc=par[5]; pf(x-loc,df1,df2)}
df.1<-ifelse(inv.mean>1,2*inv.mean/(inv.mean-1),1)
df.2<-ifelse(mean.mydata>1,2*mean.mydata/(mean.mydata-1),1)
starts<-c(1,1,df.1,df.2,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){df1=par[3]; df2=par[4]; df(x,df1,df2)}
cum=function(par,x){df1=par[3]; df2=par[4]; pf(x,df1,df2)}
starts<-c(1,1,df.1,df.2)
}}
if(g=="chen"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; a*b*(x-loc)^(a-1)*exp((x-loc)^a)*exp(-b*(exp((x-loc)^a)-1))}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1-exp(-b*(exp((x-loc)^a)-1))}
standard.mydata<-(y-mean.mydata)/std.mydata
s.hat<-ifelse(max(standard.mydata)<=1,-log(1-n/(n+.4))/(exp(1)-1),-
log(1-length(standard.mydata[standard.mydata<=1])/length(standard.mydata))/(exp(1)-1))
r.hat<-abs(log(log(1-log(1-0.5)/s.hat))/log(median.mydata))
starts<-c(1,1,r.hat,s.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; a*b*(x)^(a-1)*exp((x)^a)*exp(-b*(exp((x)^a)-1))}
cum=function(par,x){a=par[3]; b=par[4]; 1-exp(-b*(exp((x)^a)-1))}
starts<-c(1,1,r.hat,s.hat)
}}
if(g=="lfr"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (a+b*(x-loc))*exp(-a*(x-loc)-(b*(x-loc)^2)/2)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1-exp(-a*(x-loc)-(b*(x-loc)^2)/2)}
ii<-seq(1,(n-1))
yy<--log(1-(ii+0.3)/(n-1+0.4))
res<-summary(lm(yy ~ -1+y+I(y^2)))$coefficient
coeff<-c(abs(res[1,1]),abs(res[2,1]))
z1=NULL;
lfr.log<-function(p) {
z1<--log(sum(p[1]+p[2]*y))+p[1]*sum(y)+p[2]/2*sum(y^2)
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(abs(coeff[1]),2*abs(coeff[2])),lfr.log)$par)
starts<-c(1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (a+b*(x))*exp(-a*(x)-(b*(x)^2)/2)}
cum=function(par,x){a=par[3]; b=par[4]; 1-exp(-a*(x)-(b*(x)^2)/2)}
starts<-c(1,1,p.hat[1],p.hat[2])
}}
if(g=="log-logistic"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (a*b^(-a)*(x-loc)^(a-1))/((((x-loc)/b)^a +1)^2)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1/(((x-loc)/b)^(-a)+1)}
starts<-c(1,1,log(0.75/(1-0.75))/log(qp3/median.mydata),median.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (a*b^(-a)*(x)^(a-1))/((((x)/b)^a +1)^2)}
cum=function(par,x){a=par[3]; b=par[4]; 1/(((x)/b)^(-a)+1)}
starts<-c(1,1,log(0.75/(1-0.75))/log(qp3/median.mydata),median.mydata)
}}
if(g=="lomax"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (a*b)/((1+a*(x-loc))^(b+1))}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1-(1+a*(x-loc))^(-b)}
a.hat<-1
b.hat<-(.5^(-1/a.hat)-1)/median.mydata
z1=NULL;
lomax.log<-function(p) {
z1<--n*log(p[1])-n*log(p[2])+(p[2]+1)*sum(log(1+p[1]*y))
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(1,b.hat),lomax.log, method)$par)
starts<-c(1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (a*b)/((1+a*(x))^(b+1))}
cum=function(par,x){a=par[3]; b=par[4]; 1-(1+a*(x))^(-b)}
starts<-c(1,1,p.hat[1],p.hat[2])
}}
if(g=="gompertz"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; b*exp(a*(x-loc))*exp(-(exp(a*(x-loc))-1)*b/a)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1-exp(-(exp(a*(x-loc))-1)*b/a)}
a.hat<-log(log(1-.75)/log(1-.5))/(qp3-median.mydata)
b.hat<--a.hat*log(0.5)/(exp(a.hat*median.mydata)-1)
z1=NULL;
gompertz.log<-function(p){
z1<--n*log(p[2])-p[1]*sum(y)+p[2]/p[1]*sum(exp(p[1]*y)-1)
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(a.hat,b.hat),gompertz.log, method)$par)
starts<-c(1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; b*exp(a*(x))*exp(-(exp(a*(x))-1)*b/a)}
cum=function(par,x){a=par[3]; b=par[4]; 1-exp(-(exp(a*(x))-1)*b/a)}
starts<-c(1,1,p.hat[1],p.hat[2])
}}
if(g=="exp"){
den=function(par,x){rate=par[3]; loc=par[4]; dexp(x-loc,rate)}
cum=function(par,x){rate=par[3]; loc=par[4]; pexp(x-loc,rate)}
starts<-c(1,1,1/mean.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){rate=par[3]; dexp(x,rate)}
cum=function(par,x){rate=par[3]; pexp(x,rate)}
starts<-c(1,1,1/mean.mydata)
}}
if(g=="rayleigh"){
den=function(par,x){a=par[3]; loc=par[4]; 2*(x-loc)/a*exp(-((x-loc)/a)^2)}
cum=function(par,x){a=par[3]; loc=par[4]; 1-exp(-((x-loc)/a)^2)}
starts<-c(1,1,log(2)/(median.mydata)^2,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; 2*x/a*exp(-(x/a)^2)}
cum=function(par,x){a=par[3]; 1-exp(-(x/a)^2)}
starts<-c(1,1,log(2)/(median.mydata)^2)
}}
if(g=="burrxii"){
den=function(par,x){a=par[3]; d=par[4]; loc=par[5]; d*a*(1+(x-loc)^d)^(-a-1)*((x-loc)^(d-1))}
cum=function(par,x){a=par[3]; d=par[4]; loc=par[5]; 1-(1+(x-loc)^d)^(-a)}
z1=NULL;
burrxii.log<-function(p){
z1<--n*log(p[1])-n*log(p[2])-(p[2]-1)*sum(log(y))+(p[1]+1)*sum(log(1+y^p[2]))
}
p.hat<-suppressWarnings(optim(c(1,1),burrxii.log, method="BFGS")$par)
starts<-c(1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; d=par[4]; d*a*(1+(x)^d)^(-a-1)*((x)^(d-1))}
cum=function(par,x){a=par[3]; d=par[4]; 1-(1+(x)^d)^(-a)}
starts<-c(1,1,p.hat[1],p.hat[2])
}}
if(g=="frechet"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (a*exp(-((x-loc)/b)^(-a))*((x-loc)/b)^(-a-1))/(b)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; exp(-((x-loc)/b)^(-a))}
cons<-cor(inv.mydata,rank(inv.mydata))*(sd(inv.mydata)/mean(inv.mydata))*sqrt((n+1)/(n-1))/sqrt(3)
a.hat<-abs(ifelse(cons<1,-log(2)/log(1-cons),-log(2)/log(1-.98)))
b.hat<-median(inv.mydata)/(log(2))^(1/a.hat)
starts<-c(1,1,a.hat,1/b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (a*exp(-((x)/b)^(-a))*((x)/b)^(-a-1))/(b)}
cum=function(par,x){a=par[3]; b=par[4]; exp(-((x)/b)^(-a))}
starts<-c(1,1,a.hat,1/b.hat)
}}
if(g=="birnbaum-saunders"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (sqrt((x-loc)/b)+sqrt(b/(x-loc)))/(2*a*(x-loc))*dnorm((sqrt((x-loc)/b)-sqrt(b/(x-loc)))/a)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; pnorm((sqrt((x-loc)/b)-sqrt(b/(x-loc)))/a)}
r<-(mean(inv.mydata))^(-1)
b.hat<-sqrt(mean.mydata*r)
a.hat<-sqrt(2*(sqrt(mean.mydata/r)))
starts<-c(1,1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (sqrt((x)/b)+sqrt(b/(x)))/(2*a*(x))*dnorm((sqrt((x)/b)-sqrt(b/(x)))/a)}
cum=function(par,x){a=par[3]; b=par[4]; pnorm((sqrt((x)/b)-sqrt(b/(x)))/a)}
starts<-c(1,1,a.hat,b.hat)
}}
if(g=="weibull"){
den=function(par,x){shape=par[3]; scale=par[4]; loc=par[5]; dweibull(x-loc,shape,scale)}
cum=function(par,x){shape=par[3]; scale=par[4]; loc=par[5]; pweibull(x-loc,shape,scale)}
cons<-cor(mydata,rank(mydata))*(std.mydata/mean.mydata)*sqrt((n+1)/(n-1))/sqrt(3)
a.hat<-abs(ifelse(cons<1,-log(2)/log(1-cons),-log(2)/log(1-.98)))
b.hat<-median.mydata/(log(2))^(1/a.hat)
starts<-c(1,1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){shape=par[3]; scale=par[4]; dweibull(x,shape,scale)}
cum=function(par,x){shape=par[3]; scale=par[4]; pweibull(x,shape,scale)}
starts<-c(1,1,a.hat,b.hat)
}}
if(g=="gamma"){
den=function(par,x){shape=par[3]; scale=par[4]; loc=par[5]; dgamma(x-loc,shape,scale)}
cum=function(par,x){shape=par[3]; scale=par[4]; loc=par[5]; pgamma(x-loc,shape,scale)}
a.hat<-uniroot(function(ss) trigamma(ss)-var(log(y)[is.finite(log(y))]),c(0,1000000))$root
b.hat<-mean.mydata/a.hat
starts<-c(1,1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){shape=par[3]; scale=par[4]; dgamma(x,shape,scale)}
cum=function(par,x){shape=par[3]; scale=par[4]; pgamma(x,shape,scale)}
starts<-c(1,1,a.hat,b.hat)
}}
pdf0<-function(par,x){
f0=den(par,x)
c0=cum(par,x)
a=par[1]
b=par[2]
f0*beta(a,b)^(-1)*c0^(a-1)*(1-c0)^(b-1)}
cdf0<-function(par,x){
c0=cum(par,x)
a=par[1]
b=par[2]
pbeta(c0,shape1=a,shape2=b)}
mpsw<-function(par,x){
z=NULL
z=diff(c(0,cdf0(par,x),1))
for (j in 2:n){if((x[j]-x[j-1])==0){z[j]=pdf0(par,x[j])
}}
z[z<1e-16]=1e-16
-sum(log(z))
}
out<-suppressWarnings(optim(starts,fn=mpsw,x=sort.mydata,method=method))
n.p<-length(out$par)
gammam<-(n+1)*(log(n+1)-digamma(1))-1/2-1/(12*(n+1))
sigma2m<-(n+1)*(pi^2/6-1)-1/2-1/(6*(n+1))
c1<-gammam-sqrt(n*sigma2m/2)
c2<-sqrt(sigma2m/(2*n))
stat.chisquare<-(mpsw(out$par,sort.mydata)+n.p/2-c1)/c2
pvalue.chisquare<-pchisq(stat.chisquare,df=n,lower.tail=FALSE)
Moran<-out$value
log.likelihood=sum(log(pdf0(out$par,sort.mydata)))
u=cdf0(out$par,sort.mydata)
von<-c()
anderson<-c()
for(i in 1:n){
u[i]<-ifelse(u[i]==1,0.999999999,u[i])
von[i]=(u[i]-(2*i-1)/(2*n))^2
anderson[i]=(2*i-1)*log(u[i])+(2*n+1-2*i)*log(1-u[i])
}
anderson.stat=-n-mean(anderson)
von.stat=sum(von)+1/(12*n)
CAIC=-2*log.likelihood + 2*n.p + 2*(n.p*(n.p+1))/(n-n.p-1)
AIC=-2*log.likelihood + 2*n.p
BIC=-2*log.likelihood + n.p*log(n)
HQIC=-2*log.likelihood + 2*log(log(n))*n.p
ks.stat=suppressWarnings(ks.test(mydata, "cdf0", par=out$par))
aux2=cbind(AIC, CAIC, BIC, HQIC, von.stat, anderson.stat,log.likelihood,Moran)
colnames(aux2)=c("AIC","CAIC","BIC","HQIC","CM","AD", "log", "Moran")
rownames(aux2)=c("")
aux3=cbind(ks.stat$statistic,ks.stat$p.value)
colnames(aux3)=c("statistic","p-value")
rownames(aux3)=c("")
aux4=cbind(stat.chisquare,qchisq(sig.level,df=n,lower.tail=FALSE),1-pchisq(stat.chisquare,df=n))
colnames(aux4)=c("statistic","chi-value","p-value")
rownames(aux4)=c("")
aux5=cbind(if(out$convergence==0){"Algorithm Converged"} else {"Algorithm Not Converged"})
list("MPS"=out$par,"Measures"=aux2,"KS"=aux3,"chi-square"=aux4,"Convergence Status"=aux5)
}
mpsexpexppg<-function(mydata, g, location=TRUE, method , sig.level){
min.mydata<-min(mydata)
sort.mydata<-sort(mydata)
n<-length(mydata)
y<-(sort.mydata[2:n]-min.mydata)
y<-(sort.mydata[2:n]-min.mydata)[is.finite(log(sort.mydata[2:n]-min.mydata))]
n<-length(mydata)
qp1<-y[floor(.25*n)]
qp3<-y[floor(.75*n)]
median.mydata<-y[floor(.5*n)]
mean.mydata<-mean(y)
inv.mydata<-1/y[is.finite(1/y)]
inv.mean<-mean(inv.mydata)
std.mydata<-sd(y)
if(g!="birnbaum-saunders" & g!="exp" & g!="rayleigh" & g!="weibull" & g!="gompertz" & g!="gamma"
& g!="log-normal" & g!="chisq" & g!="f" & g!="burrxii" & g!="frechet"
& g!="lomax" & g!="log-logistic" & g!="lfr" & g!="chen")
{ stop ("Baseline distribution not implemented or misspelled. Please check the manual for guidelines.") }
if(g=="log-normal"){
den=function(par,x){meanlog=par[3]; sdlog=par[4]; loc=par[5]; dlnorm(x-loc,meanlog,sdlog)}
cum=function(par,x){meanlog=par[3]; sdlog=par[4]; loc=par[5]; plnorm(x-loc,meanlog,sdlog)}
starts<-c(1,1,log(median.mydata),sqrt(2*abs(log(mean.mydata/median.mydata))),(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){meanlog=par[3]; sdlog=par[4]; dlnorm(x,meanlog,sdlog)}
cum=function(par,x){meanlog=par[3]; sdlog=par[4]; plnorm(x,meanlog,sdlog)}
starts<-c(1,1,log(median.mydata),sqrt(2*abs(log(mean.mydata/median.mydata))))
}}
if(g=="chisq"){
den=function(par,x){df=par[3]; loc=par[4]; dchisq(x-loc,df)}
cum=function(par,x){df=par[3]; loc=par[4]; pchisq(x-loc,df)}
starts<-c(1,1,mean.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){df=par[3]; dchisq(x,df)}
cum=function(par,x){df=par[3]; pchisq(x,df)}
starts<-c(1,1,mean.mydata)
}}
if(g=="f"){
den=function(par,x){df1=par[3]; df2=par[4]; loc=par[5]; df(x-loc,df1,df2)}
cum=function(par,x){df1=par[3]; df2=par[4]; loc=par[5]; pf(x-loc,df1,df2)}
df.1<-ifelse(inv.mean>1,2*inv.mean/(inv.mean-1),1)
df.2<-ifelse(mean.mydata>1,2*mean.mydata/(mean.mydata-1),1)
starts<-c(1,1,df.1,df.2,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){df1=par[3]; df2=par[4]; df(x,df1,df2)}
cum=function(par,x){df1=par[3]; df2=par[4]; pf(x,df1,df2)}
starts<-c(1,1,df.1,df.2)
}}
if(g=="chen"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; a*b*(x-loc)^(a-1)*exp((x-loc)^a)*exp(-b*(exp((x-loc)^a)-1))}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1-exp(-b*(exp((x-loc)^a)-1))}
standard.mydata<-(y-mean.mydata)/std.mydata
s.hat<-ifelse(max(standard.mydata)<=1,-log(1-n/(n+.4))/(exp(1)-1),-
log(1-length(standard.mydata[standard.mydata<=1])/length(standard.mydata))/(exp(1)-1))
r.hat<-abs(log(log(1-log(1-0.5)/s.hat))/log(median.mydata))
starts<-c(1,1,r.hat,s.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; a*b*(x)^(a-1)*exp((x)^a)*exp(-b*(exp((x)^a)-1))}
cum=function(par,x){a=par[3]; b=par[4]; 1-exp(-b*(exp((x)^a)-1))}
starts<-c(1,1,r.hat,s.hat)
}}
if(g=="lfr"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (a+b*(x-loc))*exp(-a*(x-loc)-(b*(x-loc)^2)/2)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1-exp(-a*(x-loc)-(b*(x-loc)^2)/2)}
ii<-seq(1,(n-1))
yy<--log(1-(ii+0.3)/(n-1+0.4))
res<-summary(lm(yy ~ -1+y+I(y^2)))$coefficient
coeff<-c(abs(res[1,1]),abs(res[2,1]))
z1=NULL;
lfr.log<-function(p) {
z1<--log(sum(p[1]+p[2]*y))+p[1]*sum(y)+p[2]/2*sum(y^2)
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(abs(coeff[1]),2*abs(coeff[2])),lfr.log)$par)
starts<-c(1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (a+b*(x))*exp(-a*(x)-(b*(x)^2)/2)}
cum=function(par,x){a=par[3]; b=par[4]; 1-exp(-a*(x)-(b*(x)^2)/2)}
starts<-c(1,1,p.hat[1],p.hat[2])
}}
if(g=="log-logistic"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (a*b^(-a)*(x-loc)^(a-1))/((((x-loc)/b)^a +1)^2)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1/(((x-loc)/b)^(-a)+1)}
starts<-c(1,1,log(0.75/(1-0.75))/log(qp3/median.mydata),median.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (a*b^(-a)*(x)^(a-1))/((((x)/b)^a +1)^2)}
cum=function(par,x){a=par[3]; b=par[4]; 1/(((x)/b)^(-a)+1)}
starts<-c(1,1,log(0.75/(1-0.75))/log(qp3/median.mydata),median.mydata)
}}
if(g=="lomax"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (a*b)/((1+a*(x-loc))^(b+1))}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1-(1+a*(x-loc))^(-b)}
a.hat<-1
b.hat<-(.5^(-1/a.hat)-1)/median.mydata
z1=NULL;
lomax.log<-function(p) {
z1<--n*log(p[1])-n*log(p[2])+(p[2]+1)*sum(log(1+p[1]*y))
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(1,b.hat),lomax.log, method)$par)
starts<-c(1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (a*b)/((1+a*(x))^(b+1))}
cum=function(par,x){a=par[3]; b=par[4]; 1-(1+a*(x))^(-b)}
starts<-c(1,1,p.hat[1],p.hat[2])
}}
if(g=="gompertz"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; b*exp(a*(x-loc))*exp(-(exp(a*(x-loc))-1)*b/a)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1-exp(-(exp(a*(x-loc))-1)*b/a)}
a.hat<-log(log(1-.75)/log(1-.5))/(qp3-median.mydata)
b.hat<--a.hat*log(0.5)/(exp(a.hat*median.mydata)-1)
z1=NULL;
gompertz.log<-function(p){
z1<--n*log(p[2])-p[1]*sum(y)+p[2]/p[1]*sum(exp(p[1]*y)-1)
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(a.hat,b.hat),gompertz.log, method)$par)
starts<-c(1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; b*exp(a*(x))*exp(-(exp(a*(x))-1)*b/a)}
cum=function(par,x){a=par[3]; b=par[4]; 1-exp(-(exp(a*(x))-1)*b/a)}
starts<-c(1,1,p.hat[1],p.hat[2])
}}
if(g=="exp"){
den=function(par,x){rate=par[3]; loc=par[4]; dexp(x-loc,rate)}
cum=function(par,x){rate=par[3]; loc=par[4]; pexp(x-loc,rate)}
starts<-c(1,1,1/mean.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){rate=par[3]; dexp(x,rate)}
cum=function(par,x){rate=par[3]; pexp(x,rate)}
starts<-c(1,1,1/mean.mydata)
}}
if(g=="rayleigh"){
den=function(par,x){a=par[3]; loc=par[4]; 2*(x-loc)/a*exp(-((x-loc)/a)^2)}
cum=function(par,x){a=par[3]; loc=par[4]; 1-exp(-((x-loc)/a)^2)}
starts<-c(1,1,log(2)/(median.mydata)^2,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; 2*x/a*exp(-(x/a)^2)}
cum=function(par,x){a=par[3]; 1-exp(-(x/a)^2)}
starts<-c(1,1,log(2)/(median.mydata)^2)
}}
if(g=="burrxii"){
den=function(par,x){a=par[3]; d=par[4]; loc=par[5]; d*a*(1+(x-loc)^d)^(-a-1)*((x-loc)^(d-1))}
cum=function(par,x){a=par[3]; d=par[4]; loc=par[5]; 1-(1+(x-loc)^d)^(-a)}
z1=NULL;
burrxii.log<-function(p){
z1<--n*log(p[1])-n*log(p[2])-(p[2]-1)*sum(log(y))+(p[1]+1)*sum(log(1+y^p[2]))
}
p.hat<-suppressWarnings(optim(c(1,1),burrxii.log, method="BFGS")$par)
starts<-c(1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; d=par[4]; d*a*(1+(x)^d)^(-a-1)*((x)^(d-1))}
cum=function(par,x){a=par[3]; d=par[4]; 1-(1+(x)^d)^(-a)}
starts<-c(1,1,p.hat[1],p.hat[2])
}}
if(g=="frechet"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (a*exp(-((x-loc)/b)^(-a))*((x-loc)/b)^(-a-1))/(b)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; exp(-((x-loc)/b)^(-a))}
cons<-cor(inv.mydata,rank(inv.mydata))*(sd(inv.mydata)/mean(inv.mydata))*sqrt((n+1)/(n-1))/sqrt(3)
a.hat<-abs(ifelse(cons<1,-log(2)/log(1-cons),-log(2)/log(1-.98)))
b.hat<-median(inv.mydata)/(log(2))^(1/a.hat)
starts<-c(1,1,a.hat,1/b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (a*exp(-((x)/b)^(-a))*((x)/b)^(-a-1))/(b)}
cum=function(par,x){a=par[3]; b=par[4]; exp(-((x)/b)^(-a))}
starts<-c(1,1,a.hat,1/b.hat)
}}
if(g=="birnbaum-saunders"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (sqrt((x-loc)/b)+sqrt(b/(x-loc)))/(2*a*(x-loc))*dnorm((sqrt((x-loc)/b)-sqrt(b/(x-loc)))/a)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; pnorm((sqrt((x-loc)/b)-sqrt(b/(x-loc)))/a)}
r<-(mean(inv.mydata))^(-1)
b.hat<-sqrt(mean.mydata*r)
a.hat<-sqrt(2*(sqrt(mean.mydata/r)))
starts<-c(1,1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (sqrt((x)/b)+sqrt(b/(x)))/(2*a*(x))*dnorm((sqrt((x)/b)-sqrt(b/(x)))/a)}
cum=function(par,x){a=par[3]; b=par[4]; pnorm((sqrt((x)/b)-sqrt(b/(x)))/a)}
starts<-c(1,1,a.hat,b.hat)
}}
if(g=="weibull"){
den=function(par,x){shape=par[3]; scale=par[4]; loc=par[5]; dweibull(x-loc,shape,scale)}
cum=function(par,x){shape=par[3]; scale=par[4]; loc=par[5]; pweibull(x-loc,shape,scale)}
cons<-cor(mydata,rank(mydata))*(std.mydata/mean.mydata)*sqrt((n+1)/(n-1))/sqrt(3)
a.hat<-abs(ifelse(cons<1,-log(2)/log(1-cons),-log(2)/log(1-.98)))
b.hat<-median.mydata/(log(2))^(1/a.hat)
starts<-c(1,1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){shape=par[3]; scale=par[4]; dweibull(x,shape,scale)}
cum=function(par,x){shape=par[3]; scale=par[4]; pweibull(x,shape,scale)}
starts<-c(1,1,a.hat,b.hat)
}}
if(g=="gamma"){
den=function(par,x){shape=par[3]; scale=par[4]; loc=par[5]; dgamma(x-loc,shape,scale)}
cum=function(par,x){shape=par[3]; scale=par[4]; loc=par[5]; pgamma(x-loc,shape,scale)}
a.hat<-uniroot(function(ss) trigamma(ss)-var(log(y)[is.finite(log(y))]),c(0,1000000))$root
b.hat<-mean.mydata/a.hat
starts<-c(1,1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){shape=par[3]; scale=par[4]; dgamma(x,shape,scale)}
cum=function(par,x){shape=par[3]; scale=par[4]; pgamma(x,shape,scale)}
starts<-c(1,1,a.hat,b.hat)
}}
pdf0<-function(par,x){
f0=den(par,x)
c0=cum(par,x)
a=par[1]
b=par[2]
a*b*(1-exp(-b))^(-1)*f0*c0^(a-1)*exp(-b*c0^a)}
cdf0<-function(par,x){
f0=den(par,x)
c0=cum(par,x)
a=par[1]
b=par[2]
(1-exp(-b))^(-1)*(1-exp(-b*c0^a))}
mpsw<-function(par,x){
z=NULL
z=diff(c(0,cdf0(par,x),1))
for (j in 2:n){if((x[j]-x[j-1])==0){z[j]=pdf0(par,x[j])
}}
z[z<1e-16]=1e-16
-sum(log(z))
}
out<-suppressWarnings(optim(starts,fn=mpsw,x=sort.mydata,method=method))
n.p<-length(out$par)
gammam<-(n+1)*(log(n+1)-digamma(1))-1/2-1/(12*(n+1))
sigma2m<-(n+1)*(pi^2/6-1)-1/2-1/(6*(n+1))
c1<-gammam-sqrt(n*sigma2m/2)
c2<-sqrt(sigma2m/(2*n))
stat.chisquare<-(mpsw(out$par,sort.mydata)+n.p/2-c1)/c2
pvalue.chisquare<-pchisq(stat.chisquare,df=n,lower.tail=FALSE)
Moran<-out$value
log.likelihood=sum(log(pdf0(out$par,sort.mydata)))
u=cdf0(out$par,sort.mydata)
von<-c()
anderson<-c()
for(i in 1:n){
u[i]<-ifelse(u[i]==1,0.999999999,u[i])
von[i]=(u[i]-(2*i-1)/(2*n))^2
anderson[i]=(2*i-1)*log(u[i])+(2*n+1-2*i)*log(1-u[i])
}
anderson.stat=-n-mean(anderson)
von.stat=sum(von)+1/(12*n)
CAIC=-2*log.likelihood + 2*n.p + 2*(n.p*(n.p+1))/(n-n.p-1)
AIC=-2*log.likelihood + 2*n.p
BIC=-2*log.likelihood + n.p*log(n)
HQIC=-2*log.likelihood + 2*log(log(n))*n.p
ks.stat=suppressWarnings(ks.test(mydata, "cdf0", par=out$par))
aux2=cbind(AIC, CAIC, BIC, HQIC, von.stat, anderson.stat,log.likelihood,Moran)
colnames(aux2)=c("AIC","CAIC","BIC","HQIC","CM","AD", "log",
"Moran")
rownames(aux2)=c("")
aux3=cbind(ks.stat$statistic,ks.stat$p.value)
colnames(aux3)=c("statistic","p-value")
rownames(aux3)=c("")
aux4=cbind(stat.chisquare,qchisq(sig.level,df=n,lower.tail=FALSE),1-pchisq(stat.chisquare,df=n))
colnames(aux4)=c("statistic","chi-value","p-value")
rownames(aux4)=c("")
aux5=cbind(if(out$convergence==0){"Algorithm Converged"} else {"Algorithm Not Converged"})
list("MPS"=out$par,"Measures"=aux2,"KS"=aux3,"chi-square"=aux4,"Convergence Status"=aux5)
}
mpsexpgg<-function(mydata, g, location=TRUE, method , sig.level){
min.mydata<-min(mydata)
sort.mydata<-sort(mydata)
n<-length(mydata)
y<-(sort.mydata[2:n]-min.mydata)
y<-(sort.mydata[2:n]-min.mydata)[is.finite(log(sort.mydata[2:n]-min.mydata))]
n<-length(mydata)
qp1<-y[floor(.25*n)]
qp3<-y[floor(.75*n)]
median.mydata<-y[floor(.5*n)]
mean.mydata<-mean(y)
inv.mydata<-1/y[is.finite(1/y)]
inv.mean<-mean(inv.mydata)
std.mydata<-sd(y)
if(g!="birnbaum-saunders" & g!="exp" & g!="rayleigh" & g!="weibull" & g!="gompertz" & g!="gamma"
& g!="log-normal" & g!="chisq" & g!="f" & g!="burrxii" & g!="frechet"
& g!="lomax" & g!="log-logistic" & g!="lfr" & g!="chen")
{ stop ("Baseline distribution not implemented or misspelled. Please check the manual for guidelines.") }
if(g=="log-normal"){
den=function(par,x){meanlog=par[3]; sdlog=par[4]; loc=par[5]; dlnorm(x-loc,meanlog,sdlog)}
cum=function(par,x){meanlog=par[3]; sdlog=par[4]; loc=par[5]; plnorm(x-loc,meanlog,sdlog)}
starts<-c(1,1,log(median.mydata),sqrt(2*abs(log(mean.mydata/median.mydata))),(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){meanlog=par[3]; sdlog=par[4]; dlnorm(x,meanlog,sdlog)}
cum=function(par,x){meanlog=par[3]; sdlog=par[4]; plnorm(x,meanlog,sdlog)}
starts<-c(1,1,log(median.mydata),sqrt(2*abs(log(mean.mydata/median.mydata))))
}}
if(g=="chisq"){
den=function(par,x){df=par[3]; loc=par[4]; dchisq(x-loc,df)}
cum=function(par,x){df=par[3]; loc=par[4]; pchisq(x-loc,df)}
starts<-c(1,1,mean.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){df=par[3]; dchisq(x,df)}
cum=function(par,x){df=par[3]; pchisq(x,df)}
starts<-c(1,1,mean.mydata)
}}
if(g=="f"){
den=function(par,x){df1=par[3]; df2=par[4]; loc=par[5]; df(x-loc,df1,df2)}
cum=function(par,x){df1=par[3]; df2=par[4]; loc=par[5]; pf(x-loc,df1,df2)}
df.1<-ifelse(inv.mean>1,2*inv.mean/(inv.mean-1),1)
df.2<-ifelse(mean.mydata>1,2*mean.mydata/(mean.mydata-1),1)
starts<-c(1,1,df.1,df.2,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){df1=par[3]; df2=par[4]; df(x,df1,df2)}
cum=function(par,x){df1=par[3]; df2=par[4]; pf(x,df1,df2)}
starts<-c(1,1,df.1,df.2)
}}
if(g=="chen"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; a*b*(x-loc)^(a-1)*exp((x-loc)^a)*exp(-b*(exp((x-loc)^a)-1))}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1-exp(-b*(exp((x-loc)^a)-1))}
standard.mydata<-(y-mean.mydata)/std.mydata
s.hat<-ifelse(max(standard.mydata)<=1,-log(1-n/(n+.4))/(exp(1)-1),-
log(1-length(standard.mydata[standard.mydata<=1])/length(standard.mydata))/(exp(1)-1))
r.hat<-abs(log(log(1-log(1-0.5)/s.hat))/log(median.mydata))
starts<-c(1,1,r.hat,s.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; a*b*(x)^(a-1)*exp((x)^a)*exp(-b*(exp((x)^a)-1))}
cum=function(par,x){a=par[3]; b=par[4]; 1-exp(-b*(exp((x)^a)-1))}
starts<-c(1,1,r.hat,s.hat)
}}
if(g=="lfr"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (a+b*(x-loc))*exp(-a*(x-loc)-(b*(x-loc)^2)/2)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1-exp(-a*(x-loc)-(b*(x-loc)^2)/2)}
ii<-seq(1,(n-1))
yy<--log(1-(ii+0.3)/(n-1+0.4))
res<-summary(lm(yy ~ -1+y+I(y^2)))$coefficient
coeff<-c(abs(res[1,1]),abs(res[2,1]))
z1=NULL;
lfr.log<-function(p) {
z1<--log(sum(p[1]+p[2]*y))+p[1]*sum(y)+p[2]/2*sum(y^2)
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(abs(coeff[1]),2*abs(coeff[2])),lfr.log)$par)
starts<-c(1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (a+b*(x))*exp(-a*(x)-(b*(x)^2)/2)}
cum=function(par,x){a=par[3]; b=par[4]; 1-exp(-a*(x)-(b*(x)^2)/2)}
starts<-c(1,1,p.hat[1],p.hat[2])
}}
if(g=="log-logistic"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (a*b^(-a)*(x-loc)^(a-1))/((((x-loc)/b)^a +1)^2)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1/(((x-loc)/b)^(-a)+1)}
starts<-c(1,1,log(0.75/(1-0.75))/log(qp3/median.mydata),median.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (a*b^(-a)*(x)^(a-1))/((((x)/b)^a +1)^2)}
cum=function(par,x){a=par[3]; b=par[4]; 1/(((x)/b)^(-a)+1)}
starts<-c(1,1,log(0.75/(1-0.75))/log(qp3/median.mydata),median.mydata)
}}
if(g=="lomax"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (a*b)/((1+a*(x-loc))^(b+1))}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1-(1+a*(x-loc))^(-b)}
a.hat<-1
b.hat<-(.5^(-1/a.hat)-1)/median.mydata
z1=NULL;
lomax.log<-function(p) {
z1<--n*log(p[1])-n*log(p[2])+(p[2]+1)*sum(log(1+p[1]*y))
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(1,b.hat),lomax.log, method)$par)
starts<-c(1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (a*b)/((1+a*(x))^(b+1))}
cum=function(par,x){a=par[3]; b=par[4]; 1-(1+a*(x))^(-b)}
starts<-c(1,1,p.hat[1],p.hat[2])
}}
if(g=="gompertz"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; b*exp(a*(x-loc))*exp(-(exp(a*(x-loc))-1)*b/a)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; 1-exp(-(exp(a*(x-loc))-1)*b/a)}
a.hat<-log(log(1-.75)/log(1-.5))/(qp3-median.mydata)
b.hat<--a.hat*log(0.5)/(exp(a.hat*median.mydata)-1)
z1=NULL;
gompertz.log<-function(p){
z1<--n*log(p[2])-p[1]*sum(y)+p[2]/p[1]*sum(exp(p[1]*y)-1)
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(a.hat,b.hat),gompertz.log, method)$par)
starts<-c(1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; b*exp(a*(x))*exp(-(exp(a*(x))-1)*b/a)}
cum=function(par,x){a=par[3]; b=par[4]; 1-exp(-(exp(a*(x))-1)*b/a)}
starts<-c(1,1,p.hat[1],p.hat[2])
}}
if(g=="exp"){
den=function(par,x){rate=par[3]; loc=par[4]; dexp(x-loc,rate)}
cum=function(par,x){rate=par[3]; loc=par[4]; pexp(x-loc,rate)}
starts<-c(1,1,1/mean.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){rate=par[3]; dexp(x,rate)}
cum=function(par,x){rate=par[3]; pexp(x,rate)}
starts<-c(1,1,1/mean.mydata)
}}
if(g=="rayleigh"){
den=function(par,x){a=par[3]; loc=par[4]; 2*(x-loc)/a*exp(-((x-loc)/a)^2)}
cum=function(par,x){a=par[3]; loc=par[4]; 1-exp(-((x-loc)/a)^2)}
starts<-c(1,1,log(2)/(median.mydata)^2,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; 2*x/a*exp(-(x/a)^2)}
cum=function(par,x){a=par[3]; 1-exp(-(x/a)^2)}
starts<-c(1,1,log(2)/(median.mydata)^2)
}}
if(g=="burrxii"){
den=function(par,x){a=par[3]; d=par[4]; loc=par[5]; d*a*(1+(x-loc)^d)^(-a-1)*((x-loc)^(d-1))}
cum=function(par,x){a=par[3]; d=par[4]; loc=par[5]; 1-(1+(x-loc)^d)^(-a)}
z1=NULL;
burrxii.log<-function(p){
z1<--n*log(p[1])-n*log(p[2])-(p[2]-1)*sum(log(y))+(p[1]+1)*sum(log(1+y^p[2]))
}
p.hat<-suppressWarnings(optim(c(1,1),burrxii.log, method="BFGS")$par)
starts<-c(1,1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; d=par[4]; d*a*(1+(x)^d)^(-a-1)*((x)^(d-1))}
cum=function(par,x){a=par[3]; d=par[4]; 1-(1+(x)^d)^(-a)}
starts<-c(1,1,p.hat[1],p.hat[2])
}}
if(g=="frechet"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (a*exp(-((x-loc)/b)^(-a))*((x-loc)/b)^(-a-1))/(b)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; exp(-((x-loc)/b)^(-a))}
cons<-cor(inv.mydata,rank(inv.mydata))*(sd(inv.mydata)/mean(inv.mydata))*sqrt((n+1)/(n-1))/sqrt(3)
a.hat<-abs(ifelse(cons<1,-log(2)/log(1-cons),-log(2)/log(1-.98)))
b.hat<-median(inv.mydata)/(log(2))^(1/a.hat)
starts<-c(1,1,a.hat,1/b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (a*exp(-((x)/b)^(-a))*((x)/b)^(-a-1))/(b)}
cum=function(par,x){a=par[3]; b=par[4]; exp(-((x)/b)^(-a))}
starts<-c(1,1,a.hat,1/b.hat)
}}
if(g=="birnbaum-saunders"){
den=function(par,x){a=par[3]; b=par[4]; loc=par[5]; (sqrt((x-loc)/b)+sqrt(b/(x-loc)))/(2*a*(x-loc))*dnorm((sqrt((x-loc)/b)-sqrt(b/(x-loc)))/a)}
cum=function(par,x){a=par[3]; b=par[4]; loc=par[5]; pnorm((sqrt((x-loc)/b)-sqrt(b/(x-loc)))/a)}
r<-(mean(inv.mydata))^(-1)
b.hat<-sqrt(mean.mydata*r)
a.hat<-sqrt(2*(sqrt(mean.mydata/r)))
starts<-c(1,1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[3]; b=par[4]; (sqrt((x)/b)+sqrt(b/(x)))/(2*a*(x))*dnorm((sqrt((x)/b)-sqrt(b/(x)))/a)}
cum=function(par,x){a=par[3]; b=par[4]; pnorm((sqrt((x)/b)-sqrt(b/(x)))/a)}
starts<-c(1,1,a.hat,b.hat)
}}
if(g=="weibull"){
den=function(par,x){shape=par[3]; scale=par[4]; loc=par[5]; dweibull(x-loc,shape,scale)}
cum=function(par,x){shape=par[3]; scale=par[4]; loc=par[5]; pweibull(x-loc,shape,scale)}
cons<-cor(mydata,rank(mydata))*(std.mydata/mean.mydata)*sqrt((n+1)/(n-1))/sqrt(3)
a.hat<-abs(ifelse(cons<1,-log(2)/log(1-cons),-log(2)/log(1-.98)))
b.hat<-median.mydata/(log(2))^(1/a.hat)
starts<-c(1,1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){shape=par[3]; scale=par[4]; dweibull(x,shape,scale)}
cum=function(par,x){shape=par[3]; scale=par[4]; pweibull(x,shape,scale)}
starts<-c(1,1,a.hat,b.hat)
}}
if(g=="gamma"){
den=function(par,x){shape=par[3]; scale=par[4]; loc=par[5]; dgamma(x-loc,shape,scale)}
cum=function(par,x){shape=par[3]; scale=par[4]; loc=par[5]; pgamma(x-loc,shape,scale)}
a.hat<-uniroot(function(ss) trigamma(ss)-var(log(y)[is.finite(log(y))]),c(0,1000000))$root
b.hat<-mean.mydata/a.hat
starts<-c(1,1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){shape=par[3]; scale=par[4]; dgamma(x,shape,scale)}
cum=function(par,x){shape=par[3]; scale=par[4]; pgamma(x,shape,scale)}
starts<-c(1,1,a.hat,b.hat)
}}
pdf0<-function(par,x){
f0=den(par,x)
c0=cum(par,x)
r=par[1]
u=par[2]
r*u*f0*((1-c0)^(r-1))*(1-(1-c0)^r)^(u-1)
}
cdf0<-function(par,x){
f0=den(par,x)
c0=cum(par,x)
r=par[1]
u=par[2]
(1-(1-c0)^r)^u
}
mpsw<-function(par,x){
z=NULL
z=diff(c(0,cdf0(par,x),1))
for (j in 2:n){if((x[j]-x[j-1])==0){z[j]=pdf0(par,x[j])
}}
z[z<1e-16]=1e-16
-sum(log(z))
}
out<-suppressWarnings(optim(starts,fn=mpsw,x=sort.mydata,method=method))
n.p<-length(out$par)
gammam<-(n+1)*(log(n+1)-digamma(1))-1/2-1/(12*(n+1))
sigma2m<-(n+1)*(pi^2/6-1)-1/2-1/(6*(n+1))
c1<-gammam-sqrt(n*sigma2m/2)
c2<-sqrt(sigma2m/(2*n))
stat.chisquare<-(mpsw(out$par,sort.mydata)+n.p/2-c1)/c2
pvalue.chisquare<-pchisq(stat.chisquare,df=n,lower.tail=FALSE)
Moran<-out$value
log.likelihood=sum(log(pdf0(out$par,sort.mydata)))
u=cdf0(out$par,sort.mydata)
von<-c()
anderson<-c()
for(i in 1:n){
u[i]<-ifelse(u[i]==1,0.999999999,u[i])
von[i]=(u[i]-(2*i-1)/(2*n))^2
anderson[i]=(2*i-1)*log(u[i])+(2*n+1-2*i)*log(1-u[i])
}
anderson.stat=-n-mean(anderson)
von.stat=sum(von)+1/(12*n)
CAIC=-2*log.likelihood + 2*n.p + 2*(n.p*(n.p+1))/(n-n.p-1)
AIC=-2*log.likelihood + 2*n.p
BIC=-2*log.likelihood + n.p*log(n)
HQIC=-2*log.likelihood + 2*log(log(n))*n.p
ks.stat=suppressWarnings(ks.test(mydata, "cdf0", par=out$par))
aux2=cbind(AIC, CAIC, BIC, HQIC, von.stat, anderson.stat,log.likelihood,Moran)
colnames(aux2)=c("AIC","CAIC","BIC","HQIC","CM","AD", "log",
"Moran")
rownames(aux2)=c("")
aux3=cbind(ks.stat$statistic,ks.stat$p.value)
colnames(aux3)=c("statistic","p-value")
rownames(aux3)=c("")
aux4=cbind(stat.chisquare,qchisq(sig.level,df=n,lower.tail=FALSE),1-pchisq(stat.chisquare,df=n))
colnames(aux4)=c("statistic","chi-value","p-value")
rownames(aux4)=c("")
aux5=cbind(if(out$convergence==0){"Algorithm Converged"} else {"Algorithm Not Converged"})
list("MPS"=out$par,"Measures"=aux2,"KS"=aux3,"chi-square"=aux4,"Convergence Status"=aux5)
}
mpsexpg<-function(mydata, g, location=TRUE, method , sig.level){
min.mydata<-min(mydata)
sort.mydata<-sort(mydata)
n<-length(mydata)
y<-(sort.mydata[2:n]-min.mydata)
y<-(sort.mydata[2:n]-min.mydata)[is.finite(log(sort.mydata[2:n]-min.mydata))]
n<-length(mydata)
qp1<-y[floor(.25*n)]
qp3<-y[floor(.75*n)]
median.mydata<-y[floor(.5*n)]
mean.mydata<-mean(y)
inv.mydata<-1/y[is.finite(1/y)]
inv.mean<-mean(inv.mydata)
std.mydata<-sd(y)
if(g!="birnbaum-saunders" & g!="exp" & g!="rayleigh" & g!="weibull" & g!="gompertz" & g!="gamma"
& g!="log-normal" & g!="chisq" & g!="f" & g!="burrxii" & g!="frechet"
& g!="lomax" & g!="log-logistic" & g!="lfr" & g!="chen")
{ stop ("Baseline distribution not implemented or misspelled. Please check the manual for guidelines.") }
if(g=="log-normal"){
den=function(par,x){meanlog=par[2]; sdlog=par[3]; loc=par[4]; dlnorm(x-loc,meanlog,sdlog)}
cum=function(par,x){meanlog=par[2]; sdlog=par[3]; loc=par[4]; plnorm(x-loc,meanlog,sdlog)}
starts<-c(1,log(median.mydata),sqrt(2*abs(log(mean.mydata/median.mydata))),(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){meanlog=par[2]; sdlog=par[3]; dlnorm(x,meanlog,sdlog)}
cum=function(par,x){meanlog=par[2]; sdlog=par[3]; plnorm(x,meanlog,sdlog)}
starts<-c(1,log(median.mydata),sqrt(2*abs(log(mean.mydata/median.mydata))))
}}
if(g=="chisq"){
den=function(par,x){df=par[2]; loc=par[3]; dchisq(x-loc,df)}
cum=function(par,x){df=par[2]; loc=par[3]; pchisq(x-loc,df)}
starts<-c(1,mean.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){df=par[2]; dchisq(x,df)}
cum=function(par,x){df=par[2]; pchisq(x,df)}
starts<-c(1,mean.mydata)
}}
if(g=="f"){
den=function(par,x){df1=par[2]; df2=par[3]; loc=par[4]; df(x-loc,df1,df2)}
cum=function(par,x){df1=par[2]; df2=par[3]; loc=par[4]; pf(x-loc,df1,df2)}
df.1<-ifelse(inv.mean>1,2*inv.mean/(inv.mean-1),1)
df.2<-ifelse(mean.mydata>1,2*mean.mydata/(mean.mydata-1),1)
starts<-c(1,df.1,df.2,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){df1=par[2]; df2=par[3]; df(x,df1,df2)}
cum=function(par,x){df1=par[2]; df2=par[3]; pf(x,df1,df2)}
starts<-c(1,df.1,df.2)
}}
if(g=="chen"){
den=function(par,x){a=par[2]; b=par[3]; loc=par[4]; a*b*(x-loc)^(a-1)*exp((x-loc)^a)*exp(-b*(exp((x-loc)^a)-1))}
cum=function(par,x){a=par[2]; b=par[3]; loc=par[4]; 1-exp(-b*(exp((x-loc)^a)-1))}
standard.mydata<-(y-mean.mydata)/std.mydata
s.hat<-ifelse(max(standard.mydata)<=1,-log(1-n/(n+.4))/(exp(1)-1),-
log(1-length(standard.mydata[standard.mydata<=1])/length(standard.mydata))/(exp(1)-1))
r.hat<-abs(log(log(1-log(1-0.5)/s.hat))/log(median.mydata))
starts<-c(1,r.hat,s.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[2]; b=par[3]; a*b*(x)^(a-1)*exp((x)^a)*exp(-b*(exp((x)^a)-1))}
cum=function(par,x){a=par[2]; b=par[3]; 1-exp(-b*(exp((x)^a)-1))}
starts<-c(1,r.hat,s.hat)
}}
if(g=="lfr"){
den=function(par,x){a=par[2]; b=par[3]; loc=par[4]; (a+b*(x-loc))*exp(-a*(x-loc)-(b*(x-loc)^2)/2)}
cum=function(par,x){a=par[2]; b=par[3]; loc=par[4]; 1-exp(-a*(x-loc)-(b*(x-loc)^2)/2)}
ii<-seq(1,(n-1))
yy<--log(1-(ii+0.3)/(n-1+0.4))
res<-summary(lm(yy ~ -1+y+I(y^2)))$coefficient
coeff<-c(abs(res[1,1]),abs(res[2,1]))
z1=NULL;
lfr.log<-function(p) {
z1<--log(sum(p[1]+p[2]*y))+p[1]*sum(y)+p[2]/2*sum(y^2)
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(abs(coeff[1]),2*abs(coeff[2])),lfr.log)$par)
starts<-c(1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[2]; b=par[3]; (a+b*(x))*exp(-a*(x)-(b*(x)^2)/2)}
cum=function(par,x){a=par[2]; b=par[3]; 1-exp(-a*(x)-(b*(x)^2)/2)}
starts<-c(1,p.hat[1],p.hat[2])
}}
if(g=="log-logistic"){
den=function(par,x){a=par[2]; b=par[3]; loc=par[4]; (a*b^(-a)*(x-loc)^(a-1))/((((x-loc)/b)^a +1)^2)}
cum=function(par,x){a=par[2]; b=par[3]; loc=par[4]; 1/(((x-loc)/b)^(-a)+1)}
starts<-c(1,log(0.75/(1-0.75))/log(qp3/median.mydata),median.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[2]; b=par[3]; (a*b^(-a)*(x)^(a-1))/((((x)/b)^a +1)^2)}
cum=function(par,x){a=par[2]; b=par[3]; 1/(((x)/b)^(-a)+1)}
starts<-c(1,log(0.75/(1-0.75))/log(qp3/median.mydata),median.mydata)
}}
if(g=="lomax"){
den=function(par,x){a=par[2]; b=par[3]; loc=par[4]; (a*b)/((1+a*(x-loc))^(b+1))}
cum=function(par,x){a=par[2]; b=par[3]; loc=par[4]; 1-(1+a*(x-loc))^(-b)}
a.hat<-1
b.hat<-(.5^(-1/a.hat)-1)/median.mydata
z1=NULL;
lomax.log<-function(p) {
z1<--n*log(p[1])-n*log(p[2])+(p[2]+1)*sum(log(1+p[1]*y))
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(1,b.hat),lomax.log, method)$par)
starts<-c(1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[2]; b=par[3]; (a*b)/((1+a*(x))^(b+1))}
cum=function(par,x){a=par[2]; b=par[3]; 1-(1+a*(x))^(-b)}
starts<-c(1,p.hat[1],p.hat[2])
}}
if(g=="gompertz"){
den=function(par,x){a=par[2]; b=par[3]; loc=par[4]; b*exp(a*(x-loc))*exp(-(exp(a*(x-loc))-1)*b/a)}
cum=function(par,x){a=par[2]; b=par[3]; loc=par[4]; 1-exp(-(exp(a*(x-loc))-1)*b/a)}
a.hat<-log(log(1-.75)/log(1-.5))/(qp3-median.mydata)
b.hat<--a.hat*log(0.5)/(exp(a.hat*median.mydata)-1)
z1=NULL;
gompertz.log<-function(p){
z1<--n*log(p[2])-p[1]*sum(y)+p[2]/p[1]*sum(exp(p[1]*y)-1)
z1[z1<1e-16]<-1e-16
}
p.hat<-suppressWarnings(optim(c(a.hat,b.hat),gompertz.log, method)$par)
starts<-c(1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[2]; b=par[3]; b*exp(a*(x))*exp(-(exp(a*(x))-1)*b/a)}
cum=function(par,x){a=par[2]; b=par[3]; 1-exp(-(exp(a*(x))-1)*b/a)}
starts<-c(1,p.hat[1],p.hat[2])
}}
if(g=="exp"){
den=function(par,x){rate=par[2]; loc=par[3]; dexp(x-loc,rate)}
cum=function(par,x){rate=par[2]; loc=par[3]; pexp(x-loc,rate)}
starts<-c(1,1/mean.mydata,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){rate=par[2]; dexp(x,rate)}
cum=function(par,x){rate=par[2]; pexp(x,rate)}
starts<-c(1,1/mean.mydata)
}}
if(g=="rayleigh"){
den=function(par,x){a=par[2]; loc=par[3]; 2*(x-loc)/a*exp(-((x-loc)/a)^2)}
cum=function(par,x){a=par[2]; loc=par[3]; 1-exp(-((x-loc)/a)^2)}
starts<-c(1,log(2)/(median.mydata)^2,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[2]; 2*x/a*exp(-(x/a)^2)}
cum=function(par,x){a=par[2]; 1-exp(-(x/a)^2)}
starts<-c(1,log(2)/(median.mydata)^2)
}}
if(g=="burrxii"){
den=function(par,x){a=par[2]; d=par[3]; loc=par[4]; d*a*(1+(x-loc)^d)^(-a-1)*((x-loc)^(d-1))}
cum=function(par,x){a=par[2]; d=par[3]; loc=par[4]; 1-(1+(x-loc)^d)^(-a)}
z1=NULL;
burrxii.log<-function(p){
z1<--n*log(p[1])-n*log(p[2])-(p[2]-1)*sum(log(y))+(p[1]+1)*sum(log(1+y^p[2]))
}
p.hat<-suppressWarnings(optim(c(1,1),burrxii.log, method="BFGS")$par)
starts<-c(1,p.hat[1],p.hat[2],(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[2]; d=par[3]; d*a*(1+(x)^d)^(-a-1)*((x)^(d-1))}
cum=function(par,x){a=par[2]; d=par[3]; 1-(1+(x)^d)^(-a)}
starts<-c(1,p.hat[1],p.hat[2])
}}
if(g=="frechet"){
den=function(par,x){a=par[2]; b=par[3]; loc=par[4]; (a*exp(-((x-loc)/b)^(-a))*((x-loc)/b)^(-a-1))/(b)}
cum=function(par,x){a=par[2]; b=par[3]; loc=par[4]; exp(-((x-loc)/b)^(-a))}
cons<-cor(inv.mydata,rank(inv.mydata))*(sd(inv.mydata)/mean(inv.mydata))*sqrt((n+1)/(n-1))/sqrt(3)
a.hat<-abs(ifelse(cons<1,-log(2)/log(1-cons),-log(2)/log(1-.98)))
b.hat<-median(inv.mydata)/(log(2))^(1/a.hat)
starts<-c(1,a.hat,1/b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[2]; b=par[3]; (a*exp(-((x)/b)^(-a))*((x)/b)^(-a-1))/(b)}
cum=function(par,x){a=par[2]; b=par[3]; exp(-((x)/b)^(-a))}
starts<-c(1,a.hat,1/b.hat)
}}
if(g=="birnbaum-saunders"){
den=function(par,x){a=par[2]; b=par[3]; loc=par[4]; (sqrt((x-loc)/b)+sqrt(b/(x-loc)))/(2*a*(x-loc))*dnorm((sqrt((x-loc)/b)-sqrt(b/(x-loc)))/a)}
cum=function(par,x){a=par[2]; b=par[3]; loc=par[4]; pnorm((sqrt((x-loc)/b)-sqrt(b/(x-loc)))/a)}
r<-(mean(inv.mydata))^(-1)
b.hat<-sqrt(mean.mydata*r)
a.hat<-sqrt(2*(sqrt(mean.mydata/r)))
starts<-c(1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){a=par[2]; b=par[3]; (sqrt((x)/b)+sqrt(b/(x)))/(2*a*(x))*dnorm((sqrt((x)/b)-sqrt(b/(x)))/a)}
cum=function(par,x){a=par[2]; b=par[3]; pnorm((sqrt((x)/b)-sqrt(b/(x)))/a)}
starts<-c(1,a.hat,b.hat)
}}
if(g=="weibull"){
den=function(par,x){shape=par[2]; scale=par[3]; loc=par[4]; dweibull(x-loc,shape,scale)}
cum=function(par,x){shape=par[2]; scale=par[3]; loc=par[4]; pweibull(x-loc,shape,scale)}
cons<-cor(mydata,rank(mydata))*(std.mydata/mean.mydata)*sqrt((n+1)/(n-1))/sqrt(3)
a.hat<-abs(ifelse(cons<1,-log(2)/log(1-cons),-log(2)/log(1-.98)))
b.hat<-median.mydata/(log(2))^(1/a.hat)
starts<-c(1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){shape=par[2]; scale=par[3]; dweibull(x,shape,scale)}
cum=function(par,x){shape=par[2]; scale=par[3]; pweibull(x,shape,scale)}
starts<-c(1,a.hat,b.hat)
}}
if(g=="gamma"){
den=function(par,x){shape=par[2]; scale=par[3]; loc=par[4]; dgamma(x-loc,shape,scale)}
cum=function(par,x){shape=par[2]; scale=par[3]; loc=par[4]; pgamma(x-loc,shape,scale)}
a.hat<-uniroot(function(ss) trigamma(ss)-var(log(y)[is.finite(log(y))]),c(0,1000000))$root
b.hat<-mean.mydata/a.hat
starts<-c(1,a.hat,b.hat,(min.mydata-1/length(y)))
if (location==FALSE){
den=function(par,x){shape=par[2]; scale=par[3]; dgamma(x,shape,scale)}
cum=function(par,x){shape=par[2]; scale=par[3]; pgamma(x,shape,scale)}
starts<-c(1,a.hat,b.hat)
}}
pdf0<-function(par,x){
f0=den(par,x)
c0=cum(par,x)
a=par[1]
a*f0*c0^(a-1)
}
cdf0<-function(par,x){
f0=den(par,x)
c0=cum(par,x)
a=par[1]
c0^a
}
mpsw<-function(par,x){
z=NULL
z=diff(c(0,cdf0(par,x),1))
for (j in 2:n){if((x[j]-x[j-1])==0){z[j]=pdf0(par,x[j])
}}
z[z<1e-16]=1e-16
-sum(log(z))
}
out<-suppressWarnings(optim(starts,fn=mpsw,x=sort.mydata,method=method))
n.p<-length(out$pa