Modeling complex systems, like neural networks, simple liquids or flocks of birds, often works in reverse to textbook approaches: given data for which averages and correlations are known, we try to find the parameters of a given model consistent with it. In general, no exact calculation directly from the model is available and we are left with expensive numerical approaches. A particular situation is that of a perturbed Gaussian model with polynomial corrections for continuous degrees of freedom. Indeed perturbation expansions for this case have been implemented in the last 60 years. However, there are models for which the exactly solvable part is non-Gaussian, such as independent Ising spins in a field, or an ideal gas of particles. We implement a diagrammatic perturbative scheme in weak correlations around a non-Gaussian yet solvable probability weight. This applies in particular to spin models (Ising, Potts, Heisenberg) with weak couplings, or to a simple liquid with a weak interaction potential. Our method casts systems with discrete degrees of freedom and those with continuous ones within the same theoretical framework. When the core theory is Gaussian it reduces to the well-known Feynman diagrammatics.