# ConReg by Shott and Habtzghi (2016, 2019) # The ConRge function uses Quadpgrog package in R to fit Constrained Regression # Objective: Computes the estimated proportions of the four stages # Input : know standards C's values and Y # This function accommodates 2,3 or 4 number of size intervals ConReg <- function(Y,C1,C2,C3=NA,C4=NA) { X1=cbind(C1,C2,C3,C4) k=4 if(is.na(C4[1])|is.na(C3[1])) { X1=cbind(C1,C2,C3) k=3 } if(is.na(C3[1])& is.na(C4[1])) { X1=cbind(C1,C2) k=2 } # design matrix combination of the known values X <- matrix(X1, ncol=k) Dmat <- crossprod(X, X) # H=XtX The symmetric matrix H is often called the Hessian. dvec <- crossprod(Y, X) # D=YtX Amat <- cbind(1, diag(NROW(Dmat))) # constraint matrix under which we want to #minimize the quadratic function bvec <- c(1, rep(0,NROW(Dmat))) # vector holding the values of b_0 # An indicator for inequality constraints meq <- 1 # The Quadpgrog package in R library(quadprog) #The solve.QP function implements the dual method of Goldfarb and Idnani (1982, 1983) #for solving quadratic programming problems # of the form min(-d^T b + 1/2 b^T H b) with the constraints A^T b >= b_0. res <- solve.QP(Dmat, dvec, Amat, bvec, meq) # Estimates coefficient coef1=zapsmall(res\$solution) phat1=coef1[1] phat2=coef1[2] phat3=coef1[3] phat4=coef1[4] estimated=c(phat1,phat2,phat3,phat4) if(is.na(phat4)|is.na(phat3)) { estimated=c(phat1,phat2,phat3) } if(is.na(phat4)& is.na(phat3)) { estimated=c(phat1,phat2) } tt=round(estimated, 3) return(tt) } #Example #Known standards C1=c(0.7059,0.8759,0.9232,0.9591,0.9713,0.9813,0.9878,0.9950,0.9993,1.0000) C2=c(0.7075,0.8834,0.9330,0.9737,0.9841,0.9918,0.9944,0.9996,1.0000,1.0000) C3=c(0.7437,0.9246,0.9663,0.9894,0.9952,0.9979,0.9992,1.0000,1.0000,1.0000) C4=c(0.7876,0.9556,0.9800,0.9946,0.9990,1.0000,1.0000,1.0000,1.0000,1.0000) #case y1=c(0.7643,0.9000,0.9143,0.9429,0.9714,0.9714,0.9714,0.9929,1.000,1.000) y2=c(0.7500,0.8929,0.9286,0.9643,1.0000,1.0000,1.000,1.000,1.000,1.000) y3=c(0.6807,0.8675,0.9036,0.9438,0.9900,0.9799,0.9598,1.0000,1.000,1.000) y4=c(0.6916,0.8703,0.9290,0.9670,0.9927,0.9914,0.9816,0.9988,1.000,1.00) y5=c(0.7340,0.8951,0.9417,0.9864,0.9942,0.9922,0.9922,1.00,1.000,1.000) y6=c(0.7343,0.9004,0.9520,0.9889,1.000,1.000,0.9963,1.000,1.000,1.000) y7=c(0.8133,0.9533,0.9800,0.9933,1.000,1.000,1.00,1.00,1.000,1.000) y8=c(0.7660,0.9503,0.9834,0.9959,1.000,1.000,0.9979,1.000,1.00, 1.000) #1a, 11b y9=c(0.70219,0.87461,0.92685,0.96343,0.98015,0.98851,0.98955,0.99791,1.00000,1.00000) #2a, 11b y10=c(0.69349,0.87101,0.92899,0.96686,0.98225,0.99172,0.99290,0.99882,1.00000,1.00000) #1a, 11b, 9c y11=c(0.71332,0.88179,0.93207,0.97147,0.98438,0.98981,0.99117,0.99864,1.00000,1.00000) #1a, 11b, 12c y12=c(0.70928,0.88029,0.93241,0.96906,0.98371,0.99104,.99186,0.99837,1.00000,1.00000) out=ConReg(y1,C1,C2,C3,C4) out #Example #Known standards C1=c(0.7059,0.8759,0.9232,0.9591,0.9713,0.9813,0.9878,0.9950,0.9993,1.0000) C2=c(0.7075,0.8834,0.9330,0.9737,0.9841,0.9918,0.9944,0.9996,1.0000,1.0000) C3=c(0.7437,0.9246,0.9663,0.9894,0.9952,0.9979,0.9992,1.0000,1.0000,1.0000) y1=c(0.7643,0.9000,0.9143,0.9429,0.9714,0.9714,0.9714,0.9929,1.000,1.000) out=ConReg(y1,C1,C2,C3) out # Weibull Distribution d=c(.125,.25,.375,.5,.625,.75,.875,1.375) C211=c(.692,.87,.929,.967,.982,.991,.993,.999) d=c(.125,.25,.375,.5,.625) C1=c(0.7082,0.8755,0.9225,0.9587,0.9710,0.9812,0.9877,0.9949) C2=c(0.7073,0.8832,0.9329,0.9735,0.9839,0.9916,0.9942,0.9994) C3=c(0.7437,0.9246,0.9663,0.9894,0.9952,0.9979,0.9992) C4=c(0.7874,0.9572,0.9791,0.9944,0.9990) F=C1 x=log(d) y=log(log(1/(1-F))) n=length(x) B1=(n*sum(x*y)-sum(x)*sum(y))/(n*sum(x^2)-sum(x)^2) B0=(sum(x^2)*sum(y)-sum(x)*sum(x*y))/(n*sum(x^2)-sum(x)^2) lambda=exp(-B0/B1) t=c(.125,.25,.375,.5,.625,.75,.875,1.375) beta=B1 F1=1-exp(-(t/lambda)^beta) cbind(d,F, F1) Rsq=cor(F,F1)^2