What is the probability distribution for a subsystem in canonical ensemble? Suppose we have a 3d Ising Model(NN interaction) in simple cubic lattice, if we define a subsystem of it to be a 2d plane of spins(for example all sites with z = L/2, L being the linear system size) and also define this subsystem energy just like a 2d Ising model(only 4 neighbor). 
What is the probability distribution of this subsystem states? 
For the 3d system itself this probability distribution is proportional to exp(-E/T) but I do not think for the 2d subsystem this holds.
 A: The total density operator for the 3D system is:
$$\rho = \frac{1}{\mathcal Z} e^{-\beta H}$$
where $H \equiv - J \sum_{\langle i,j\rangle} \sigma_i^z \sigma_j^z$ is the Ising Hamiltonian with nearest neighbor interactions. Now if you want to describe the density operator of the 2D subsystem (and hence the probability distribution), which I'll call $A$, you need to take a partial trace over all other degrees of freedom ($\bar A$), i.e. all spins outside that 2D plane. So the density operator of the subsystem is:
$$\rho_A \equiv \mathrm{Tr}_{\bar A} (\rho) = \frac{1}{\mathcal Z} \mathrm {Tr}_{\bar A}(e^{\beta J \sum_{\langle i,j\rangle} \sigma_i^z \sigma_j^z})$$
I don't think there is an analytic way to exactly evaluate such an expression; so you would probably need to resort to numerical methods for it. 
For reference, even calculating simple quantities like just the partition function of the Ising model is a hard problem (except in 1D and 2D where exact expressions exist). People usually resort to approximation methods such as mean field theory to calculate such simple quantities. So actually calculating stuff in the Ising model by hand is usually pretty hard.
