GOF=function(x, tipo, a){
#' @title Goodness-of-fit-test
#' @description This function estimates the parameters and performs goodness-of-fit-test for the common density functions and probability mass functions.
#' @param x Data.
#' @param tipo Data Type: discreto (discrete) or continuo (continuous).
#' @param a Number of simulations for the p-value, 4000 for default.
#' @return Results of parameter estimates and the Goodness of fit tests
#' @export GOF
#' @author Israel Hernandez Luna
#' @references Hernandez, Israel (2018). SoluciĆ³n en R para las pruebas de bondad de ajuste. UNAM.
#' @examples
#' x=c(rep(0,109),rep(1,65),rep(2,22),rep(3,3),4) #muestra
#' GOF(x,"discreto")
#' GOF(x,"discreto",10000)
#'
if (missing(x) || mode(x) != "numeric")
stop("'x' deberia ser un vector numerico no vacio")
if (any(!is.finite(x)))
stop("'x' contiene valores perdidos o infinitos")
if (!is.character(tipo) || missing(tipo))
stop("Falta el tipo de dato (continuo o discreto)")
if (missing(a))
a=4000
if(a < 1000)
stop("El numero de simulaciones para el p-valor debe de ser mayor a 1000")
a=floor(a)
fun=tolower(tipo)
if (fun == "discreto" || fun =="d")
{
if (any(x < 0))
stop("existen valores menores a 0, las distribuciones Binomial,Poisson, Bin Negativa, Geometrica tienen valores positivos")
else
{
#Estadistica W^2 Cramer-von Mises
W=function(x,pj)
{
n=length(x)
a=sort(unique(x))
k=length(a)
ej=pj*n
oj=numeric(k) #cuantos hay de cada xi
for( i in 1:k)
{
oj[i]=sum(x==a[i])
}
zj=cumsum(oj-ej)
w=sum((zj^2)*pj)/n
return(w)
}
#Estadistica U^2 Watson
U=function(x,pj)
{
n=length(x)
a=sort(unique(x))
k=length(a)
ej=pj*n
oj=numeric(k)
for( i in 1:k)
{
oj[i]=sum(x==a[i])
}
zj=cumsum(oj-ej)
z=sum(zj*pj)
u=sum(((zj-z)^2)*pj)/n
return(u)
}
#Estadistica A^2 Anderson-Darling
A=function(x,pj)
{
n=length(x)
a=sort(unique(x))
k=length(a)
ej=pj*n
hj=cumsum(pj)
oj=numeric(k)
for( i in 1:k)
{
oj[i]=sum(x==a[i])
}
zj=cumsum(oj-ej)
aa=(zj^2)*pj/(hj*(1-hj))
aa[aa == 1/0] = 0
Am=sum(aa)/n
if(Am == 1/0) return(1) else return(Am)
}
parametros=rep(0,7)
names(parametros)=c("n.Uniforme" ,"size.Binom", "p.Binom" ,"l.Poisson" , "p.Geom", "size.BN", "p.BN")
Wm=rep(0,5)
Um=rep(0,5)
Am=rep(0,5)
names(Wm)=c("Uniforme", "Binomial" ,"Poisson", "Geometrica", "Binomial.Neg")
Wp.valor=rep(0,5)
Up.valor=rep(0,5)
Ap.valor=rep(0,5)
n=length(x)
m=mean(x)
v=var(x)
#Distribucion Uniforme
dis.unif=function(x,N) c(rep(1,length(x)))/N
N.u=max(x)
parametros[1]=N.u
prob0=dis.unif(sort(unique(x)),N.u)
Wm[1]=W(x,prob0)
Um[1]=U(x,prob0)
Am[1]=A(x,prob0)
xx=unique(x)
#p.valor
W0=numeric(length(a))
U0=numeric(length(a))
A0=numeric(length(a))
for(i in 1:a)
{
ru=sample(xx,n,replace=TRUE)
W0[i]=W(ru,dis.unif(sort(unique(ru)),N.u))
U0[i]=U(ru,dis.unif(sort(unique(ru)),N.u))
A0[i]=A(ru,dis.unif(sort(unique(ru)),N.u))
}
Wp.valor[1]=sum(W0>Wm[1])/a
Up.valor[1]=sum(U0>Um[1])/a
Ap.valor[1]=sum(A0>Am[1])/a
#Distribucion Binomial
nb=0
nbb=floor((max(x)^1.95)*(v^.95)/((m^.95)*(max(x)-m)^.95)+.5)
if (max(x) > nbb) nb=max(x) else nb=nbb
pb=m/nb
parametros[2]=nb
parametros[3]=pb
prob1=dbinom(sort(unique(x)),nb,pb)
Wm[2]=W(x,prob1)
Um[2]=U(x,prob1)
Am[2]=A(x,prob1)
#p.valor
W0=numeric(length(a))
U0=numeric(length(a))
A0=numeric(length(a))
for(i in 1:a)
{
rb=rbinom(n,nb,pb)
W0[i]=W(rb,dbinom(sort(unique(rb)),nb,pb))
U0[i]=W(rb,dbinom(sort(unique(rb)),nb,pb))
A0[i]=W(rb,dbinom(sort(unique(rb)),nb,pb))
}
Wp.valor[2]=sum(W0>Wm[2])/a
Up.valor[2]=sum(U0>Um[2])/a
Ap.valor[2]=sum(A0>Am[2])/a
#Distribucion Poisson
l=m
prob2=dpois(sort(unique(x)),l)
Wm[3]=W(x,prob2)
Um[3]=U(x,prob2)
Am[3]=A(x,prob2)
parametros[4]=l
#p.valor
W0=numeric(length(a))
U0=numeric(length(a))
A0=numeric(length(a))
for(i in 1:a)
{
rp=rpois(n,l)
W0[i]=W(rp,dpois(sort(unique(rp)),l))
U0[i]=U(rp,dpois(sort(unique(rp)),l))
A0[i]=A(rp,dpois(sort(unique(rp)),l))
}
Wp.valor[3]=sum(W0>Wm[3])/a
Up.valor[3]=sum(U0>Um[3])/a
Ap.valor[3]=sum(A0>Am[3])/a
#Distribucion Geometrica
pg= 1/(1+mean(x))
prob3=dgeom(sort(unique(x)),pg)
Wm[4]=W(x,prob3)
Um[4]=U(x,prob3)
Am[4]=A(x,prob3)
parametros[5]=pg
#p.valor
W0=numeric(length(a))
U0=numeric(length(a))
A0=numeric(length(a))
for(i in 1:a)
{
rg=rgeom(n,pg)
W0[i]=W(rg,dgeom(sort(unique(rg)),pg))
U0[i]=U(rg,dgeom(sort(unique(rg)),pg))
A0[i]=A(rg,dgeom(sort(unique(rg)),pg))
}
Wp.valor[4]=sum(W0>Wm[4])/a
Up.valor[4]=sum(U0>Um[4])/a
Ap.valor[4]=sum(A0>Am[4])/a
#Distribucion Binomial Negativa
if( v > m)
{
size=m^2/(v-m)
pbn=size/(size+m)
prob4=dnbinom(sort(unique(x)),size,pbn)
Wm[5]=W(x,prob4)
Um[5]=U(x,prob4)
Am[5]=A(x,prob4)
parametros[6]=size
parametros[7]=pbn
#p.valor
W0=numeric(length(a))
U0=numeric(length(a))
A0=numeric(length(a))
for(i in 1:a)
{
rbn=rnbinom(n,size,pbn)
W0[i]=W(rbn,dnbinom(sort(unique(rbn)),size,pbn))
U0[i]=U(rbn,dnbinom(sort(unique(rbn)),size,pbn))
A0[i]=A(rbn,dnbinom(sort(unique(rbn)),size,pbn))
}
Wp.valor[5]=sum(W0>Wm[5])/a
Up.valor[5]=sum(U0>Um[5])/a
Ap.valor[5]=sum(A0>Am[5])/a
}
else
{
Wm[4]=0
Um[4]=0
Am[4]=0
Wp.valor[4]=0
Up.valor[4]=0
Ap.valor[4]=0
}
}
return(list(rbind(Wm,Wp.valor,Um,Up.valor,Am,Ap.valor),parametros))
}
if(fun == "continuo" || fun == "c")
{
#Estadistica Dn de Kolmogorov-Smirnov
Dn=function(pr)
{
pr=sort(pr)
n=length(pr)
a=numeric(n)
aa=numeric(n)
b=seq(1:n)
for( i in 1:n)
{
a[i]=(b[i]/n)-pr[i]
aa[i]=pr[i]-((b[i]-1)/n)
}
d=max(a)
dd=max(aa)
return(max(c(d,dd)))
}
#Estadistica Vn de Kuiper
Vn=function(pr)
{
pr=sort(pr)
n=length(pr)
a=numeric(n)
aa=numeric(n)
b=seq(1:n)
for( i in 1:n)
{
a[i]=(b[i]/n)-pr[i]
aa[i]=pr[i]-((b[i]-1)/n)
}
return(max(a)+max(aa))
}
#Estadistica Wn^2 de Cramer-von Mises
Wn=function(pr)
{
pr=sort(pr)
n=length(pr)
a=numeric(n)
b=seq(1:n)
for( i in 1:n)
{
a[i]=(pr[i]-(2*b[i]-1)/(2*n))^2
}
return(sum(a)+1/(12*n))
}
#Estadistica An^2 Anderson-Darling
An=function(pr)
{
pr=sort(pr)
n=length(pr)
a=numeric(n)
b=seq(1:n)
for(i in 1:n){
a[i]=(1/n)*((2*b[i]-1)*log(pr[i])+(2*n+1-2*b[i])*log(1-pr[i]))
}
a=a[a!=-Inf]
n=length(a)
return(-n-sum(a))
}
#Estadistica Un^2 de Watson
Un=function(pr)
{
pr=sort(pr)
n=length(pr)
a=numeric(n)
b=seq(1:n)
for( i in 1:n)
{
a[i]=(pr[i]-(2*b[i]-1)/(2*n))^2
}
c=(sum(a)+1/(12*n))
return(c-n*(mean(pr)-.5)^2)
}
parametros=rep(0,17)
names(parametros)=c("min.U", "max.U", "l.Exp", "mu.Normal", "sd.Normal" , "mu.LogN", "sd.LogN", "a.Gamma" , "b.Gamma", "a.Cauchy" , "b.Cauchy" , "a.Beta" , "b.Beta" , "a.Weibull" , "b.Weibull" , "k.Ji^2" , "k.T" )
Dm=rep(0,10)
Vm=rep(0,10)
Wm=rep(0,10)
Am=rep(0,10)
Um=rep(0,10)
names(Dm)=c("Uniforme", "Exponencial", "Normal" , "Log-Normal" , "Gamma" , "Cauchy" , "Beta", "Weibull" , "Ji^2" , "T" )
Dp.valor=rep(0,10)
Vp.valor=rep(0,10)
Wp.valor=rep(0,10)
Ap.valor=rep(0,10)
Up.valor=rep(0,10)
n=length(x)
m=mean(x)
v=var(x)
#Distribucion Uniforme
mx=ceiling(max(x))
mn=floor(min(x))
prob1=punif(sort(x),mn,mx)
Dm[1]=Dn(prob1)
Vm[1]=Vn(prob1)
Wm[1]=Wn(prob1)
Am[1]=An(prob1)
Um[1]=Un(prob1)
parametros[1]=mn
parametros[2]=mx
#p.valor
D0=numeric(length(a))
V0=numeric(length(a))
W0=numeric(length(a))
A0=numeric(length(a))
U0=numeric(length(a))
for(i in 1:a)
{
ru=runif(n,mn,mx)
p=punif(sort(ru),mn,mx)
D0[i]=Dn(p)
V0[i]=Vn(p)
W0[i]=Wn(p)
A0[i]=An(p)
U0[i]=Un(p)
}
Dp.valor[1]=sum(D0>Dm[1])/a
Vp.valor[1]=sum(V0>Vm[1])/a
Wp.valor[1]=sum(W0>Wm[1])/a
Ap.valor[1]=sum(A0>Am[1])/a
Up.valor[1]=sum(U0>Um[1])/a
#Distribucion Exponencial
if (any(x <= 0))
{
Dm[2]=0
Vm[2]=0
Wm[2]=0
Am[2]=0
Um[2]=0
parametros[3]=0
Dp.valor[2]=0
Vp.valor[2]=0
Wp.valor[2]=0
Ap.valor[2]=0
Up.valor[2]=0
}
else
{
lambda.e= 1/m
prob2=pexp(sort(x),lambda.e)
Dm[2]=Dn(prob2)
Vm[2]=Vn(prob2)
Wm[2]=Wn(prob2)
Am[2]=An(prob2)
Um[2]=Un(prob2)
parametros[3]=lambda.e
#p.valor
D0=numeric(length(a))
V0=numeric(length(a))
W0=numeric(length(a))
A0=numeric(length(a))
U0=numeric(length(a))
for(i in 1:a)
{
rex=rexp(n,lambda.e)
p=pexp(sort(rex),lambda.e)
D0[i]=Dn(p)
V0[i]=Vn(p)
W0[i]=Wn(p)
A0[i]=An(p)
U0[i]=Un(p)
}
Dp.valor[2]=sum(D0>Dm[2])/a
Vp.valor[2]=sum(V0>Vm[2])/a
Wp.valor[2]=sum(W0>Wm[2])/a
Ap.valor[2]=sum(A0>Am[2])/a
Up.valor[2]=sum(U0>Um[2])/a
}
#Distribucion Normal
xx=(x-m)^2
sd0=sqrt((1/n)*sum(xx))
prob3=pnorm(sort(x),m,sd0)
Dm[3]=Dn(prob3)
Vm[3]=Vn(prob3)
Wm[3]=Wn(prob3)
Am[3]=An(prob3)
Um[3]=Un(prob3)
parametros[4]=m
parametros[5]=sd0
#p.valor
D0=numeric(length(a))
V0=numeric(length(a))
W0=numeric(length(a))
A0=numeric(length(a))
U0=numeric(length(a))
for(i in 1:a)
{
rn=rnorm(n,m,sd0)
p=pnorm(sort(rn),m,sd0)
D0[i]=Dn(p)
V0[i]=Vn(p)
W0[i]=Wn(p)
A0[i]=An(p)
U0[i]=Un(p)
}
Dp.valor[3]=sum(D0>Dm[3])/a
Vp.valor[3]=sum(V0>Vm[3])/a
Wp.valor[3]=sum(W0>Wm[3])/a
Ap.valor[3]=sum(A0>Am[3])/a
Up.valor[3]=sum(U0>Um[3])/a
#Distribucion Lognormal
if (any(x <= 0))
{
Dm[4]=0
Vm[4]=0
Wm[4]=0
Am[4]=0
Um[4]=0
parametros[6]=0
parametros[7]=0
Dp.valor[4]=0
Vp.valor[4]=0
Wp.valor[4]=0
Ap.valor[4]=0
Up.valor[4]=0
}
else
{
mu.l=sum(log(x))/n
lxx=(log(x)-mu.l)^2
sd.l=sqrt((1/n)*sum(lxx))
prob4=plnorm(sort(x),mu.l,sd.l)
Dm[4]=Dn(prob4)
Vm[4]=Vn(prob4)
Wm[4]=Wn(prob4)
Am[4]=An(prob4)
Um[4]=Un(prob4)
parametros[6]=mu.l
parametros[7]=sd.l
#p.valor
D0=numeric(length(a))
V0=numeric(length(a))
W0=numeric(length(a))
A0=numeric(length(a))
U0=numeric(length(a))
for(i in 1:a)
{
rl=rlnorm(n,mu.l,sd.l)
p=plnorm(sort(rl),mu.l,sd.l)
D0[i]=Dn(p)
V0[i]=Vn(p)
W0[i]=Wn(p)
A0[i]=An(p)
U0[i]=Un(p)
}
Dp.valor[4]=sum(D0>Dm[4])/a
Vp.valor[4]=sum(V0>Vm[4])/a
Wp.valor[4]=sum(W0>Wm[4])/a
Ap.valor[4]=sum(A0>Am[4])/a
Up.valor[4]=sum(U0>Um[4])/a
}
#Distribucion Gamma
if (any(x <= 0))
{
Dm[5]=0
Vm[5]=0
Wm[5]=0
Am[5]=0
Um[5]=0
parametros[8]=0
parametros[9]=0
Dp.valor[5]=0
Vp.valor[5]=0
Wp.valor[5]=0
Ap.valor[5]=0
Up.valor[5]=0
}
else
{
Newton=function(x,a0){
m=mean(x)
ml=mean(log(x))
a1=a0*(1+(digamma(a0)+log(m/a0)-ml)/(1-a0*trigamma(a0)))
}
alpha=function(x){
a00=(mean(x)^2)/var(x)
error=1
while(error>0.0001){
a1=Newton(x,a00)
error=abs(a1-a00)
a00=a1
}
a1
}
a.g=alpha(x)
b.g=a.g/m
prob5=pgamma(sort(x),a.g,b.g)
Dm[5]=Dn(prob5)
Vm[5]=Vn(prob5)
Wm[5]=Wn(prob5)
Am[5]=An(prob5)
Um[5]=Un(prob5)
parametros[8]=a.g
parametros[9]=b.g
#p.valor
D0=numeric(length(a))
V0=numeric(length(a))
W0=numeric(length(a))
A0=numeric(length(a))
U0=numeric(length(a))
for(i in 1:a)
{
rg=rgamma(n,a.g,b.g)
p=pgamma(sort(rg),a.g,b.g)
D0[i]=Dn(p)
V0[i]=Vn(p)
W0[i]=Wn(p)
A0[i]=An(p)
U0[i]=Un(p)
}
Dp.valor[5]=sum(D0>Dm[5])/a
Vp.valor[5]=sum(V0>Vm[5])/a
Wp.valor[5]=sum(W0>Wm[5])/a
Ap.valor[5]=sum(A0>Am[5])/a
Up.valor[5]=sum(U0>Um[5])/a
}
#Distribucion Cauchy
a.c=median(x)
b.c=IQR(x)/2
prob6=pcauchy(sort(x),a.c,b.c)
Dm[6]=Dn(prob6)
Vm[6]=Vn(prob6)
Wm[6]=Wn(prob6)
Am[6]=An(prob6)
Um[6]=Un(prob6)
parametros[10]=a.c
parametros[11]=b.c
#p.valor
D0=numeric(length(a))
V0=numeric(length(a))
W0=numeric(length(a))
A0=numeric(length(a))
U0=numeric(length(a))
for(i in 1:5000)
{
rc=rcauchy(n,a.c,b.c)
p=pcauchy(sort(rc),a.c,b.c)
D0[i]=Dn(p)
V0[i]=Vn(p)
W0[i]=Wn(p)
A0[i]=An(p)
U0[i]=Un(p)
}
Dp.valor[6]=sum(D0>Dm[6])/a
Vp.valor[6]=sum(V0>Vm[6])/a
Wp.valor[6]=sum(W0>Wm[6])/a
Ap.valor[6]=sum(A0>Am[6])/a
Up.valor[6]=sum(U0>Um[6])/a
#Distribucion Beta
#Beta
if (any(x > 1) || any(x < 0))
{
Dm[7]=0
Vm[7]=0
Wm[7]=0
Am[7]=0
Um[7]=0
parametros[12]=0
parametros[13]=0
Dp.valor[7]=0
Vp.valor[7]=0
Wp.valor[7]=0
Ap.valor[7]=0
Up.valor[7]=0
}
else
{
a.beta=-(m*v + m^3 - m^2)/v
b.beta=((m-1)*v + m^3 - 2*m^2 + m)/v
prob7=pbeta(sort(x),a.beta,b.beta)
Dm[7]=Dn(prob7)
Vm[7]=Vn(prob7)
Wm[7]=Wn(prob7)
Am[7]=An(prob7)
Um[7]=Un(prob7)
parametros[12]=a.beta
parametros[13]=b.beta
#p.valor
D0=numeric(length(a))
V0=numeric(length(a))
W0=numeric(length(a))
A0=numeric(length(a))
U0=numeric(length(a))
for(i in 1:a)
{
rb=rbeta(n,a.beta,b.beta)
p=pbeta(sort(rb),a.beta,b.beta)
D0[i]=Dn(p)
V0[i]=Vn(p)
W0[i]=Wn(p)
A0[i]=An(p)
U0[i]=Un(p)
}
Dp.valor[7]=sum(D0>Dm[7])/a
Vp.valor[7]=sum(V0>Vm[7])/a
Wp.valor[7]=sum(W0>Wm[7])/a
Ap.valor[7]=sum(A0>Am[7])/a
Up.valor[7]=sum(U0>Um[7])/a
}
#Distribucion Weibull
Newton.w=function(x,a0){
n=length(x)
a1= a0 - ( sum((x^a0)*log(x))/sum(x^a0) - 1/a0 - sum(log(x))/n )/
( sum((x^a0)*(log(x)^2))/sum(x^a0) - ((sum((x^a0)*log(x)))^2)/((sum(x^a0))^2) +1/a0^2 )
}
a.w=function(x){
a00=median(x)
error=1
while(error>0.0001){
a1=Newton.w(x,a00)
error=abs(a1-a00)
a00=a1
}
a1
}
if ( any(x < 0) || is.na(Newton.w(x,median(x)))==TRUE)
{
Dm[8]=0
Vm[8]=0
Wm[8]=0
Am[8]=0
Um[8]=0
parametros[14]=0
parametros[15]=0
Dp.valor[8]=0
Vp.valor[8]=0
Wp.valor[8]=0
Ap.valor[8]=0
Up.valor[8]=0
}
else
{
a.we=a.w(x)
b.we=(sum(x^a.we)/n)^(1/a.we)
prob8=pweibull(sort(x),a.we,b.we)
Dm[8]=Dn(prob8)
Vm[8]=Vn(prob8)
Wm[8]=Wn(prob8)
Am[8]=An(prob8)
Um[8]=Un(prob8)
parametros[14]=a.we
parametros[15]=b.we
#p.valor
D0=numeric(length(a))
V0=numeric(length(a))
W0=numeric(length(a))
A0=numeric(length(a))
U0=numeric(length(a))
for(i in 1:a)
{
rw=rweibull(n,a.we,b.we)
p=pweibull(sort(rw),a.we,b.we)
D0[i]=Dn(p)
V0[i]=Vn(p)
W0[i]=Wn(p)
A0[i]=An(p)
U0[i]=Un(p)
}
Dp.valor[8]=sum(D0>Dm[8])/a
Vp.valor[8]=sum(V0>Vm[8])/a
Wp.valor[8]=sum(W0>Wm[8])/a
Ap.valor[8]=sum(A0>Am[8])/a
Up.valor[8]=sum(U0>Um[8])/a
}
#Distribucion Ji-cuadrada
if ( any(x < 0) )
{
Dm[9]=0
Vm[9]=0
Wm[9]=0
Am[9]=0
Um[9]=0
parametros[16]=0
Dp.valor[9]=0
Vp.valor[9]=0
Wp.valor[9]=0
Ap.valor[9]=0
Up.valor[9]=0
}
else
{
k=floor(mean(x)+.5)
prob9=pchisq(sort(x),k)
Dm[9]=Dn(prob9)
Vm[9]=Vn(prob9)
Wm[9]=Wn(prob9)
Am[9]=An(prob9)
Um[9]=Un(prob9)
parametros[16]=k
#p.valor
D0=numeric(length(a))
V0=numeric(length(a))
W0=numeric(length(a))
A0=numeric(length(a))
U0=numeric(length(a))
for(i in 1:a)
{
rq=rchisq(n,k)
p=pchisq(sort(rq),k)
D0[i]=Dn(p)
V0[i]=Vn(p)
W0[i]=Wn(p)
A0[i]=An(p)
U0[i]=Un(p)
}
Dp.valor[9]=sum(D0>Dm[9])/a
Vp.valor[9]=sum(V0>Vm[9])/a
Wp.valor[9]=sum(W0>Wm[9])/a
Ap.valor[9]=sum(A0>Am[9])/a
Up.valor[9]=sum(U0>Um[9])/a
}
#Distribucion t
k.t=floor(2*v/(v-1)+.5)
if (k.t<=.5)
{
Dm[9]=0
Vm[9]=0
Wm[9]=0
Am[9]=0
Um[9]=0
parametros[16]=0
Dp.valor[9]=0
Vp.valor[9]=0
Wp.valor[9]=0
Ap.valor[9]=0
Up.valor[9]=0
}
else
{
prob10=pt(sort(x),k.t)
Dm[10]=Dn(prob10)
Vm[10]=Vn(prob10)
Wm[10]=Wn(prob10)
Am[10]=An(prob10)
Um[10]=Un(prob10)
parametros[17]=k.t
#p.valor
D0=numeric(length(a))
V0=numeric(length(a))
W0=numeric(length(a))
A0=numeric(length(a))
U0=numeric(length(a))
for(i in 1:a)
{
r.t=rt(n,k.t)
p=pt(sort(r.t),k.t)
D0[i]=Dn(p)
V0[i]=Vn(p)
W0[i]=Wn(p)
A0[i]=An(p)
U0[i]=Un(p)
}
Dp.valor[10]=sum(D0>Dm[10])/a
Vp.valor[10]=sum(V0>Vm[10])/a
Wp.valor[10]=sum(W0>Wm[10])/a
Ap.valor[10]=sum(A0>Am[10])/a
Up.valor[10]=sum(U0>Um[10])/a
}
return(list(rbind(Dm,Dp.valor, Vm,Vp.valor, Wm,Wp.valor, Am,Ap.valor, Um,Up.valor ),parametros))
}
else {stop("en tipo debe de ir 'continuo' o 'discreto'")}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.