iqp1mean <- function(x,mu,error=1e-8,w0=rep(1/length(x),length(x))){
###########  by Bin Huang, Mai Zhou  5/20/2000 version 1.0
###########  acquire original vales

  n <- length(x)            # smaple size
  xbar <- sum(x*w0)         # weighted sample mean

    mat0a <- diag(1/w0^2)          # original first entry
    mat0b <- 1/w0                    # original second entry
    mat0c <- cbind(rep(1,n), x)      # the constrain matrix, third entry
    mat0d <- c(0, mu-xbar)     # original constrain values, forth entry

    value0 <- solve.QP(mat0a, mat0b, mat0c, mat0d, meq=2) 
                                  # first calculation
    w <- w0 + value0$solution     # the weight vector

########## iteration

     matd <- c(0, 0)          # the new constrain values, forth entry
     diff <- 10               # original valus of accuracy

     while(diff > error) {               
        mata <- diag(1/w^2)      # the new first entry
        matb <- 1/w              # the new second entry
                                 # third entry do not need to be changed
        value <- solve.QP(mata, matb, mat0c, matd, meq=2)
        w <- w + value$solution   # update the weight function
        diff <- sum(value$solution^2)
     }

     ratio <- exp(sum(log(n*w)))
     weight <- w
     likevalue <- sum(log(w)) + value$value
 
  list(Ratio = ratio, Weight = weight, Likelihood = likevalue) 

#  return(weight, likevalue)
}
x<- c(28,-44,29,30,26,27,22,23,33,16,24,29,24,40,21,31,34,-2,25,19)
mu <- 20  ### let mu be a scalar not vector==seq(10, 26, 0.16)