Si je comprends bien votre question, c'est assez facile. Il vous suffit de décider de la distribution que vous souhaitez que vos erreurs aient et d'utiliser la fonction de génération aléatoire correspondante.
Il existe un certain nombre de distributions asymétriques, vous devez donc déterminer celle que vous aimez. De plus, la plupart des distributions asymétriques (par exemple, log normal, chi carré, Gamma, Weibull, etc.) sont asymétriques à droite, donc quelques adaptations mineures seraient nécessaires (par exemple, multiplier par ). - 1
Voici un exemple de modification de votre code:
set.seed(5840) # this makes the example exactly reproducible
N <- 100
x <- rnorm(N)
beta <- 0.4
errors <- rlnorm(N, meanlog=0, sdlog=1)
errors <- -1*errors # this makes them left skewed
errors <- errors - 1 # this centers the error distribution on 0
y <- 1 + x*beta + errors
Je dois noter à ce stade que la régression ne fait aucune hypothèse sur les distributions de ou , seulement sur les erreurs, (voir ici: Que faire si les résidus sont normalement distribués, mais y ne l'est pas? ). C'était donc l'objet de ma réponse ci-dessus. XOuiε
Mise à jour: Voici une version asymétrique avec les erreurs distribuées comme Weibull:
set.seed(5840) # this makes the example exactly reproducible
N <- 100
x <- rnorm(N)
beta <- 0.4
errors <- rweibull(N, shape=1.5, scale=1)
# errors <- -1*errors # this makes them left skewed
errors <- errors - factorial(1/1.5) # this centers the error distribution on 0
y <- 1 + x*beta + errors
Les données Weibull sont déjà faussées à droite, nous n'avons donc pas besoin de changer de direction (c'est-à-dire que nous supprimons la -1*errors
pièce). De plus, à partir de la page Wikipedia pour la distribution de Weibull, nous voyons que la moyenne d'un Weibull devrait être:. Nous voulons soustraire cette valeur de chacune des erreurs afin que la distribution des erreurs résultante soit centrée sur . Cela permet à la partie structurelle (c'est-à-dire ) de votre code de refléter avec précision la partie structurelle du processus de génération de données. E[ W] = ( 1 / s h a p e ) !01 + x*beta
La distribution exgaussienne est la somme d'une normale et d'une exponentielle. Il y a une fonction ? RexGAUS dans le package gamlss.dist pour les générer. Je n'ai pas ce package, mais vous devriez pouvoir adapter mon code ci-dessus sans trop de difficulté. Vous pouvez également générer une variable normale aléatoire (via rnorm()
) et une exponentielle (via rexp()
) et les additionner assez facilement. N'oubliez pas de soustraire la moyenne de la population, , de chaque erreur avant d'ajouter les erreurs à la partie structurelle du processus de génération de données. (Attention cependant à ne pas soustraire la moyenne de l' échantillon !) μ + 1 / λmean(errors)
Quelques commentaires finaux sans rapport: Votre exemple de code dans la question est quelque peu confus (ce qui signifie aucune infraction). Parce que rnorm(N)
génère des données avec mean=0
et sd=1
par défaut, 0.4*rnorm(N)
générera rnorm(N, mean=0, sd=0.4)
. Votre code (et éventuellement votre réflexion) sera beaucoup plus clair si vous utilisez cette dernière formulation. De plus, votre code pour beta
semble confus. Nous pensons généralement à laβdans un modèle de type régression en tant que paramètre, pas une variable aléatoire. Autrement dit, c'est une constante inconnue qui régit le comportement du processus de génération de données, mais la nature stochastique du processus est encapsulée par les erreurs. Ce n'est pas ainsi que nous pensons à ce sujet lorsque nous travaillons avec des modèles à plusieurs niveaux, et votre code semble à mi-chemin entre un modèle de régression standard et le code d'un modèle de régression à plusieurs niveaux. La spécification de vos bêtas séparément est une bonne idée pour maintenir la clarté conceptuelle du code, mais pour un modèle de régression standard, vous attribuez simplement un numéro unique à chaque version bêta (par exemple, beta0 <- 1; beta1 <- .04
).