Ich möchte Stichproben mit präzisen Momenten aus einer Dichtefunktion (PDF) ziehen. Es geht also um die Standardfehler von Kurtosis, Schiefe, Standardabschweichung und Mittel, und zwar insbesondere in dieser Priorität.
Dazu habe ich - nach den wie stelle ich eine gute Frage Richtlinien des Forums - natürlich auch ein kleines Programm geschrieben.
Code: Alles auswählen
#################################
# Standardfehler der Momenten (EW, STD, Schiefe, Kurtosis) von Stichproben (SH)
##################################
# Some Setup
skewness_man= function(x) {
n= length(x)
n/(n-1)/(n-2) * sum( ((x- mean(x))/sd(x))^3)
}
kurtosis_man= function(x) {
n= length(x)
n*(n+1)/(n-1)/(n-2)/(n-3) * sum( ((x- mean(x))/sd(x))^4)
}
# Function to construct a measurement population in data
pop2data= function(x, dd, m) {
medpos= round(length(dd)/2)
ffloor= floor(dd *m)
frest= dd* m - ffloor
frestcum= cumsum(frest)
flow= rep(0, length(x))
for (i in 1:medpos) {
if (frestcum[i]>= 0.5) {
flow[i]= 1
frestcum[(i+1):medpos]= frestcum[(i+1):medpos] -1
}
else {
flow[i]= 0
}
}
frestcumrev= rev(cumsum(rev(frest)))
for (i in length(x):medpos) {
if (frestcumrev[i]>= 0.5) {
flow[i]= 1
frestcumrev[medpos:(i-1)]= frestcumrev[medpos:(i-1)] -1
}
else {
flow[i]= 0
}
}
f= ffloor + flow
rep(x, f)
}
# pop2data Multiplier
m= 1000
# Sample Size
n= 10000
# Number of Samplings done
# o= 100
o= 1000
# Random Variable
u= seq(-100, 100, by= 0.01)
length(u)
##################################
# Main
# Construction mu
cmu= 0
# Construction deviation
csigma= 1
# Asymmetry paramater
a= 5
# probability density function (PDF)
f= function(u) { 1/(2 *sqrt(2*pi)) *exp(-0.5*((u-cmu)/csigma)^2)* (2/(1 +exp(-a*((u-cmu)/csigma)))) }
norm= integrate(f, -Inf, Inf)$value
# density data
dd= f(u)
# direct draw of the population with PDF f
ud= pop2data(u, dd, m)
length(ud)
# calculating its moments
u_mean= mean(ud)
u_s= sd(ud)
u_skew= skewness_man(ud)
u_kurt= kurtosis_man(ud)
# drawing sample from u with the density dd
sd= sample(u, n, replace= T, prob= dd)
# calculating its moments
s_mean= mean(sd)
s_s= sd(sd)
s_skew= skewness_man(sd)
s_kurt= kurtosis_man(sd)
# True moments of the poulation and random moments of the random sample
u_mean
u_s
u_skew
u_kurt
s_mean
s_s
s_skew
s_kurt
# Calculating the moments' of o multiple random samples
results= data.frame()
for (i in 1:o) {
# drawing sample from u with the density dd
sd= sample(u, n, replace= T, prob= dd)
# calculating its moments
s_mean= mean(sd)
s_s= sd(sd)
s_skew= skewness_man(sd)
s_kurt= kurtosis_man(sd)
results= rbind(results, t(c(s_mean, s_s, s_skew, s_kurt)) )
}
colnames(results)= c("mean", "s", "skew", "kurt")
# Moments' numerical expected value of o multiple random samples
ev_mean= mean(results$mean)
ev_s= mean(results$s)
ev_skew= mean(results$skew)
ev_kurt= mean(results$kurt)
ev_mean
ev_s
ev_skew
ev_kurt
# Moments' numerical standard errors of o multiple random samples
se_mean= sd(results$mean)
se_s= sd(results$s)
se_skew= sd(results$skew)
se_kurt= sd(results$kurt)
se_mean
se_s
se_skew
se_kurt
Fragen
1. Gibt es in R andere Stichprobenfunktionen?
2. Könnte ich das stratified Bootstrapping aus boot für meinen Zweck missbrauchen?
3. Was ist der Quellcode von sample.int ?
4. Wie sampled man manuell mit erhöhter Präzision?
Bereits von mir durchgeführte Lösungsversuche
1. "reject/pursuit" Sampling. Ich ziehe mit sample() solange eine Stichprobe, berechne aus dem Vergleich mit pop2data-Momenten ihre Momentenfehler und verwerfe sie, bis eine Ziehung meine Momentenfehlererwartungen erfüllt. Was sich zunächst toll anhört, ist aber nur ein Sturm im Wasserglass. Denn bei meinen idealen Fehlerwartungen, muss ich ungefähr 1000 bis 10000 mal ziehen, bis ich eine passende Stichprobe erhalte. Ich hab das dann irgendwann abgebrochen.
2. Ich ziehe einfach aus meiner Messpopulation aus pop2data. Das Problem daran ist, dass die Messpopulation nur eine mögliche Messpopulation von vielen möglichen der gegebenen PDF ist. Immerhin konnte ich auf diesem Weg feststellen, dass symmetrisches Ziehen (q-> -q) die Standardabweichung der Schiefe verrringert (aber die Standardfabweichung der Kurtosis steigen lässt.)
Das ist derzeit alles, was ich zur Lösung beitragen kann.
Insbesondere würde mir auch ein Verschieben der Präzision auf die hohen Momente (Kurtosis, Schiefe) sehr weiterhelfen, selbst wenn die Momente Standardabweichung und Erwartungswert dadurch grottenschlecht würden.