dist = MultinormalDistribution[IdentityMatrix]
pts = RandomVariate[dist, 500];
Graphics[Point[pts], Axes -> True]
You may chose a different variance var by e.g.:
dist = MultinormalDistribution[ var IdentityMatrix]
I don't know why you change notation and produce what seems to be a more complicated than necessary version of $g(x)$. If you just add assumptions to the integration, then you'll get an answer with your original code. I would think that defining $g(x)$ as
g[x_]:= a (b/(1 + Exp[-μ Ps x/r^α + ϕ]) - 1)
would be more straightforward. Here you're using the ...
$P_v$ is a mixture of a continuous random variable and a discrete random variable (with a probability mass at 0). Mathematica's MixtureDistribution requires random variables in the mixture distribution to be all discrete or all continuous (and of the same dimension). The point is that the determination of distributions and associated summaries usually can'...