function [Sigma] = coupled_model_Bayesian(theta,xdata)

A = theta(1);
B = theta(2);
Pnum = xdata;

%fit energy to a polynomial
Sigma = A/2*(Pnum).^2 + B/4*(Pnum).^4; %Landau energy model