In R how can I fit multivariative distribution to data and sample from it? I collected data on post parcells delivery, each object is a combination of 3 variables:


*

*Weight (continuous >0)  

*Destination city (categorial - factor of hundreds)

*Delivery type (categorial - factor of 3)


I want to fit a multivariative distribution to my data. Later I need to sample data from it for simulation purposes. I use R.
If I undestand things correctly in total continous case one might try gaussian mix or copula approach (using mixtools or copula R packages)
But my the only continous variable Weight is NOT normally distributed and two others are categorial.
How to deal with this kind of problem?
Will appresiate any help: general directions in statistics theory or R packages.
UPDATE
I have approximately 1 million data rows of each Delivery Type.
Also Kernel density estimation of Weight distribution was made for two major types. So shape of distribution is heavy-tail and slightly depends on Type, for ex. upper chart for Type 1 is a bit steeper:

UPDATE 2
I found R package np very useful for nonparametric conditional/joint distribution estimation.
 A: I think it depends on how much data you have.  For simplicity let's call weight, city and type $W, C$ and $T$ respectively.  If you have a good amount of data you could fully estimate the dependence structure between these three variables by first estimating the joint probability mass function of $C$ and $T$ (that is, $P(T = t, C = c)$ for every combination of $t$ and $c$), and then the conditional distribution of $W$ for each possible value of $C$ and $T$.  This could be problematic since $C$ can take on hundreds of values and the number of parameters you'd need to estimate is proportional to the number of cities.
The other extreme would be the naive approach of assuming $W, C$ and $T$ are independent, in which case you would estimate all the marginal distributions separately and when you do your simulation you simply draw from each of these marginals individually and combine the results.
My suggestion would be somewhere in between and would again depend on how much data you have as well as which variables you can reasonably expect to have an influence on the others.  I would suspect that $W$ depends more on $T$ than on $C$, so try treating $W$ and $C$ independent but model the distribution of $W$ conditional of the three values of $T$.  You may be able to get away with treating $T$ and $C$ as independent for simplicity's sake.  For the distribution of $W$ some nonparametric approaches would be to use the empirical distribution function (when you sample you just draw from the original distribution with replacement) or a kernel density estimate.  These are nice approaches and quite flexible but require that you "remember" the original data.  To reduce the complexity of your model there's essentially an unlimited number of parametric approaches you can take to density estimation.  Since $W$ is nonnegative you could consider treating it as exponential, log-normal, Pareto, and so on.  It depends on how your data are structured.  Also, does the shape of the distribution depend on $T$?  These are just a few ideas and it's hard to say more without seeing the actual data.
Regarding your question about how to test for dependence between a categorical and continuous variable, this could be approached using an analysis of variance, but unfortunately the ANOVA test itself is unlikely to tell you much about what's going on except in the rare case that you fail to reject the null.  Goodness of fit tests could also be used to test for equality of distributions, but these would suffer from the same shortcomings of the ANOVA test, namely that the null is generally false anyways (but not necessarily in important ways).  It's always a good idea to look at the data, so consider some parallel boxplots to try to visualize how $W$ might depend on $T$, or again check out some kernel density estimates (which I see from your update you've already done).
