We report a systematic study of the collective effect of thermal plumes on the probability density function (p.d.f.) $P(\delta T)$ of temperature fluctuations $\delta T(t)$ in turbulent Rayleigh–Bénard convection. By decomposing $\delta T(t)$ into four basic fluctuation modes associated with single and multiple warm and cold plumes and a turbulent background, we derive an analytic form of $P(\delta T)$ based on the convolutions of the five independent modes. To test the derived form of $P(\delta T)$ in the multiple-plume regions, where the thermal plumes are heavily populated, we conduct time series measurements of temperature fluctuations in two convection cells; one is a vertical thin disk and the other is an upright cylinder of aspect ratio unity. For a given normalized position in most regions of the convection cell, all of the measured p.d.f.s $P(\delta T)$ for different Rayleigh numbers fall onto a single master curve, once $\delta T$ is normalized by its root-mean-square (r.m.s.) value $\sigma _T$. It is found that the measured $P(\delta T/\sigma _T)$ at different locations along the symmetric horizontal and vertical axes of the convection cells can all be well described by the derived form of $P(\delta T/\sigma _T)$. The fitted values of the parameters associated with the number of plumes in multiple plume clusters and their relative strengths and degrees of intermittency are closely linked to the spatial distribution of thermal plumes and local dynamics of the large-scale circulation in a closed convection cell. Our work thus provides a unified theoretical approach for understanding scalar p.d.f.s in a turbulent field, which is very useful not only for the present study but also for the study of many turbulent mixing problems of practical interest.