目的

中央値付近に同値がある場合に比例配分によりより妥当な中央値を求める

使用法

my.median(x)

引数

x          データベクトルまたは級限界のベクトル（級限界のベクトルの要素数は，階級数より一つ多いことに注意）
f          x がデータベクトルの場合は NULL，
           x が級限界のベクトルの場合は度数ベクトル
accuracy   測定精度。省略時は普通のメディアン
 
ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/my_median.R", encoding="euc-jp")

# 中央値付近に同値がある場合に比例配分によりより妥当な中央値を求める
my.median <- function (      x,              # データベクトルまたは級限界のベクトル 
                        y=NULL,         # x がデータベクトルの場合は NULL
                                        # x が級限界のベクトルの場合は度数ベクトル
                        accuracy=0)     # 測定精度。省略時は普通のメディアン
{
        median.sub <- function(x)
        {
                n <- length(x)
                half <- (n + 1)/2
                if (n%%2 == 1) {
                        sort(x, partial = half)[half]
                }
                else {
                        sum(sort(x, partial = c(half, half + 1))[c(half, half + 1)])/2
                }
        }

        if (is.null(y)) {
                x <- x[!is.na(x)]
                median <- median.sub(x)
                if ((ntie <- length(x[x == median])) != 1 && accuracy > 0) {
                        x <- c(x[x != median], (median-(ntie+1)*accuracy/2/ntie)+1:ntie/ntie)
                        median <- median.sub(x)
                }
                return(median)
        }
        else {
                stopifnot(length(x)-1 == length(y))
                k <- length(y)
                csum <- cumsum(y)
                n <- csum[k]
                for (i in 1:k) {
                        if (csum[i] >= n/2) break
                }
                return(x[i]-accuracy/2+(n/2-csum[i-1])/y[i]*(x[i+1]-x[i]))
        }
}


使用例

x <- c(1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5)

median(x)	# R に用意されている関数を使う場合
my.median(x, accuracy=1)

出力結果例

> median(x)
[1] 3
> my.median(x, accuracy=1)
[1] 3.125

> x <- 1:6
> f <- c(2, 4, 8, 6, 2)
> my.median(x, f, accuracy=1)
[1] 3.125

 解説ページ