# Exponential distribution and numerical errors?

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

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:

h=-0.3;
g=1;
n=100;
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?

• 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? – Roman May 15 '19 at 13:37