I am working on the measure of a Curie Weiss model over a lattice

$$ \mu[m] = \frac{e^{\beta N (hm + m^2)}}{\sum_z e^{\beta N (hz + z^2)}}\binom{N}{\frac{(1-m)N}{2}}2^{-N} $$

where $m\in\{-1,-1+2/N,\ldots,1-2/N,1\}$.

I would like to plot moments of the distribution for some largish, fixed, $N$ between $50$ and $200$ vs. $\beta$ between $1.1$ and $150$, in particular I would like to compute $ Cov(X,X^2) $ and $ V(X) $ (i.e. a measure of skewness of this distribution and a measure of variability).

I have the following code

g = 1;
\[Iota][m_, n_] := Binomial[n, n*(1 - m)/2]*2^(-n);
f[m_, h_, b_, g_, n_] := (h*m + g/2*m^2) + 
   1/(b*n)*Log[\[Iota][m, n]];
\[Mu][m_, h_, b_, g_, n_] := 
  Exp[-b*n*f[m, h, b, g, n]]/
   Sum[Exp[-b*n*f[x, h, b, g, n]], {x, -1 + 2/n, 1 - 2/n, 2/n}];
moment[h_, x_, b_, g_, n_] := 
  moment[h, x, b, g, n] = 
   Sum[m^x*\[Mu][m, h, b, g, n], {m, -1 + 2/n, 1 - 2/n, 2/n}];
var[h_, b_, g_, n_] := moment[h, 2, b, g, n] - moment[h, 1, b, g, n]^2;
cov[h_, b_, g_, n_] := 
  moment[h, 3, b, g, n] - moment[h, 1, b, g, n]*moment[h, 2, b, g, n];

but looking at various plot or values of

$$ Cov(X,X^2)+2V(X) $$ either the plot has sometimes erratic behavior (changes in behavior for one specific value of N) and trying to compute the value throws an error

enter image description here

it seems the exponential in this computation are too large.

1) Plotting throws no error, can I trust the plot in anyway?

2) How can I work with this in mathematica?

EDIT 2: Here's an example with parameter numbers:

Table[-b*(cov[h, b, g, n] + 2 var[h, b, g, n])], {b,   1.1, 100}]

EDIT 1: One option I considered is to use the following measure

$$ \mu[m] = \frac{e^{\beta N (hm + m^2)-\beta N \bar{f}}}{\sum_z e^{\beta N (hz + z^2)-\beta N \bar{f}}}\binom{N}{\frac{(1-m)N}{2}}2^{-N} $$

with an appropriate $\bar{f}$, which gives me the same probability but should be more stable numerically. Doing some experiments it seems Mathematica stopped throwing precision error.

But I would like some reassurances: can I rely on mathematica simplifying the exponent first and carry out the calculation later?

  • 1
    $\begingroup$ Can you please give example parameters and a concrete call to your functions that results in this warning, so that the reader does not need to speculate? $\endgroup$
    – Roman
    May 15, 2019 at 13:37


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.