Prawidłowe jest porównanie kilku podejść, ale nie w celu wybrania tego, które faworyzuje nasze pragnienia / przekonania.
Moja odpowiedź na twoje pytanie brzmi: możliwe jest, że dwie dystrybucje pokrywają się, gdy mają różne środki, co wydaje się być twoim przypadkiem (ale musielibyśmy zobaczyć twoje dane i kontekst, aby podać bardziej precyzyjną odpowiedź).
Zilustruję to za pomocą kilku podejść do porównania normalnych środków .
1. testt
Rozważ dwie symulowane próbki o rozmiarze z i , a następnie wartość wynosi około jak w twoim przypadku (patrz kod R poniżej).70N(10,1)N(12,1)t10
rm(list=ls())
# Simulated data
dat1 = rnorm(70,10,1)
dat2 = rnorm(70,12,1)
set.seed(77)
# Smoothed densities
plot(density(dat1),ylim=c(0,0.5),xlim=c(6,16))
points(density(dat2),type="l",col="red")
# Normality tests
shapiro.test(dat1)
shapiro.test(dat2)
# t test
t.test(dat1,dat2)
Jednak gęstości wykazują znaczne nakładanie się. Pamiętaj jednak, że testujesz hipotezę o środkach, które w tym przypadku są wyraźnie różne, ale ze względu na wartość nakładają się na siebie gęstości.σ
2. Prawdopodobieństwo profiluμ
Definicja prawdopodobieństwa i prawdopodobieństwa profilu znajduje się w punktach 1 i 2 .
W tym przypadku prawdopodobieństwo profilu próbki o wielkości i średniej próbki wynosi po prostu .μnx¯Rp(μ)=exp[−n(x¯−μ)2]
W przypadku danych symulowanych można je obliczyć w R w następujący sposób
# Profile likelihood of mu
Rp1 = function(mu){
n = length(dat1)
md = mean(dat1)
return( exp(-n*(md-mu)^2) )
}
Rp2 = function(mu){
n = length(dat2)
md = mean(dat2)
return( exp(-n*(md-mu)^2) )
}
vec=seq(9.5,12.5,0.001)
rvec1 = lapply(vec,Rp1)
rvec2 = lapply(vec,Rp2)
# Plot of the profile likelihood of mu1 and mu2
plot(vec,rvec1,type="l")
points(vec,rvec2,type="l",col="red")
Jak widać przedziały prawdopodobieństwa i nie pokrywają się na żadnym rozsądnym poziomie.μ1μ2
3. tylnego używając Jeffreys przedμ
Rozważmy Jeffreys przed o(μ,σ)
π(μ,σ)∝1σ2
Tylne dla każdego zestawu danych można obliczyć w następujący sposóbμ
# Posterior of mu
library(mcmc)
lp1 = function(par){
n=length(dat1)
if(par[2]>0) return(sum(log(dnorm((dat1-par[1])/par[2])))- (n+2)*log(par[2]))
else return(-Inf)
}
lp2 = function(par){
n=length(dat2)
if(par[2]>0) return(sum(log(dnorm((dat2-par[1])/par[2])))- (n+2)*log(par[2]))
else return(-Inf)
}
NMH = 35000
mup1 = metrop(lp1, scale = 0.25, initial = c(10,1), nbatch = NMH)$batch[,1][seq(5000,NMH,25)]
mup2 = metrop(lp2, scale = 0.25, initial = c(12,1), nbatch = NMH)$batch[,1][seq(5000,NMH,25)]
# Smoothed posterior densities
plot(density(mup1),ylim=c(0,4),xlim=c(9,13))
points(density(mup2),type="l",col="red")
Ponownie, przedziały wiarygodności środków nie pokrywają się na żadnym rozsądnym poziomie.
Podsumowując, można zobaczyć, w jaki sposób wszystkie te podejścia wskazują na znaczną różnicę środków (co jest głównym przedmiotem zainteresowania), pomimo nakładania się rozkładów.
⋆ podejście porównawcze
Sądząc z twoich obaw związanych z nakładaniem się gęstości, kolejną interesującą wartością może być , prawdopodobieństwo, że pierwsza zmienna losowa jest mniejsza niż druga zmienna. Ilość tę można oszacować nieparametrycznie, jak w tej odpowiedzi . Pamiętaj, że nie ma tu żadnych założeń dystrybucyjnych. W przypadku danych symulowanych estymator ten wynosi , co pokazuje pewne nakładanie się w tym sensie, podczas gdy średnie są znacznie różne. Proszę spojrzeć na kod R pokazany poniżej.0,8823825P(X<Y)0.8823825
# Optimal bandwidth
h = function(x){
n = length(x)
return((4*sqrt(var(x))^5/(3*n))^(1/5))
}
# Kernel estimators of the density and the distribution
kg = function(x,data){
hb = h(data)
k = r = length(x)
for(i in 1:k) r[i] = mean(dnorm((x[i]-data)/hb))/hb
return(r )
}
KG = function(x,data){
hb = h(data)
k = r = length(x)
for(i in 1:k) r[i] = mean(pnorm((x[i]-data)/hb))
return(r )
}
# Baklizi and Eidous (2006) estimator
nonpest = function(dat1B,dat2B){
return( as.numeric(integrate(function(x) KG(x,dat1B)*kg(x,dat2B),-Inf,Inf)$value))
}
nonpest(dat1,dat2)
Mam nadzieję, że to pomoże.