The study of functional networks in the brain is essential in order to gain a better insight into its diverse set of operations and to characterise the associated normal and abnormal behaviours. Present methods of analysing fMRI data to obtain functional connectivity are largely limited to approaches such as correlation, regression and independent component analysis, which give simple point estimates. By contrast, we propose a stochastic linear model in a Bayesian setting and employ Markov Chain Monte Carlo methods to approximate posterior distributions of full connectivity and covariance matrices. Through the use of a Bayesian probabilistic framework, distributional estimates of the linkage strengths are obtained as opposed to point estimates, and the uncertainty of the existence of such links is accounted for. We decompose the connectivity matrix as the Hadamard product of binary indicators and real-valued variables, and formulate an efficient joint-sampling scheme to infer them. The well-characterised somato-motor network is examined in a self-paced, right-handed finger opposition task based experiment, while nodes from the visual network are used for contrast during the same experiment. Unlike for the visual network, significant changes in connectivity are found in the motor network during the task. Our work provides a distributional metric for functional connectivity along with causality information, and contributes to the collection of network level descriptors of brain functions.