%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 画置信椭圆
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ellipse <-
function(hlaxa = 1, hlaxb = 1, theta = 0, xc = 0, yc = 0, npoints = 100)
{
        xp <- NULL
        yp <- NULL
        for(i in 0:npoints) {
                a <- (2 * pi * i)/npoints
                x <- hlaxa * cos(a)
                y <- hlaxb * sin(a)
                if(theta != 0) {
                        alpha <- angle(x, y)
                        rad <- sqrt(x^2 + y^2)
                        x <- rad * cos(alpha + theta)
                        y <- rad * sin(alpha + theta)
                }
                xp <- c(xp, x + xc)
                yp <- c(yp, y + yc)
        }
        lines(xp, yp, type = "l")
        invisible()
}

> angle <-
function(x, y)
{
        if(x > 0) {
                atan(y/x)
        }
        else {
                if(x < 0 & y != 0) {
                        atan(y/x) + sign(y) * pi
                }
                else {
                        if(x < 0 & y == 0) {
                                pi
                        }
                        else {
                                if(y != 0) {
                                        (sign(y) * pi)/2
                                }
                                else {
                                        NA
                                }
                        }
                }
        }
}

> ellipsem <-
function (mu, amat, c2, npoints = 100) 
{
    if (dim(amat) == c(2, 2)) {
        eamat <- eigen(amat)
        hlen <- sqrt(c2/eamat$val)
        theta <- angle(eamat$vec[1, 1], eamat$vec[2, 1])
        ellipse(hlen[1], hlen[2], theta, mu[1], mu[2], npoints = npoints)
    }
    invisible()
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> A<-read.csv("Hemophilia.csv")
> A1<-A[1:30,2:3]
> plot(A1)
> points(mean(A1[,1]),mean(A1[,2]), pch="+")
> n<-30
> p<-2
> ellipsem(mean(A1), solve(var(A1)), ((n-1)*p/(n-p))*qf(0.95,p,n-p)/n）