1. A method of relational analysis of data, comprising: receiving a data set {{F (j) } j=1 m , {S (j) } j=1 m , {R (ij) } i,j=1 m } having a plurality of types of data objects which are interrelated to each other, the data objects having attributes of respective data objects {Θ (j) } j=1 m , homogeneous relations {Γ (j) } j=1 m between objects of the same type, and heterogeneous relations {γ (ij) } i,j=1 m between objects of different types, each type of object having latent classes; formulating a relational mixed membership model {Λ (j) } j=1 m representing data objects and respective relationships of data objects based on latent classifications, the model comprising an expectation having a latent indicator {C (j) } j=1 m for each respective object generated based on a multinomial distribution, with the class memberships of the respective data objects as parameters; automatically initializing the parameters of a likelihood function 
 {tilde over (Ω)}={{{tilde over (Λ)} (j) } j=1 m ,{{tilde over (Θ)} (j) } j=1 m ,{{tilde over (Γ)} (j) } j=1 m ,{{tilde over (γ)} (ij) } i,j=1 m , and computing the posterior 
 Pr ({ C (j) }|F (j) } j=1 m ,{S (j) } j=1 m ,{R (ij) } i,j=1 m ,{tilde over (Ω)}) using the Gibbs sampler, and maximizing the likelihood function by iteratively updating {tilde over (Ω)}=Ω over the range of i and j, until convergence is achieved, where the update is computed by at least one of: Method I Λ hp ( 1 ) = Pr ⁡ ( C hp ( 1 ) = 1 ❘ F , S , R , Ω ~ ) , Θ · g = ∑ p = 1 n 1 ⁢ ⁢ F · p ⁢ Pr ⁡ ( C gp ( 1 ) = 1 ❘ F , S , R , Ω ~ ) ∑ p = 1 n 1 ⁢ ⁢ Pr ⁡ ( C gp ( 1 ) = 1 ❘ F , S , R , Ω ~ ) Γ gh = ∑ p , q = 1 n 1 ⁢ ⁢ S pq ⁢ Pr ⁡ ( C gp ( 1 ) = 1 , C hq ( 1 ) = 1 ❘ F , S , R , Ω ~ ) ∑ p , q = 1 n 1 ⁢ ⁢ Pr ⁡ ( C gp ( 1 ) = 1 , C hq ( 1 ) = 1 ❘ F , S , R , Ω ~ ) , and Method II update {C (j) } j=1 m , to maximize the objective function log ⁢ ⌊ ( Ω ❘ Ψ ) = ∑ p = 1 n 1 ⁢ ⁢ ∑ g = 1 k 1 ⁢ ⁢ C gp ( 1 ) ⁢ log ⁢ ⁢ Λ gp ( 1 ) + ∑ q = 1 n 2 ⁢ ⁢ ∑ h = 1 k 2 ⁢ ⁢ C hq ( 2 ) ⁢ log ⁢ ⁢ Λ gp ( 2 ) - D ϕ 1 ⁡ ( F , Θ ⁢ ⁢ C ( 1 ) ) - D ϕ 2 ⁡ ( S , ( C ( 1 ) ) T ⁢ Γ ⁢ ⁢ C ( 1 ) ) - D ϕ 3 ⁡ ( R , C ( 1 ) ) T ⁢ Υ ⁢ ⁢ C ( 2 ) + τ , wherein Θ · g ( j ) = ∑ p : C gp ( j ) = 1 ⁢ ⁢ F · p ( j ) ∑ p = 1 n j ⁢ ⁢ C gp ( j ) , Γ gh ( j ) = ∑ p : C gp ( i ) = 1 , q : C hq ( j ) = 1 ⁢ ⁢ S pq ( j ) ∑ p = 1 n j ⁢ ⁢ C gp ( j ) ⁢ ∑ q = 1 n j ⁢ ⁢ C hq ( j ) , and Υ gh ( ij ) = ∑ p : C gp ( i ) = 1 , q : C hq ( j ) = 1 ⁢ ⁢ R pq ( ij ) ∑ p = 1 n i ⁢ ⁢ C gp ( i ) ⁢ ∑ q = 1 n j ⁢ ⁢ C hq ( j ) , to estimate unknown parameters of a joint probability distribution equation over the observations and latent indicators; and storing in a memory at least the maximized expectations associated with the model.