R циркуляр - Центр распределения - Среднее

Question

У меня есть набор данных с количеством вхождений по месяцам, и я хотел бы рассчитать центр распределения / средний месяц. Если возможно, я также хотел бы иметь доверительные интервалы.

Я прочитал руководства для циркулярных и CircStats и посмотрел на подобные вопросы здесь и здесь.

Мне удалось получить, казалось бы, разумный результат в некоторых случаях, но не в других, и я еще не выяснил, как рассчитать доверительный интервал.

Чтобы проиллюстрировать мою точку зрения, вот некоторые фиктивные данные:

library(CircStats)

# The number of observations by month (Jan-Dec):
obsMonths1 <- c(12,15,1,2,3,1,1,4,1,2,7,1) 
obsMonths2 <- c(1,1,1,1,2,10,11,2,1,1,2,1)

# Convert data to radians:
obsRadians1 <- (obsMonths1/12*2)*pi
obsRadians2 <-(obsMonths2/12*2)*pi

# Calculate circular mean:
mean1 <- circ.mean(obsRadians1-1)#assume January is 0
mean2 <- circ.mean(obsRadians2-1)#assume January is 0

# Convert radians to months:
mean1*12/(2*pi)+12 
mean2*12/(2*pi)+12 

Для первого набора наблюдений ответ кажется разумным, но для второго набора наблюдений это должен быть июль-август.

Answer 1

Похоже, ваш код не имеет проблем. Вы спросили "доверительный интервал", но не похоже, чтобы вы проводили какие-либо статистические тесты или процедуры начальной загрузки. Я не уверен, как бы вы рассчитали доверительный интервал без этих результатов в первую очередь. Если вы уверены, что поступили правильно, возможно, вы захотите поднять этот вопрос в сообществе StackExchange, специализирующемся на математике и статистике.

Answer 2

Я согласен, что ваши расчеты в основном верны; Я думаю, что это твоя интуиция не так, как я иллюстрирую ниже с некоторыми картинками.

library(CircStats)
# The number of observations by month (Jan-Dec):
obsMonths1 <- c(12,15,1,2,3,1,1,4,1,2,7,1) 
obsMonths2 <- c(1,1,1,1,2,10,11,2,1,1,2,1)

Я сделал преобразование в и из радиан немного по-другому; по модулю (%%) оператор автоматически конвертирует 12 в ноль. Я добавил 11, чтобы сделать январь ==0, но сохранить все позитивно...

to_rad <- function(x) (x+11 %% 12)/12*2*pi
## check results
stopifnot(to_rad(1)==0,to_rad(7)==pi,to_rad(4)==pi/2,to_rad(10)==3*pi/2)

И преобразовать обратно:

from_rad <- function(x) (12/(2*pi)*x)+1
## check round-trip with an arbitrary number
stopifnot(isTRUE(all.equal(from_rad(to_rad(7.931)),7.931)))

Конверсия:

(m1 <- from_rad(circ.mean(to_rad(obsMonths1))))  ## 1.77
(m2 <- from_rad(circ.mean(to_rad(obsMonths2))))  ## 0.93

Самозагрузочный код:

bootquant <- function(x,n=1000,alpha=0.05) {
    bootsamp  <- replicate(n,
                           from_rad(circ.mean(to_rad(sample(x,replace=TRUE)))))
    qq <- quantile(bootsamp,c(alpha/2,1-alpha/2))
    names(qq) <- c("lwr","upr")
    return(qq)
}
(bq1 <- bootquant(obsMonths1))
##      lwr      upr 
## 1.076794 2.670130 
(bq2 <- bootquant(obsMonths2))
##      lwr      upr 
## 0.231873 1.414766 

Я не уверен, что доверял бы начальной загрузке для таких небольших наборов данных; Вы также можете проверить ?circ.disp функция от CircStats...

library(ggplot2)
dd <- data.frame(OM1=obsMonths1,OM2=obsMonths2)
ggplot(dd,aes(x=OM1,y=1))+stat_sum()+coord_polar()+
    scale_x_continuous(limits=c(0,12), breaks=c(0,3,6,9,12))+
    annotate(geom="point",y=1,x=m1,colour="red")+
    annotate(geom="segment",x=bq1[["lwr"]],xend=bq1[["upr"]],y=1,yend=1,colour="red")

ggplot(dd,aes(x=OM2,y=1))+stat_sum()+coord_polar()+
    scale_x_continuous(limits=c(0,12), breaks=c(0,3,6,9,12))+
    annotate(geom="point",y=1,x=m2,colour="red")+
    annotate(geom="segment",x=bq2[["lwr"]],xend=bq2[["upr"]],y=1,yend=1,colour="red")