Nothing
Joe.Markov.MLE <-
function(Y,k=3,D=1,plot=TRUE,GOF=FALSE, method = "nlm"){
n=length(Y) ##sample size##
log.L = function(par){
mu = par[1]
sigma = exp(par[2])
alpha = exp(par[3])+1
g_yt = dnorm(Y, mean = mu, sd = sigma)
G_yt = pnorm(Y, mean = mu, sd = sigma)
A = function(u, v){
(1-u)^alpha + (1-v)^alpha - (1-u)^alpha*(1-v)^alpha
}
C_11 = function(u,v){
alpha * A(u,v)^(1/alpha-1) * (1-u)^(alpha-1) * (1-v)^(alpha-1) +
(alpha-1) * A(u,v)^(1/alpha-2) * (1-u)^(alpha-1) * (1-v)^(alpha-1) *
(1-(1-v)^alpha) * (1-(1-u)^alpha)
}
res = sum( log(g_yt) ) + sum( log( C_11( G_yt[-n], G_yt[-1] ) ) )
return(-res)
}
if (method == "Newton"){
G=function(y,mu,sigma){pnorm((y-mu)/sigma)} #G function
g=function(y,mu,sigma){dnorm((y-mu)/sigma)} #g function
##############Log-likelihood function #################
L_function=function(theta){
mu=theta[1]
sigma=theta[2]
alpha=theta[3]
U_t_1=G(Y[1:n-1],mu,sigma);U_t=G(Y[2:n],mu,sigma)
u_t_1=g(Y[1:n-1],mu,sigma);u_t=g(Y[2:n],mu,sigma)
A=(1-U_t_1)^alpha+(1-U_t)^alpha-((1-U_t_1)^alpha)*((1-U_t)^alpha)
Z=log(alpha-1+A)+(alpha-1)*log(1-U_t_1)+(alpha-1)*log(1-U_t)+(1/alpha-2)*log(A)
ZZ=log(g(Y[1:n],mu,sigma)/sigma)
-(sum(Z)+sum(ZZ))/n
}
#res=nlm(L_function,c(mean(Y),sd(Y),1),hessian=TRUE)
#est=res$estimate
##############dL/dmu############################
dL_dmu=function(mu,sigma,alpha){
U_t_1=G(Y[1:n-1],mu,sigma);U_t=G(Y[2:n],mu,sigma)
u_t_1=g(Y[1:n-1],mu,sigma);u_t=g(Y[2:n],mu,sigma)
A=(1-U_t_1)^alpha+(1-U_t)^alpha-((1-U_t_1)^alpha)*((1-U_t)^alpha)
A_mu1=alpha*(1-U_t_1)^(alpha-1)*(1-(1-U_t)^alpha)*u_t_1/sigma
A_mu2=alpha*(1-U_t)^(alpha-1)*(1-(1-U_t_1)^alpha)*u_t/sigma
A_mu=A_mu1+A_mu2
L1=(Y[1:n]-mu)/sigma^2
L2=A_mu/(alpha-1+A)
L3=(alpha-1)/sigma*( u_t_1/(1-U_t_1)+u_t/(1-U_t) )
L4=(1/alpha-2)*A_mu/A
return((sum(L1)+sum(L2)+sum(L3)+sum(L4))/n)
}
############dL/dsigma###########################
dL_dsigma=function(mu,sigma,alpha){
U_t_1=G(Y[1:n-1],mu,sigma);U_t=G(Y[2:n],mu,sigma)
u_t_1=g(Y[1:n-1],mu,sigma);u_t=g(Y[2:n],mu,sigma)
A=(1-U_t_1)^alpha+(1-U_t)^alpha-((1-U_t_1)^alpha)*((1-U_t)^alpha)
A_s1=alpha*(1-U_t_1)^(alpha-1)*(1-(1-U_t)^alpha)*(Y[1:n-1]-mu)*u_t_1/sigma^2
A_s2=alpha*(1-U_t)^(alpha-1)*(1-(1-U_t_1)^alpha)*(Y[2:n]-mu)*u_t/sigma^2
A_s=A_s1+A_s2
L1=(Y[1:n]-mu)^2/sigma^3
L2=A_s/(alpha-1+A)
L3=(alpha-1)/sigma^2*( (Y[1:n-1]-mu)*u_t_1/(1-U_t_1)+(Y[2:n]-mu)*u_t/(1-U_t) )
L4=(1/alpha-2)*A_s/A
return(-1/sigma+(sum(L1)+sum(L2)+sum(L3)+sum(L4))/n)
}
############dL/dalpha#############################
dL_dalpha=function(mu,sigma,alpha){
U_t_1=G(Y[1:n-1],mu,sigma);U_t=G(Y[2:n],mu,sigma)
u_t_1=g(Y[1:n-1],mu,sigma);u_t=g(Y[2:n],mu,sigma)
A=(1-U_t_1)^alpha+(1-U_t)^alpha-((1-U_t_1)^alpha)*((1-U_t)^alpha)
A_a1=(1-U_t_1)^alpha*log(1-U_t_1)+(1-U_t)^alpha*log(1-U_t)
A_a2=(1-U_t_1)^alpha*(1-U_t)^alpha*log((1-U_t_1)*(1-U_t))
A_a=A_a1-A_a2
L1=log(1-U_t_1)+log(1-U_t)
L2=(1+A_a)/(alpha-1+A)
L3=-log(A)/alpha^2
L4=(1/alpha-2)*A_a/A
return((sum(L1)+sum(L2)+sum(L3)+sum(L4))/n)
}
############F function##############################
F=function(mu,sigma,alpha){
c(dL_dmu(mu,sigma,alpha),dL_dsigma(mu,sigma,alpha),dL_dalpha(mu,sigma,alpha))
}
#F(est[1],est[2],est[3])
############d^2L/dmu^2 (in progress) ###########################
d2L_dmu2=function(mu,sigma,alpha){
U_t_1=G(Y[1:n-1],mu,sigma);U_t=G(Y[2:n],mu,sigma)
u_t_1=g(Y[1:n-1],mu,sigma);u_t=g(Y[2:n],mu,sigma)
h_t_1=u_t_1/(1-U_t_1);h_t=u_t/(1-U_t)
A=(1-U_t_1)^alpha+(1-U_t)^alpha-((1-U_t_1)^alpha)*((1-U_t)^alpha)
A_mu1=alpha*(1-U_t_1)^(alpha-1)*(1-(1-U_t)^alpha)*u_t_1/sigma
A_mu2=alpha*(1-U_t)^(alpha-1)*(1-(1-U_t_1)^alpha)*u_t/sigma
A_mu=A_mu1+A_mu2
B_mm1=(alpha-1)/(1-U_t_1)*u_t_1*(1-(1-U_t)^alpha)
B_mm2=-alpha*(1-U_t)^(alpha-1)*u_t
B_mm3=(1-(1-U_t)^alpha)*(Y[1:n-1]-mu)/sigma
A_mm1=alpha*(1-U_t_1)^(alpha-1)*(B_mm1+B_mm2+B_mm3)*u_t_1/sigma^2
C_mm1=(alpha-1)/(1-U_t)*u_t*(1-(1-U_t_1)^alpha)
C_mm2=-alpha*(1-U_t_1)^(alpha-1)*u_t_1
C_mm3=(1-(1-U_t_1)^alpha)*(Y[2:n]-mu)/sigma
A_mm2=alpha*(1-U_t)^(alpha-1)*(C_mm1+C_mm2+C_mm3)*u_t/sigma^2
A_mm=A_mm1+A_mm2
L1=A_mm/(alpha-1+A)-( A_mu/(alpha-1+A) )^2
L2=(alpha-1)/sigma^2*( (Y[1:n-1]-mu)*h_t_1/sigma-h_t_1^2 )
L3=(alpha-1)/sigma^2*( (Y[2:n]-mu)*h_t/sigma-h_t^2 )
L4=(1/alpha-2)*( A_mm/A-(A_mu/A)^2 )
return((sum(L1)+sum(L2)+sum(L3)+sum(L4))/n-1/sigma^2)
}
############d^2L/dsigma^2 (in progress)#####################
d2L_dsigma2=function(mu,sigma,alpha){
U_t_1=G(Y[1:n-1],mu,sigma);U_t=G(Y[2:n],mu,sigma)
u_t_1=g(Y[1:n-1],mu,sigma);u_t=g(Y[2:n],mu,sigma)
h_t_1=u_t_1/(1-U_t_1);h_t=u_t/(1-U_t)
A=(1-U_t_1)^alpha+(1-U_t)^alpha-((1-U_t_1)^alpha)*((1-U_t)^alpha)
A_mu1=alpha*(1-U_t_1)^(alpha-1)*(1-(1-U_t)^alpha)*u_t_1/sigma
A_mu2=alpha*(1-U_t)^(alpha-1)*(1-(1-U_t_1)^alpha)*u_t/sigma
A_mu=A_mu1+A_mu2
A_s1=alpha*(1-U_t_1)^(alpha-1)*(1-(1-U_t)^alpha)*(Y[1:n-1]-mu)*u_t_1/sigma^2
A_s2=alpha*(1-U_t)^(alpha-1)*(1-(1-U_t_1)^alpha)*(Y[2:n]-mu)*u_t/sigma^2
A_s=A_s1+A_s2
B_ms1=(alpha-1)/(1-U_t_1)*u_t_1*(1-(1-U_t)^alpha)*(Y[1:n-1]-mu)
B_ms2=-alpha*(1-U_t)^(alpha-1)*u_t*(Y[1:n-1]-mu)-1+(1-U_t)^alpha
B_ms3=(1-(1-U_t)^alpha)*(Y[1:n-1]-mu)^2/sigma
A_ms1=alpha*(1-U_t_1)^(alpha-1)*(B_ms1+B_ms2+B_ms3)*u_t_1/sigma^3
C_ms1=(alpha-1)/(1-U_t)*u_t*(1-(1-U_t_1)^alpha)*(Y[2:n]-mu)
C_ms2=-alpha*(1-U_t_1)^(alpha-1)*u_t_1*(Y[2:n]-mu)-1+(1-U_t_1)^alpha
C_ms3=(1-(1-U_t_1)^alpha)*(Y[2:n]-mu)^2/sigma
A_ms2=alpha*(1-U_t)^(alpha-1)*(C_ms1+C_ms2+C_ms3)*u_t/sigma^3
A_ms=A_ms1+A_ms2
L1=-2*(Y[1:n]-mu)/sigma^3
L2=A_ms/(alpha-1+A)-A_mu*A_s/(alpha-1+A)^2
L3=-(alpha-1)*( h_t_1/sigma^2+h_t/sigma^2 )
L4=(alpha-1)*( (Y[1:n-1]-mu)/sigma^3*( (Y[1:n-1]-mu)/sigma*h_t_1-u_t_1^2 ) )
L5=(alpha-1)*( (Y[2:n]-mu)/sigma^3*( (Y[2:n]-mu)/sigma*h_t-u_t^2 ) )
#L6=(1/alpha-2)*( A_ss/A-(A_s/A)^2 )
return(1/sigma^2+( sum(L1)+sum(L2)+sum(L3)+sum(L4)+sum(L5) )/n)
}
d2L_dsigma2=function(mu,sigma,alpha){
h=0.0000001
(dL_dsigma(mu,sigma+h,alpha)-dL_dsigma(mu,sigma,alpha))/h
}
#d2L_dsigma2(est[1],est[2],est[3])
#-res$hessian[2,2]
#h=0.0000001
#-(L_function(est+2*c(0,h,0))-2*L_function(est+c(0,h,0))+L_function(est))/h^2
############d^2L/dalpha^2#################################
d2L_dalpha2=function(mu,sigma,alpha){
U_t_1=G(Y[1:n-1],mu,sigma);U_t=G(Y[2:n],mu,sigma)
u_t_1=g(Y[1:n-1],mu,sigma);u_t=g(Y[2:n],mu,sigma)
A=(1-U_t_1)^alpha+(1-U_t)^alpha-((1-U_t_1)^alpha)*((1-U_t)^alpha)
A_a1=(1-U_t_1)^alpha*log(1-U_t_1)+(1-U_t)^alpha*log(1-U_t)
A_a2=(1-U_t_1)^alpha*(1-U_t)^alpha*log((1-U_t_1)*(1-U_t))
A_a=A_a1-A_a2
A_aa1=(1-U_t_1)^alpha*(log(1-U_t_1))^2+(1-U_t)^alpha*(log(1-U_t))^2
A_aa2=(1-U_t_1)^alpha*(1-U_t)^alpha*(log((1-U_t_1)*(1-U_t)))^2
A_aa=A_aa1-A_aa2
L1=A_aa/(alpha-1+A)-((1+A_a)/(alpha-1+A))^2
L2=2*log(A)/alpha^3-2*A_a/A/alpha^2
L3=(1/alpha-2)*(A_aa/A-(A_a/A)^2)
return((sum(L1)+sum(L2)+sum(L3))/n)
}
############d^2L/dmudsigma (in progress) #################################
d2L_dmudsigma=function(mu,sigma,alpha){
U_t_1=G(Y[1:n-1],mu,sigma);U_t=G(Y[2:n],mu,sigma)
u_t_1=g(Y[1:n-1],mu,sigma);u_t=g(Y[2:n],mu,sigma)
h_t_1=u_t_1/(1-U_t_1);h_t=u_t/(1-U_t)
A=(1-U_t_1)^alpha+(1-U_t)^alpha-((1-U_t_1)^alpha)*((1-U_t)^alpha)
A_s1=alpha*(1-U_t_1)^(alpha-1)*(1-(1-U_t)^alpha)*(Y[1:n-1]-mu)*u_t_1/sigma^2
A_s2=alpha*(1-U_t)^(alpha-1)*(1-(1-U_t_1)^alpha)*(Y[2:n]-mu)*u_t/sigma^2
A_s=A_s1+A_s2
B_ss1=(alpha-1)/(1-U_t_1)*u_t_1*(1-(1-U_t)^alpha)*(Y[1:n-1]-mu)^2
B_ss2=-alpha*(1-U_t)^(alpha-1)*u_t*(Y[2:n]-mu)*(Y[1:n-1]-mu)
B_ss3=(1-(1-U_t)^alpha)*(Y[1:n-1]-mu)^3/sigma
A_ss1=alpha*(1-U_t_1)^(alpha-1)*(B_ss1+B_ss2+B_ss3)*u_t_1/sigma^4
C_ss1=(alpha-1)/(1-U_t)*u_t*(1-(1-U_t_1)^alpha)*(Y[2:n]-mu)^2
C_ss2=-alpha*(1-U_t_1)^(alpha-1)*u_t_1*(Y[1:n-1]-mu)*(Y[2:n]-mu)
C_ss3=(1-(1-U_t_1)^alpha)*(Y[2:n]-mu)^3/sigma
A_ss2=alpha*(1-U_t)^(alpha-1)*(C_ss1+C_ss2+C_ss3)*u_t/sigma^4
A_ss=A_ss1+A_ss2
}
d2L_dmudsigma=function(mu,sigma,alpha){
h=0.0000001
(dL_dsigma(mu+h,sigma,alpha)-dL_dsigma(mu,sigma,alpha))/h
}
#d2L_dmudsigma(est[1],est[2],est[3])
#-res$hessian[1,2]
#h=0.0000001
############d^2L/dmudalpha#################################
d2L_dmudalpha=function(mu,sigma,alpha){
U_t_1=G(Y[1:n-1],mu,sigma);U_t=G(Y[2:n],mu,sigma)
u_t_1=g(Y[1:n-1],mu,sigma);u_t=g(Y[2:n],mu,sigma)
A=(1-U_t_1)^alpha+(1-U_t)^alpha-((1-U_t_1)^alpha)*((1-U_t)^alpha)
A_a1=(1-U_t_1)^alpha*log(1-U_t_1)+(1-U_t)^alpha*log(1-U_t)
A_a2=(1-U_t_1)^alpha*(1-U_t)^alpha*log((1-U_t_1)*(1-U_t))
A_a=A_a1-A_a2
A_mu1=alpha*(1-U_t_1)^(alpha-1)*(1-(1-U_t)^alpha)*u_t_1/sigma
A_mu2=alpha*(1-U_t)^(alpha-1)*(1-(1-U_t_1)^alpha)*u_t/sigma
A_mu=A_mu1+A_mu2
B1=1-(1-U_t)^alpha+alpha*log(1-U_t_1)*(1-(1-U_t)^alpha)
B2=-alpha*(1-U_t)^alpha*log(1-U_t)
A_am1=(1-U_t_1)^(alpha-1)*(B1+B2)*u_t_1/sigma
C1=1-(1-U_t_1)^alpha+alpha*log(1-U_t)*(1-(1-U_t_1)^alpha)
C2=-alpha*(1-U_t_1)^alpha*log(1-U_t_1)
A_am2=(1-U_t)^(alpha-1)*(C1+C2)*u_t/sigma
A_am=A_am1+A_am2
L1=A_am/(alpha-1+A)-(1+A_a)*A_mu/(alpha-1+A)^2
L2=1/sigma*( u_t_1/(1-U_t_1)+u_t/(1-U_t) )-A_mu/alpha^2/A
L3=(1/alpha-2)*( A_am/A-A_mu*A_a/A^2 )
return((sum(L1)+sum(L2)+sum(L3))/n)
}
############d^2L/dsigmadalpha (in progress) #################################
d2L_dsigmadalpha=function(mu,sigma,alpha){
h=0.0000001
(dL_dsigma(mu,sigma,alpha+h)-dL_dsigma(mu,sigma,alpha))/h
}
#############Jacobian function######################################
Ja=function(mu,sigma,alpha){
AA=matrix( c(d2L_dmu2(mu,sigma,alpha),d2L_dmudsigma(mu,sigma,alpha),d2L_dmudalpha(mu,sigma,alpha),d2L_dmudsigma(mu,sigma,alpha),d2L_dsigma2(mu,sigma,alpha),d2L_dsigmadalpha(mu,sigma,alpha),d2L_dmudalpha(mu,sigma,alpha),d2L_dsigmadalpha(mu,sigma,alpha),d2L_dalpha2(mu,sigma,alpha)),3,3 )
return(AA)
}
###############Multivariate Newton Raphson#####################
X=matrix(,1,3)
tau=cor(Y[1:n-1],Y[2:n],method="kendall")
alpha_est=-2*tau/(tau-1)
alpha_est=2
X[1,]=c(mean(Y),sd(Y),alpha_est) #initial value
i=2
Ran.num=1
repeat{
Z=X
X=matrix(,i,3)
X[1:i-1,]=Z[1:i-1,]
##
Aa=Ja(X[i-1,1],X[i-1,2],X[i-1,3])
Ainv11=Aa[2,2]*Aa[3,3]-Aa[3,2]*Aa[2,3]
Ainv12=Aa[1,2]*Aa[3,3]-Aa[3,2]*Aa[1,3]
Ainv13=Aa[1,2]*Aa[2,3]-Aa[2,2]*Aa[1,3]
Ainv21=Aa[2,1]*Aa[3,3]-Aa[3,1]*Aa[2,3]
Ainv22=Aa[1,1]*Aa[3,3]-Aa[3,1]*Aa[1,3]
Ainv23=Aa[1,1]*Aa[2,3]-Aa[1,3]*Aa[2,1]
Ainv31=Aa[2,1]*Aa[3,2]-Aa[3,1]*Aa[2,2]
Ainv32=Aa[1,1]*Aa[3,2]-Aa[1,2]*Aa[3,1]
Ainv33=Aa[1,1]*Aa[2,2]-Aa[1,2]*Aa[2,1]
Ainv=matrix(c(Ainv11,-Ainv21,Ainv31,-Ainv12,Ainv22,-Ainv32,Ainv13,-Ainv23,Ainv33),3,3)/det(Aa)
##
X[i,]=X[i-1,]-Ainv%*%F(X[i-1,1],X[i-1,2],X[i-1,3])
if(1*is.nan(X)[i,1]==1){
X=matrix(,2,3)
X[1,]=c(mean(Y),sd(Y),alpha_est+runif(1,-D,D)) #initial value
Ran.num=Ran.num+1
i=1
}else if(abs(X[i,1]-X[i-1,1])<0.0001&abs(X[i,2]-X[i-1,2])<0.0001&abs(X[i,3]-X[i-1,3])<0.0001&abs(X[i,3]-alpha_est)>5){
X=matrix(,2,3)
X[1,]=c(mean(Y),sd(Y),alpha_est+runif(1,-D,D)) #initial value
Ran.num=Ran.num+1
i=1
}else if(abs(X[i,1]-X[i-1,1])<0.0001&abs(X[i,2]-X[i-1,2])<0.0001&abs(X[i,3]-X[i-1,3])<0.0001&X[i,2]>0&abs(X[i,3]-alpha_est)<5){break
}else if(abs(X[i,1]-X[i-1,1])<0.0001&abs(X[i,2]-X[i-1,2])<0.0001&abs(X[i,3]-X[i-1,3])<0.0001&X[i,2]<0){
X=matrix(,2,3)
X[1,]=c(mean(Y),sd(Y),alpha_est+runif(1,-D,D)) #initial value
Ran.num=Ran.num+1
i=1
}else if(abs(X[i,1]-X[i-1,1])>10^10&abs(X[i,2]-X[i-1,2])>10^10&abs(X[i,3]-X[i-1,3])>10^10){
X=matrix(,2,3)
X[1,]=c(mean(Y),sd(Y),alpha_est+runif(1,-D,D)) #initial value
Ran.num=Ran.num+1
i=1
}
if(Ran.num>=100){break}
i=i+1
}
mle.res=X[length(X[,1]),]
if(Ran.num>=10){ mle.res=c(mean(Y),sd(Y),alpha_est) }
mu.hat=mle.res[1]
sigma.hat=mle.res[2]
alpha.hat=mle.res[3]
Gradient=n^3*F(mu.hat,sigma.hat,alpha.hat)
Hessian=n*Ja(mu.hat,sigma.hat,alpha.hat)
SE_mu = SE_sigma = SE_alpha = lower_mu = upper_mu = lower_sigma = upper_sigma =
lower_alpha = upper_alpha = NA
}
if (method == "nlm"){
tau_0 = cor(Y[-n],Y[-1],method="kendall")
initial = c(mean(Y),log(sd(Y)), log(2))
mle.res = nlm(f = log.L, p = initial, hessian = TRUE)
mu.hat=mle.res$estimate[1]
sigma.hat=exp(mle.res$estimate[2])
alpha.hat=exp(mle.res$estimate[3])+1
Gradient=-mle.res$gradient
Hessian=-mle.res$hessian
inverse_Hessian=solve(Hessian,tol=10^(-50))
SE_mu = sqrt(-inverse_Hessian[1,1])
SE_sigma = sqrt(-inverse_Hessian[2,2])*sigma.hat
SE_alpha = sqrt(-inverse_Hessian[3,3])*(alpha.hat-1)
lower_mu = mu.hat-1.96*SE_mu
upper_mu = mu.hat+1.96*SE_mu
lower_sigma=sigma.hat*exp(-1.96*SE_sigma/sigma.hat)
upper_sigma=sigma.hat*exp(1.96*SE_sigma/sigma.hat)
lower_alpha=(alpha.hat-1)*exp(-1.96*SE_alpha/(alpha.hat-1))+1
upper_alpha=(alpha.hat-1)*exp(1.96*SE_alpha/(alpha.hat-1))+1
}
UCL=mu.hat+k*sigma.hat
LCL=mu.hat-k*sigma.hat
CL = c(Center = mu.hat, Lower = LCL, Upper = UCL)
result.mu = c(estimate = mu.hat, SE = SE_mu, Lower = lower_mu, Upper = upper_mu)
result.sigma = c(estimate = sigma.hat, SE = SE_sigma, Lower = lower_sigma, Upper = upper_sigma)
result.alpha = c(estimate = alpha.hat, SE = SE_alpha, Lower = lower_alpha, Upper = upper_alpha)
####### Plot Control Chart #######
if(plot==TRUE){
Min=min(min(Y),LCL)
Max=max(max(Y),UCL)
ts.plot(Y,type="b",ylab="Y",ylim=c(Min,Max))
abline(h=mu.hat)
abline(h=UCL,lty="dotted",lwd=2)
abline(h=LCL,lty="dotted",lwd=2)
text(0,LCL+(mu.hat-LCL)*0.1,"LCL")
text(0,UCL-(UCL-mu.hat)*0.1,"UCL")
}
out_control=which( (Y<LCL)|(UCL<Y) )
if(length(out_control)==0){out_control="NONE"}
### Goodness-of-fit ###
F_par=pnorm( (sort(Y)-mu.hat)/sigma.hat )
F_emp=1:n/n
CM.test=sum( (F_emp-F_par)^2 )
KS.test=max( abs( F_emp-F_par ) )
if(GOF==TRUE){
plot(F_emp,F_par,xlab="F_empirical",ylab="F_parametric",xlim=c(0,1),ylim=c(0,1))
lines(x = c(0,1), y = c(0,1))
}
return(
list(mu=result.mu, sigma = result.sigma, alpha = result.alpha,
Control_Limit = CL, out_of_control=out_control,
Gradient=Gradient,Hessian=Hessian,Eigenvalue_Hessian=eigen(Hessian)$value,
CM.test=CM.test,KS.test=KS.test, log_likelihood = -log.L(c(mu.hat, log(sigma.hat), log(alpha.hat-1))) )
)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.