Estimating Tsallis q

From InterSciWiki

THE PURPOSE HERE IS TO FIND A MORE DIRECT Maximum Likelihood Estimation for q-exponential distributions, following the procedures that have been outlined to us by Tsallis and Borges for fitting a straight line (by MLE!) to x and $e_q^x$ . Help will be greatly appreciated.

See embedded distributions approach.

http://en.wikipedia.org/wiki/Tsallis_entropy

Sandbox has the older page

1 q exponential and q logarithm
2 Probability distribution and Cumulative probability distribution (CDF)
3 Comparison to the Shalizi MLE by Pareto II
4 Possibility of a more direct solution
5 Deriving the MLE
6 Testing the MLE 1
7 Testing the MLE 2 ???
8 Comparison to the Shalizi MLE by Pareto II
- 8.1 Another set of questions
9 References
10 See also
11 Note 1
12 Note 2: Eliminated - Curve fitting: estimating q with the q-logarithmic function

[edit] q exponential and q logarithm

The solution to $\frac{dy}{dx} = y^q \ \mathrm{(equation\ 1)}$ , y(0)=1, is the q-exponential function (Tsallis 2004:5-6). A constant such as $λ$ may be inserted before x and will carry over before x in equations 1 and 2)

$e_q^x \equiv [1+(1-q)x]^{1/(1-q)} \ \mathrm{(equation\ 2)}$ if (1 + (1 - q)x) > 0); otherwise $e_q^x=0$ , where $(e_1^x = e^x)$ ,

The inverse of the q-exponential is the q-logarithmic function (Tsallis 2004:6):

$ln_q\ x \equiv \frac{ x^{1-q} -1}{1-q} \ \mathrm{(equation\ 3)}$ where (ln₁ x = ln x)

$ln_q\ (e_q^x) \equiv x$

With an optional constant $κ$

$ln_q\ x\kappa \equiv \frac{ (x\kappa)^{1-q} -1}{1-q} \ \mathrm{(equation\ 3')}$ where $ln_1 \ x \ \kappa = ln \ x \ \kappa = ln \ x + ln \ \kappa$

[edit] Probability distribution and Cumulative probability distribution (CDF)

The probability distribution for $e_q^x$ (Tsallis 2004:6-8) normalized with $\int_{0}^{\infin} dx\ p_{e_q^x}(x) = 1$ is

$p_{e_q^x}(x) \equiv (2-q)e_q^{-x} (x \ge 0) = (2-q)[1-(1-q)x]^{1/(1-q)} \ \mathrm{(equation\ 4)}$

The cumulative probability distribution for the q-exponential function in (1) (Tsallis 2004:8) is also a q-exponential function

$P(X) \equiv \int_{X}^{\infin} dx\ p_{e_q^x}(x) = e_{q_M}^{-X/q_M} = [1 - \frac{(1-q)}{(2-q)}x]^\frac{(2-q)}{(1-q)} \ \mathrm{(equation\ 5)}$ , with

$q_M \equiv 1/(2-q)$ , where for $q \ge 2\ \ p(x)$ is not normalizable. Thus, q cannot be fitted for a range at or above 2. (Note that

κ

does not enter into this equation).

[edit] Comparison to the Shalizi MLE by Pareto II

See Cosma Shalizi, 2007 Maximum Likelihood Estimation for q-Exponential (Tsallis) Distributions. http://www.cscs.umich.edu/~crshalizi/research/tsallis-MLE These are classically known as Pareto II distributions. (Examples). Rather than using Shalizi's method, however, we procede to estimate q and $κ$ directly, as follows (the two MLE approaches will be compared).

As related to Tsallis's exposition, above, Shalizi took as his cumulative distribution

$P_{q,\lambda}(X) \equiv e_{q_M}^{-X/q_M} = [1 - (1-q_M)x/\lambda]^{1/(1-q_M)}, \lambda=1/\kappa \ \mathrm{(equation\ A\ =\ 1\ in\ Shalizi)}$ ,

Substituting to obtain Pareto II parameters, $,θ,σ$ , we may take his equations (8) and (10) for the MLE estimation of $\hat{\theta}, \hat{\sigma}$ in terms of $\hat{q}, \hat{\kappa}$

$q_M \equiv \frac{1}{(2-q)} = \frac{\theta+1}{\theta}\ \mathrm{hence} \ q_M = \frac{\theta+2}{\theta+1} \ne \frac{\hat{\theta}+1}{\hat{\theta}}\ \mathrm{as\ in\ Shalizi}\ \mathrm{and \ thus}$

$\hat{q} =\hat{q_M} = \frac{\hat{\theta}+2}{\hat{\theta}+1}$

$\hat{\kappa} = \hat{\theta}/\hat{\sigma}$

[edit] Possibility of a more direct solution

Is it possible that $ln_q\ x \equiv \frac{ (\kappa x)^{1-q} -1}{1-q} \approx \kappa x$ is a solution for q?

Or that $\sum_{i=1}^{k} ln_q\ x_i \equiv \sum_{i=1}^{k}\frac{ (\kappa x_i)^{1-q} -1}{1-q} \approx \kappa x_i$ is a solution for q?

Let $\eta \equiv -\frac{(2-q)}{(1-q)}$ , so (equation 5) becomes (following Shalizi 2007:1, and his eq. 3)

$P(X\ge x) \equiv [1 + x /\eta ]^{-\eta} \ \mathrm{(equation\ 6)}$ , with the probability density

$\boldsymbol{P}(X\ge x)=\int_a^b f(x)\,dx$ , where

f (x)

is the density with a total integral of 1.

Following Shalizi 2007:1, and his eq. 4), the function whose integral is $P(X\ge x)$ is

$p_\eta(x) \equiv (1 + x/\eta)^{-\eta-1} \ \mathrm{(equation\ 7)}$

THIS FIGURE IS FOR THE CUMULATIVE DISTRIBUTION

Figure from Thurner and Tsallis 2005:200 click for reference

$P(X_k) \equiv Z(k|x_i\ge k)= \sum_{i=1}^{k} ln_q\ x_i = y(k)$ (note that $κ$ does not enter into this equation)

[edit] Deriving the MLE

For maximum likelihood we need the correct probability density function. Following Shalizi (2007:1-2, and his eq. 5), under the q-exponential model with parameter $η$ , the log-probability density of a sequence of independent identically-distributed samples is

$log \ p_\eta(x_1^n) = - (\eta +1)\sum_{i=1}^{n} log\ (1 + x_i/\eta) = \mathfrak{k}(\eta) \ \mathrm{(equation\ 8)}$

Then (following Shalizi 2007:2, and his eq. 7): (see final Note 1 below)

$\delta \mathfrak{k}/\delta \eta = n/\eta - \sum_{i=1}^{n} log \ (1 + x_i/\eta) \ \mathrm{(equation\ 9), and \ again}$ (following Shalizi 2007:2, and his eq. 8):

To find the MLE (maximum $\mathfrak{k}(\eta)$ , take the first derivative of the log-likelihood $\ \mathfrak{k}(\eta)$ with respect to the parameter $η$ and set it equal to zero to solve for the maximum.

$\hat{\eta} = n \left [ {\sum_{i=1}^n} log(x_i/x_{min}) \right ]^{-1} \ \mathrm{(equation\ 10)}$ , where in the simplest case

x m i n = η

, but we may also take the value of

x m i n

that maximizes fit according to the Kolmogorov-Smirnov (K-S) bootstrap method.

The solution for $\hat{q} \mathrm{\ from} \ \eta \equiv -(2-q)/(1-q)$ , rearranging terms, is simply

$0 < \hat{q} = \frac{\hat{\eta}+2}{\hat{\eta}+1} < 2$ .

Let us denote by $\hat{q_W} \mathrm{\ and} \ \hat{q_S} = 1 + \frac{1}{\theta}_S$ the MLE estimates by this method (W) and by Shalizi's (S) method. Because $\hat{\eta}\ \mathrm{ and} \ \hat{\theta}$ are both derived from the same equation (W: equation 10 above) if we set $x m i n = η$ , comparison of parameter transformations show:

$\hat{q_W} = \frac{\hat{\eta}+2}{\hat{\eta}+1} = \frac{\hat{\theta}_S+2}{\hat{\theta}_S+1} = \frac{1+2\hat{q_S}-1}{1+\hat{q_S}-1} = 2 - \frac{1}{\hat{q_S}}$ .

Here, both $\hat{q_W}\ \mathrm{ and} \ \hat{q_S}$ are upper bounded by 2, but the problem is with $x m i n$ . So it seems we might need a different R program to compute the MLE than Shalizi's. It will be easier, however, to modify Clauset's qpva.m program and check the validity of this new MLE. This new formula should be

[edit] Testing the MLE 1

Cosma's 2007 equation 4, page 1, used in Aarom's qpva.m RESULTS WILL DIFFER EACH TIME BECAUSE WIDTHS OF THE INTERVALS BETWEEN x VALUES WILL VARY RANDOMLY

Not Normalized to a probability because theta > sigma

xmin = 0;
theta = 5;
sigma = 1.3;
n = 10^2;
x = (1+(1-rand(n,1)/sigma)).^(-theta - 1);
figure(1);
loglog(sort(x),(n:-1:1)./n);

Normalize to a probability

xmin = 1;
theta = 2.1;
sigma = 1.1;
n = 10^2;
x = (theta/sigma)*(1+(1-rand(n,1)/sigma)).^(-theta - 1);
figure(1);
loglog(sort(x),(n:-1:1)./n);
q=1+1/theta
kappa=theta/sigma

(note inversion of kappa by Cosma in his equation 1 and 2)

Now try $ln_q\ x\kappa \equiv \frac{ (x\kappa)^{1-q} -1}{1-q} \ \mathrm{(equation\ 3)\ where \ } q = 1 + 1/\theta, \kappa = \theta/\sigma \ne \sigma/\theta$

[edit] x not cumulative

n = 10^2;
q=1.5;
kappa=1.9
sigma=1.9
x = [(kappa*(1-rand(n,1)/sigma)).^(1-q)-1]/(1-q);
figure(1);
loglog(sort(x),(n:-1:1)./n);

n = 10^2;
q=1.5;
kappa=1.9
sigma=1.9
x = [(kappa*(1-rand(n,1)/sigma)).^(1-q)-1]/(1-q);
figure(1);
plot(sort(x),(n:-1:1)./n);

[edit] straight line?

Cumulative

q=1.47623;
kappa=1.9;
sigma=1.909;
n = 10^2;
x = sort(1-rand(n,1));
for i=1: n
 xx(i)=sum(x(i));
end
%xx=cumsum(x)
lnqx = [(kappa*xx/sigma).^(1-q)-1]/(1-q);
figure(1);
plot(1-x,lnqx);

Not cumulative

q=1.47623;
kappa=1
sigma=1.909
n = 10^2;
x = sort(1-rand(n,1));
lnqx = [(kappa*x/sigma).^(1-q)-1]/(1-q);
figure(1);
plot(x,lnqx);

lnqx=((x(n,1)/sigma)^(1-theta)-1)/(1-theta)
figure(2);
loglog(sort(lnqx),(n:-1:1)./n);

[edit] try new MLE

Then try (see if theta is returned) $\hat{\eta} = n \left [ {\sum_{i=1}^n} ln(x_i/x_{min}) \right ]^{-1} \ \mathrm{(equation\ 10)}$ ,

Or from Cosma's Equation 4 should this be

xmin = 0;
theta = 2;
sigma = 3;
n = 10^2;
x = (xmin+(1-rand(n,1)/sigma)).^(-theta-1);
figure(1);
loglog(sort(x),(n:-1:1)./n);

q=0.8;
n = 10^2;
k=2
sigma=3
x = [(k*(1-rand(n,1)/sigma)).^(1-q)-1]/(1-q);
figure(1);
loglog(sort(x),(n:-1:1)./n);

[edit] Testing the MLE 2 ???

Hi Laurent,

The best way to test whether your MLE is accurate is to try it out on synthetic data, for which you know the correct q-exponential parameters. The qpva.m code includes a way of generating synthetic data drawn from the Pareto II, formulated as by Cosma's article, so potentially you could use that code to generate synthetic data for your MLE -- just remember that you'll need Eq 2 in Cosma's article to convert the Pareto II parameters into q-exponential parameters.

Sincerely, Aaron

On Jan 15, 2008, at 2:51 PM, Laurent Tambayong wrote:

> > Hello Aaron,
> > 
> > We need to test our derivation.  http://intersci.ss.uci.edu/wiki/index.php/Estimating_Tsallis_q
> > 
> > Is that possible to adjust QPVA to be based on equation 10 (and its corresponding q underneath) ?  Currently that code is based on eq 11, 12, and 13 from Cosma.  Thanks.

> > Line 174 of qpva.m reads
> >            w2 = (xmin+sigma).*(1-rand(n2,1)).^(-1/theta) - sigma;
> > which generates n2 data points drawn from the Pareto II distribution with parameters theta and sigma, with lower-cutoff xmin. Note that xmin=0 corresponds to the "full" distribution with "no" lower cutoff.
Laurent asked

I think I don't get this correctly. I was expecting to get an output in the form of data array, plus theta and sigma, given xmin. However, I don't know what I need exactly to type in MATLAB's command window to generate such output. Please bear with me. I just started to use MATLAB last time I was at SFI. Thanks

Aaron Clauset wrote: Okay, no problem. Put the following into matlab: Then, "x" will be an array of length "n" of values distributed by a Pareto II with parameters xmin = 0, theta = 2, and sigma =3. You can visualize the distribution by the last two commands. (Note by DRW: these commands dont call an external functions, they just generate a function, and the tail will differ in each trial: why differ? Why does the tail vary so much? (a) from the body? (b) from trial to trial?)

xmin = 0;
theta = 2;
sigma = 3;
n = 10^5;
x = (xmin+sigma).*(1-rand(n,1)).^(-1/theta) - sigma;
figure(1);
loglog(sort(x),(n:-1:1)./n);

Or from Cosma's Equation 4 should this be

x = (xmin+(1-rand(n,1)/sigma).^(-theta-1);

-Aaron. Then try this very slow calculation to get the K-S probability p

x = (xmin+sigma).*(1-rand(n,1)).^(-1/theta) - sigma;
[p, gof] = qpva(x, xmin);
figure(2);
loglog(sort(x),(n:-1:1)./n);

[edit] Comparison to the Shalizi MLE by Pareto II

As noted above, in Tsallis's exposition:

The cumulative probability distribution for the q-exponential function in (1) (Tsallis 2004:8) is also a q-exponential function

$P(X) \equiv \int_{X}^{\infin} dx\ p_{e_q^x}(x) = e_{q_M}^{-X/q_M} = [1 - \frac{(1-q)}{(2-q)}x\kappa]^\frac{(2-q)}{(1-q)} \ \mathrm{(equation\ 5)}$ , with

$q_M \equiv 1/(2-q)$ , where for $q \ge 2\ \ p(x)$ is not normalizable. Thus, q cannot be fitted for a range at or above 2. (Note that

κ

does not enter into this equation).

Shalizi took as his cumulative distribution

$P(X) \equiv e_{q_M}^{-X/q_M} = [1 - (1-q_M)x/\lambda]^{1/(1-q_M)}, \lambda=1/\kappa \ \mathrm{(equation\ A)}$ ,

Substituting to obtain Pareto II parameters, $,θ,σ$ , we may take his equations (8) and (10) for the MLE estimation of $\hat{\theta}, \hat{\sigma}$ in terms of $\hat{q}, \hat{\kappa}$

$q_M \equiv \frac{1}{(2-q)} = \frac{\theta+1}{\theta}\ \mathrm{hence} \ q = \frac{\theta+2}{\theta+1} \ \mathrm{and} \ \hat{q} = \frac{\hat{\theta}+2}{\hat{\theta}+1}$

$\hat{\kappa} = \hat{\theta}/\hat{\sigma}$

[edit] Another set of questions

(Probably obsolete)

If the MLE defined above is done successfully, our question is: WILL THIS CONVERGE TO THE Shalizi Pareto II MLE?

We have a matlab program from Aaron Clauset that solves the 2-parameter Shalizi Pareto II MLE formula for $\hat{\theta}$ and $\hat{\sigma}$ (and also finds $x m i n$ for the best fitting q distribution). In Shalizi's MLE, to convert to $κ$ and q to the Pareto II, we have

$κ = θ / σ$ (equation 11) and

$q = 1 + 1 / θ$ (equation 11)

Then

$\hat{\theta} = n \left [ {\sum_{i=1}^n} ln(x_i/\sigma) \right ]^{-1}$ (equation 12)

and $\hat{\sigma}$ is estimated by the transcendental:

$\hat{\sigma} = (\hat{\theta}+1)/n \sum_{i=1}^n x_i/(1+x_i/\hat{\sigma})$ (equation 13)

[edit] References

Tsallis, Constantino. 2004. Nonextensive Statistical Mechanics: Construction and Physical Interpretation. Chapter 1, in Murray Gell-Mann and Constantino Tsallis, Nonextensive Entropy: Interdisciplinary Applications pp. 1-53.

Tsallis, Constantino. 2001. Nonextensive statistical mechanics and thermodynamics: Historical back-ground and present status. In Nonextensive Statistical Mechanics and Its Applications, eds. S. Abe e Y. Okamoto (Series Lecture Notes in Physics, Springer-Verlag, Heidelberg, 2001), pp.~3--98 [ISBN 3-540-41208-5]. The appendix gives properties of q-exponentials and q-logarithms.

Beck, C. (2002) Non-additivity of Tsallis entropies and fluctuations of temperature Europhys. Lett., 57 (3), p. 329

Straub, John E., and Ioan Andricioaei 1999 Computational methods inspired by Tsallis statistics: Monte Carlo and molecular dynamics algorithms for the simulation of classical and quantum systems Brazilian Journal of Physics vol.29 n.1

Some properties of q-logarithm and q-exponential functions in Tsallis statistics Takuya Yamano 2002 Physica A: Statistical Mechanics and its Applications 305(3-4): 486-496.

Tsallis’ entropy maximization procedure revisited

Evaluating results

[edit] Note 1

DRW: My previous equation 9 was incorrect. I had set $n / η$ =1. Natasa commented as follows:

To find the MLE (maximum $\mathfrak{k}(\eta)$ , take the first derivative of the log-likelihood $\ \mathfrak{k}(\eta)$ with respect to the parameter $η$ and set it equal to zero.

(Comment by NK: from (equation 8) you get the following equation for the derivative that has to equate to 0:

$d \mathfrak{k}/d \eta = -\sum_{i=1}^{n} log (\ 1 + x_i/\eta) + \sum_{i=1}^n (\frac{\eta+1}{\eta(\eta+x_i)}) = 0$

(I don't know how to nicely (analytically) get

η

out of that. I think the easiest way would be to use a kind of iterative algorithm to get to the solution (so something similar to that of the MLE estimation, that Shalizi did). --NK

[edit] Note 2: Eliminated - Curve fitting: estimating q with the q-logarithmic function

To estimate q it is sufficient to fit a straight line for the THE CUMULATIVE Z(k) against k, where in this illustration q ~ 1.7. This can be done by maximizing $r 2$ for Z(k) ~ k (letter k not Greek $κ$ ), where $κ$ is just the slope ( -slope more precisely) for that value of q. We can now compare maximizing $r 2$ with the MLE solution above, and derive $κ$ as the (-slope) in the figure.

Retrieved from "http://intersci.ss.uci.edu/wiki/index.php/Estimating_Tsallis_q"

Categories: Matlab software | R software | Entropy | Q exponential

Estimating Tsallis q

From InterSciWiki

Contents

[edit] q exponential and q logarithm

[edit] Probability distribution and Cumulative probability distribution (CDF)

[edit] Comparison to the Shalizi MLE by Pareto II

[edit] Possibility of a more direct solution

[edit] Deriving the MLE

[edit] Testing the MLE 1

[edit] x not cumulative

[edit] straight line?

[edit] try new MLE

[edit] Testing the MLE 2 ???

[edit] Comparison to the Shalizi MLE by Pareto II

[edit] Another set of questions

[edit] References

[edit] See also

[edit] Note 1

[edit] Note 2: Eliminated - Curve fitting: estimating q with the q-logarithmic function

Views

Personal tools

Navigation

Search

Toolbox