A late answer...
In a molecular approch to chemical reaction rates we asume that the molecules move randomly in solution and have a velocity randomly assigned according to the Boltzmann distribution and that the chance of being in one position in solution is the same as being at any other, i.e. the molecules are in thermal equilibrium at some temperature T.
The reaction rate for the reaction $\mathrm {A+B \rightarrow}$ can be found by considering the number of collisions per second and then this number reduced by the chance that a reaction occurs on collision. To arrive at a rate the chance that a randomly chosen combination of molecules of type A and B at time t will react in the next infinitesimal time interval $t+\delta t$ and has not reacted up to t has to be calculated. As a result the total number of collisions /second between $N_A$ molecules of type A and $N_B$ of type B is found to be $\sim N_AN_B\sigma_{AB}$, where $\sigma_{AB}$ is the cross sectional area for a collision, and this expression produces the rate when multiplied by the probability $p_{AB}$ of A and B reacting, i.e. rate $\sim N_AN_B\sigma_{AB}p_{AB}$. The probability of reaction $p_{AB}$ depends critically on the activation energy.
The important point is, however, that the rate is proportional to the product of the numbers of molecules of type A and B which in constant volume is the same as being proportional to the product of concentrations. That it also is proportional to the collision cross section and the chance of reacting is also physically realistic. The rate constant is proportional just to the probability of reaction and the cross section for collision.
Thus the origin of the product of concentrations in the rate expression is found in the chance that a molecule of type A will collide with that of type B in a small time interval $t, t +\delta t$. Thus if we ignore the fact that there is variation due to changes in experimental conditions, such as temperature, the electronics of the detectors used to measure the reaction etc. perfect kinetics would be measured, i.e $dA/dt$ or $ dB/dt$ vs. time would be noise free.
If we now consider the random nature of the collisions between molecules a slightly different approach has to be taken. If a small volume of the solution is considered and the rate calculated in that region, because the number of species entering and leaving in any given time period varies, the rate expression is also potentially subject to variation. This happens simply because of the random nature of the motion of molecules in a solution or gas. In this case the rate of change of a species X should be written in terms of its average (or mean) value represented as $\langle X(t)\rangle$;
$$ \frac{d\langle X(t)\rangle}{dt} = \langle f\times X(t)\rangle $$
where $\langle f\times X(t)\rangle$ represents some rate expression such as $k_1\langle X(t)\rangle + k_2\langle X(t)\rangle^2$ etc. except that the number of $X$ are represented as their average values. Clearly the $X$ are always positive values. These values have a variance (square of standard deviation) and so we would expect that different experiments give slightly different results because of this. (This is the approach taken by Gillespie, see 'Markov Processes' publ. Academic Press 1992)
In a bimolecular reaction such as $X + X \rightarrow $ with rate constant k the rate expression is
$$ \frac{d\langle X(t)\rangle}{dt} = k\langle \;X(t)(X(t)-1)\;\rangle/2 $$
The term $\langle X(t)(X(t)-1)\rangle /2$ arises becaus of the number of distinct X-X pairs in the reaction volume at time t. The question then arises as how to deal with this. Two approximations can be made; the first is to assume that X is large compared to $1$ and so $\langle X(t)(X(t)-1)\rangle \rightarrow \langle X^2(t)\rangle $, which seems to be a reasonably sensible approximation, the next step is more problematic and is $\langle X^2(t)\rangle \sim \langle X(t)\rangle^2$. This is a drastic step and assumes that the variance in the number of molecules is zero, i.e. that $X(t)=\langle X(t)\rangle$ which means that $X$ is a deterministic process, so that the randon nature of the numbers of $X$ entering into the reaction volume is always the same, which is not what was assumed at the outset and seems to fly in the face of what we know about how molecules move in solution, i.e. their random diffusional nature.
However, in any normal kinetic experiment the number of molecules we measure is vast, $\langle X(t)\rangle \sim 10^{18}$ would not be untypical, so that the fluctuations in the number of molecules is tiny and the approximation made above is valid.
$$ \frac{\sigma_X}{\langle X \rangle} \sim \frac{\sqrt{\langle X \rangle}}{\langle X \rangle} \sim \frac{1}{\sqrt{\langle X \rangle}}\sim 10^{-9}$$
Even if there are as few as 1000 molecules the realtive fluctuation is not too drastic but in experiments using small numbers of molecules, as in the case of many 'single molecule' experiments then clearly the rate equation approach needs to be modified.
Notes :
(A) To see that the variance is zero use the definition $\sigma^2_X = \langle X^2 \rangle -\langle X \rangle^2$ where $\sigma$ is the standard deviation)
(B) There are other situations where the rate equation approach breaks down in macroscopic systems and that is when they are sufficiently non-linear, one example is the Schlogl system of equations, which for certain values of concentrations and rate constants becomes bi-stable;
$$B_1+2X \rightleftharpoons 3X; \;\; B_2 \rightleftharpoons X$$
